BUG: second derivative of com_to_cdm_function diverges at the burn start time

[thc1006]

Describe the bug

Flight.u_dot_generalized builds the second time derivative of the CM offset like this:

r_CM_ddot = Vector([0, 0, r_CM_z.differentiate(t, order=2)])

com_to_cdm_function has a kink in its first derivative at the burn start time: flat before ignition, sloped after. When its source is a callable, Function.differentiate(x, dx=1e-6, order=2) falls back to (f(x+dx) - 2 f(x) + f(x-dx)) / dx**2, and across a kink of slope s that evaluates to about s / dx. It grows without bound as the step shrinks.

Whether it bites depends on how com_to_cdm_function ends up being sourced, which in turn depends on the motor:

fixture	source	r_CM'' at burn_start_time	- m * r_CM''
calisto	ndarray, linear	0	0 N
calisto_robust	ndarray, linear	0	0 N
calisto_liquid_modded	ndarray, linear	0	0 N
calisto_hybrid_modded	callable	-3.749668e+04 m/s²	+1.92e+06 N

That last row is a 51.1 kg rocket, so the spurious force is about 3800 times its weight.

The array cases are not really fine either. For a linearly interpolated array the evaluator's second derivative is identically zero, so the - total_mass * r_CM_ddot term from the documented equations is silently dropped everywhere rather than being wrong at one point.

To Reproduce

Against the fixtures:

def test_probe(calisto_hybrid_modded):
    rocket = calisto_hybrid_modded
    t0 = rocket.motor.burn_start_time
    print(rocket.com_to_cdm_function.differentiate(t0, order=2))
    # -37496.68...

And without any fixtures, showing the 1/dx growth on a bare kink:

from rocketpy.mathutils.function import Function

slope = -3.1134e-4
kink = Function(lambda t: 0.0 if t <= 1.01 else slope * (t - 1.01))
for dx in (1e-5, 1e-6, 1e-7):
    print(dx, kink.differentiate(1.01, dx=dx, order=2))

# 1e-05 -31.13...
# 1e-06 -311.3...
# 1e-07 -3113....

Expected behavior

r_CM'' stays bounded and near the analytic value on both sides of ignition, and the term contributes the same way whatever the motor happens to produce for com_to_cdm_function.

Additional context

How much this changes a flight today is modest, and I want to be straight about that. A rocket on the rail does not reach u_dot_generalized at the burn start time, since udot_rail1 covers the rail phase. Running calisto_hybrid_modded from example_plain_env with rail_length=5.2, suppressing this term alone moves apogee by 0.90 m out of 4405, about 0.02 percent.

It matters more for a flight that is integrated through the ignition instant while already off the rail, or one that starts there. Since #1085 the solver is also asked to stop exactly at burn_start_time and burn_out_time whenever burn_start > 0, which is a good fix for #411 and does mean the right hand side gets sampled precisely on the time where this term is not defined.

On what to do about it, I do not have a strong opinion and it seems like your call. The project has already solved the same shape of problem once: total_mass_ddot comes from total_mass_flow_rate.differentiate_complex_step(t) rather than from a second derivative of total_mass. The equivalent here would be to expose the first derivative of com_to_cdm_function as its own Function and differentiate that once, or to assemble r_CM' and r_CM'' from the propellant mass and propellant center of mass and their derivatives. Either way the 1/dx factor goes away.

Separately, and possibly worth its own look: Function.differentiate(order=2) returns zero for an array source with linear interpolation and s/dx for a callable, for the same underlying function. Making those two agree would remove the surprise here regardless of which fix lands.

I ran into this while looking at something else in u_dot_generalized, which I have filed as a separate issue since the two are independent.

[thc1006]

/cc @zuorenchen