diff --git a/NEWS.md b/NEWS.md
index 13958fb7..fc4d01ea 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,8 @@
+## VortexStepMethod v2.2.0 2025-08-28
+### Added
+- The kwarg `aero_coeffs` to the function `linearize`: if true the linearization will output
+    normalized coefficients instead of moments and forces.
+
 ## VortexStepMethod v2.1.0 2025-08-11
 ### Changed
 #### 1. Core New Functionality: YAML Geometry Support
diff --git a/src/solver.jl b/src/solver.jl
index f3ea1085..ecf34e3b 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -694,6 +694,7 @@ function linearize(solver::Solver, body_aero::BodyAerodynamics, y::Vector{T};
         delta_idxs=nothing,
         va_idxs=nothing,
         omega_idxs=nothing,
+        aero_coeffs=false,
         kwargs...) where T
 
     !(length(body_aero.wings) == 1) && throw(ArgumentError("Linearization only works for a body_aero with one wing"))
@@ -745,9 +746,15 @@ function linearize(solver::Solver, body_aero::BodyAerodynamics, y::Vector{T};
         end
 
         solve!(solver, body_aero; kwargs...)
-        results[1:3] .= solver.sol.force
-        results[4:6] .= solver.sol.moment
-        results[7:end] .= solver.sol.group_moment_dist
+        if !aero_coeffs
+            results[1:3] .= solver.sol.force
+            results[4:6] .= solver.sol.moment
+            results[7:end] .= solver.sol.group_moment_dist
+        else
+            results[1:3] .= solver.sol.force_coeffs
+            results[4:6] .= solver.sol.moment_coeffs
+            results[7:end] .= solver.sol.group_moment_coeff_dist
+        end
         return nothing
     end
 
diff --git a/test/body_aerodynamics/test_results.jl b/test/body_aerodynamics/test_results.jl
index 0ac6d232..2b47dc5b 100644
--- a/test/body_aerodynamics/test_results.jl
+++ b/test/body_aerodynamics/test_results.jl
@@ -70,11 +70,25 @@ end
         delta_idxs=11:14,
         moment_frac=0.1
     )
+
+    coeff_jac, coeff_lin_res = VortexStepMethod.linearize(
+        solver, 
+        body_aero, 
+        base_inputs; 
+        theta_idxs=1:4, 
+        va_idxs=5:7, 
+        omega_idxs=8:10,
+        delta_idxs=11:14,
+        moment_frac=0.1,
+        aero_coeffs=true
+    )
     
     # Verify that linearization results match nonlinear results at operating point
     baseline_res = VortexStepMethod.solve!(solver, body_aero; log=false)
     baseline_res = [solver.sol.force; solver.sol.moment; solver.sol.group_moment_dist]
+    coeff_baseline_res = [solver.sol.force_coeffs; solver.sol.moment_coeffs; solver.sol.group_moment_coeff_dist]
     @test baseline_res ≈ lin_res
+    @test coeff_baseline_res ≈ coeff_lin_res
     
     # Define test cases
     test_cases = [
@@ -88,6 +102,7 @@ end
     for (input_name, indices, magnitudes) in test_cases
         @testset "Input: $input_name" begin
             max_error_ratio = 0.0
+            coeff_max_error_ratio = 0.0
             
             # Select one index at a time for focused testing
             for idx in indices
@@ -128,18 +143,25 @@ end
                     # Get nonlinear solution
                     nonlin_res = VortexStepMethod.solve!(solver, body_aero, nothing; log=false)
                     nonlin_res = [solver.sol.force; solver.sol.moment; solver.sol.group_moment_dist]
+                    coeff_nonlin_res = [solver.sol.force_coeffs; solver.sol.moment_coeffs; solver.sol.group_moment_coeff_dist]
                     @test nonlin_res ≉ baseline_res
+                    @test coeff_nonlin_res ≉ baseline_res
                     
                     # Compute linearized prediction
                     lin_prediction = lin_res + jac * perturbation
+                    coeff_lin_prediction = coeff_lin_res + coeff_jac * perturbation
                     
                     # Calculate error ratio for this case
                     prediction_error = norm(lin_prediction - nonlin_res)
                     baseline_difference = norm(baseline_res - nonlin_res)
                     error_ratio = prediction_error / baseline_difference
+                    coeff_prediction_error = norm(coeff_lin_prediction - coeff_nonlin_res)
+                    coeff_baseline_difference = norm(coeff_baseline_res - coeff_nonlin_res)
+                    coeff_error_ratio = coeff_prediction_error / coeff_baseline_difference
 
                     # Update max error ratio
                     max_error_ratio = max(max_error_ratio, error_ratio)
+                    coeff_max_error_ratio = max(coeff_max_error_ratio, coeff_error_ratio)
                     
                     # For small perturbations, test that error ratio is small
                     if idx == first(indices) && mag == first(magnitudes)
@@ -149,6 +171,7 @@ end
             end
             
             @info "Max error ratio for $input_name: $max_error_ratio"
+            @info "Max coeff error ratio for $input_name: $coeff_max_error_ratio"
         end
     end
     
@@ -275,4 +298,4 @@ end
             end
         end
     end
-end
\ No newline at end of file
+end