Hi everyone,
while investigating the EME2000/TOD propagation discrepancy with inertial forces, I found a separate frame-dependent behaviour in NRLMSISE00.
Problem
The same physical spacecraft position produces different NRLMSISE00 densities depending on whether it is supplied in EME2000 or TOD.
As you can see, both density evaluations use the same NRLMSISE00 instance and epoch:
Frame eme2000 = FramesFactory.getEME2000();
Frame tod = FramesFactory.getTOD(IERSConventions.IERS_2010, true);
Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
// Same physical position expressed in two frames
Vector3D positionTod =
eme2000.getStaticTransformTo(tod, date)
.transformPosition(positionEme);
double rhoEme = atmosphere.getDensity(date, positionEme, eme2000);
double rhoTod = atmosphere.getDensity(date, positionTod, tod);
// Verification that both positions represent the same ITRF point
Vector3D itrfFromEme =
eme2000.getStaticTransformTo(itrf, date)
.transformPosition(positionEme);
Vector3D itrfFromTod =
tod.getStaticTransformTo(itrf, date)
.transformPosition(positionTod);
For one LEO state at 2027-01-01T00:00:00 UTC, I obtain:
|itrfFromEme - itrfFromTod| = 7.3e-9 m
rhoEme = 1.071411857531e-13 kg/m³
rhoTod = 1.071681224500e-13 kg/m³
The positions therefore agree in ITRF at nanometre level, while the relative density difference is approximately 2.5e-4.
At other epochs during a seven-day propagation, the relative density difference reached approximately 2.3e-3.
I am using Orekit 13.1 with IERS 2010 frames.
Potential source
In NRLMSISE00.getDensity(...), the geodetic coordinates are correctly calculated using the Earth body shape:
GeodeticPoint inBody = earth.transform(position, frame, date);
However, local solar time is calculated using the spacecraft and Sun coordinates directly in the caller-provided frame:
Vector3D sunPos = getSunPosition(date, frame);
double lst = FastMath.PI + FastMath.atan2(
sunPos.getX() * position.getY() -
sunPos.getY() * position.getX(),
sunPos.getX() * position.getX() +
sunPos.getY() * position.getY());
This calculation uses the XY plane of the supplied frame. EME2000 and TOD have different pole directions due to precession and nutation, so the projected angle is not invariant under the EME2000/TOD rotation.
This seems to explain why the same ITRF position produces different densities.
Propagation impact
In my EME2000/TOD propagation comparison, before correcting the atmosphere evaluation, the seven-day LEO position differences were approximately:
Without drag: 0.21 m
With drag: 15.27 m
The direct drag-acceleration mismatch reached approximately:
4.8e-10 m/s²
The physical force stack and initial physical state were otherwise equivalent in both frames.
Possible solution
As a test, I wrapped the atmosphere and always evaluated density in its body frame:
public double getDensity(AbsoluteDate date,
Vector3D position,
Frame frame) {
Frame bodyFrame = delegate.getFrame();
Vector3D bodyPosition =
frame.getStaticTransformTo(bodyFrame, date)
.transformPosition(position);
return delegate.getDensity(date, bodyPosition, bodyFrame);
}
This makes NRLMSISE00 evaluate the geodetic position, spacecraft direction and Sun direction consistently in ITRF.
After this change:
rhoEme = rhoTod // to displayed precision
drag acceleration mismatch ≈ 1e-20 m/s²
With converged integration settings (DormandPrince853 with 30 s max step and 1e-9 m position tolerance) the seven-day LEO comparison became:
Without drag: 0.0412 m
With drag: 0.0554 m
A possible internal correction in NRLMSISE00 would be to calculate local solar time after transforming both the spacecraft and Sun positions to earth.getBodyFrame(), or equivalently to derive it from the Earth-fixed longitude and subsolar longitude.
Happy to provide more details if needed.
Thank you,
Daniele
