BUG: correct the nozzle gyration tensor parallel axis term

[thc1006]

This patch corrects the parallel axis term in Rocket.evaluate_nozzle_gyration_tensor. It is currently a quarter of what the exit plane integral gives. The unit stated in the docstring is fixed too.

Pull request type

• Code changes (bugfix, features)

Checklist

• Tests for the changes have been added (if needed)
• Docs have been reviewed and added / updated
• Lint has passed locally
• All tests (pytest tests -m slow --runslow) have passed locally

Some notes so the boxes are not overstated.

The lint line still mentions black, so to be precise about what I ran: ruff check ., ruff format --check . and pylint rocketpy/ tests/ docs/, matching the Linters workflow. The only docs change is prose inside an existing numpydoc docstring, with no directives or references, so I did not build the docs locally.

I left the last box unticked rather than claim something I did not see finish. Of the 35 slow tests, 34 passed. test_hrrr_atmosphere failed once on a forecast download and passes on re-run, on this branch and on develop alike; its own docstring notes that the HRRR latest model sometimes lacks 24 hours of forecast. test_monte_carlo_simulate[True], the parallel case, does not finish inside 40 minutes on this machine, so I cannot report it either way. The sequential case of the same test passed, and the two run the same simulations.

Current behavior

evaluate_nozzle_gyration_tensor builds the tensor as:

S_noz_33 = 0.5 * self.motor.nozzle_radius**2
S_noz_11 = S_noz_22 = 0.5 * S_noz_33 + 0.25 * self.nozzle_to_cdm**2

S_noz is the second moment of the nozzle exit plane per unit area, taken about the center of dry mass. For a uniform disk of radius R centred on the body axis at a distance d from that point:

S = (1/A) * integral over A of (|x|^2 * I - x x^T) dA
  = diag(R^2/4 + d^2,  R^2/4 + d^2,  R^2/2)

The axial component already matches. The two lateral ones should carry d^2 rather than d^2/4.

The quickest way to see it is the point limit. Let R go to zero and the exit plane collapses to a single mass flux location at distance d, so S has to tend to diag(d^2, d^2, 0). Today it returns a quarter of that, and nothing about that limit depends on the rest of the model.

The older solid propulsion path already uses the full squared distance for the same quantity:

# Flight.u_dot
- mass_flow_rate_at_t * ((nozzle_radius / 2) ** 2 + (c - b * mu / rocket_dry_mass) ** 2)

One smaller thing in the same function: the docstring says the tensor is in kg*m². Since T05 = mdot * S - I_dot has to be kg·m²/s and mdot is kg/s, S is in m². I reworded that line and the tolerance comment in the existing test, which read 1e-3 * 1e-2 * 1e-2 # Equivalent to 1g * 1cm^2.

New behavior

The parallel axis term becomes self.nozzle_to_cdm**2. On the calisto fixture S_11 goes from 0.3940207 to 1.5752660 m².

S_noz reaches the trajectory only through T05 @ w in u_dot_generalized, so it is scaled by mass flow rate and body rate. How much of it shows up depends on the rocket. On flight_calisto_custom_wind, a solid motor burning for 3.9 s, it is small. At apogee:

	before	after
aerodynamic moment M1	-0.509418 N·m	-0.519917 N·m
velocity x	-14.826366 m/s	-14.593411 m/s
apogee above ground	3423.2700 m	3423.8760 m

On the Defiance example it is not small. That is a hybrid burning for 6.5 s, and the impact point moves from (1625.57, 81.79) to (1609.40, 87.03), roughly 17 m over a 1626 m range. That takes test_defiance_rocket_matches_reference_flight_metrics out of its band, and I re-recorded the two impact constants. There is more on that in the comments, including the check against the measured apogee, which the change moves closer rather than further.

Three recorded values in tests/unit/simulation/test_flight.py move past their 5e-3 tolerance and are updated. I re-recorded each tuple from the same run so it stays internally consistent, which is why the roll moment reads -1.005368e-18 where it used to read -2.089120e-14. Both are numerical zero, and that component is unchanged in substance. test_accelerations stays inside tolerance and is left alone. Everything else in tests/unit and tests/integration gives the same result it does on develop.

Breaking change

• No

No API change. Simulation results do shift slightly, as above.

Additional information

Two tests are added next to the existing one in tests/unit/rocket/test_rocket.py:

• test_evaluate_nozzle_gyration_tensor_matches_the_exit_disk is parameterized over calisto, calisto_liquid_modded and calisto_hybrid_modded, and checks all three diagonal entries against R^2/4 + d^2 and R^2/2 built from each rocket's own geometry, plus the off-diagonal zeros. It does not restate the recorded number.
• test_evaluate_nozzle_gyration_tensor_parallel_axis_term subtracts the disk contribution and checks what is left equals nozzle_to_cdm**2, which is the coefficient this patch changes.

I put the one line back to 0.25 * ... to check the tests really hold it. All five nozzle tests fail, including the parameterized ones on the liquid and hybrid motors, and the other 61 tests in that file still pass.

While reading u_dot_generalized I ran into two other things that look independent of this change, so I filed them as #1185 and #1186 rather than fold them in here.

[thc1006]

/cc @zuorenchen

[thc1006]

CI caught something I missed locally. I ran tests/unit and tests/integration but not tests/acceptance, and test_defiance_rocket_matches_reference_flight_metrics goes red.

Here is what this change does to that example:

	develop	this PR
max_speed	444.2545	444.2364
max_acceleration	10400.7640	10400.7640
x_impact	1625.5656	1609.4019
y_impact	81.7854	87.0258

y_impact lands outside its 3 m band. x_impact stayed inside, but only because pytest.approx takes the larger of rel and abs, so the band there works out at 16.26 m against a move of 16.15 m. A 0.1 m margin is not something I want to leave behind on a matrix that includes Windows and macOS.

The four REFERENCE_* constants agree with what develop produces to five significant figures, and the test's docstring says it guards the deterministic example's metrics, so I read them as recorded from the model rather than measured. The measured quantity in that file is MEASURED_APOGEE_AGL = 9308.32, and against that one the change goes the right way:

	apogee AGL	error vs measured
develop	9242.3701 m	0.709 %
this PR	9243.8350 m	0.693 %

The other 17 acceptance tests pass untouched, including the measured-flight comparisons for Bella Lui, NDRT 2020 and Prometheus.

I pushed 2953fd4 re-recording the two impact constants. Re-recording a guard someone else set is your call rather than mine, so if you would rather keep them pinned and take this a different way, say so and I will drop that commit.

While I have the numbers: on this rocket the correction is worth about 17 m of impact point, not the sub-metre it is worth on flight_calisto_custom_wind. Defiance is a hybrid burning for 6.5 s, so the mass flow rate stays large for longer and T05 @ w has more to act on. My first description understated that and I have fixed it.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 89.96%. Comparing base (7e785a6) to head (2953fd4).

Additional details and impacted files

@@           Coverage Diff            @@
##           develop    #1188   +/-   ##
========================================
  Coverage    89.96%   89.96%           
========================================
  Files          131      131           
  Lines        17527    17527           
========================================
  Hits         15769    15769           
  Misses        1758     1758           

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[thc1006]

Closing the one loose end from the description. I said test_monte_carlo_simulate[True] did not finish inside 40 minutes here and that I could not report it either way. I have since run the same test on unmodified develop on this machine and it behaves identically: the sequential case passes, the parallel one is still going at the 40 minute mark. So it is this machine rather than anything in the patch.