[POSSIBLE BUG]: Mismatch in airspeed velocity when reconstructing from rotation matrix derivative

[minervino96]

I am copying and pasting a post I started in the forum a few months ago. It is of low priority, but still worth keeping on the radar in case there is an actual bug.

While comparing the airspeed velocity as computed internally in Tudat with a manual reconstruction using rotation matrices and angular velocity, I noticed a small discrepancy that I would like to better understand.

This is all done during propagation, with spacecraft and Earth retrieved from the SystemOfBodies.

Given:

R_eci2ecef = Earth.inertial_to_body_fixed_frame
dR_eci2ecef = Earth.inertial_to_body_fixed_frame_derivative

I reconstruct the airspeed velocity in the central body–fixed (ECEF) frame as:

airspeed_velocity = np.dot(R_eci2ecef, spacecraft.velocity) + np.dot(dR_eci2ecef, spacecraft.position)

This gives:

(-35.05804206702926, 1765.1279468551625, 7246.924546200871)

This almost matches the internally computed airspeed velocity in:

spacecraft.flight_conditions.airspeed_velocity

which returns:

(-35.05804206702939, 1765.1279468551625, 7246.924546200871)

The difference appears only in the last two digits of the x-direction, likely due to floating point precision or minor implementation differences between Python and C++.

---

Now, when trying to reconstruct the airspeed velocity in the inertial frame, I do:

R_ecef2eci = Earth.body_fixed_to_inertial_frame
airspeed_velocity_inertial = spacecraft.velocity + np.dot(R_ecef2eci, np.dot(dR_eci2ecef, spacecraft.position))

which gives:

(-1718.3317034143038, -323.5535831347614, 7251.032364380492)

However, using the inertial angular velocity vector:

omega = Earth.inertial_angular_velocity
airspeed_velocity_inertial = spacecraft.velocity - np.cross(omega, spacecraft.position)

yields:

(-1718.3317137228362, -323.553528387783, 7251.0323643767315)

which is again very close, but the results differ from the 5th decimal digit onward.

---

I expected that:

np.dot(R_ecef2eci, dR_eci2ecef)

would yield the skew-symmetric matrix of omega, but this is not quite the case:

[[7.550516270495057e-12, 7.292094583056997e-05 ,-2.504434634960097e-09]
 [-7.29209458305712e-05, 7.550558840799778e-12, 1.7315942385909157e-07]
 [2.5043987757688435e-09, -1.731594243777255e-07, 4.258485558200139e-17]]

while the angular velocity vector is:

[1.7315942385909167e-07, 2.5043987757688294e-09, 7.292094583057e-05]

---

So, while the match is very close, I wonder whether the difference should be attributed to numerical differentiation artefacts, or is there something conceptually different between the two approaches?

Thanks in advance!