Potential Issue in JB2008 Implementation with Local Hour Angle Computation
I was going through the JB2008
Orekit implementation and found a potential discrepancy between the original Fortran code from Bowman and Orekit's implementation. The original code takes in the following inputs:
C AMJD : Date and Time, in modified Julian Days
C and Fraction (MJD = JD-2400000.5)
C SUN(1) : Right Ascension of Sun (radians)
C SUN(2) : Declination of Sun (radians)
C SAT(1) : Right Ascension of Position (radians)
C SAT(2) : Geocentric Latitude of Position (radians)
C SAT(3) : Height of Position (km)
C F10 : 10.7-cm Solar Flux (1.0E-22*Watt/(M**2*Hertz))
C (Tabular time 1.0 day earlier)
C F10B : 10.7-cm Solar Flux, ave.
C 81-day centered on the input time
C (Tabular time 1.0 day earlier)
C S10 : EUV index (26-34 nm) scaled to F10
C (Tabular time 1.0 day earlier)
C S10B : EUV 81-day ave. centered index
C (Tabular time 1.0 day earlier)
C XM10 : MG2 index scaled to F10
C (Tabular time 2.0 days earlier)
C XM10B : MG2 81-day ave. centered index
C (Tabular time 2.0 days earlier)
C Y10 : Solar X-Ray & Lya index scaled to F10
C (Tabular time 5.0 days earlier)
C Y10B : Solar X-Ray & Lya 81-day ave. centered index
C (Tabular time 5.0 days earlier)
C DSTDTC : Temperature change computed from Dst index
In Orekit, the Java translation of this main JB2008
Fortran subroutine is available via this method.
For this issue, what I am concerned with are the SAT
and SUN
arrays. In particular, SAT(1)
or as is called in Orekit satLon
, which is supposed to represent the right ascension of the satellite, as well as SUN(1)
or sunRA
, which is the right ascension of the Sun. This is important because internally in the JB2008 code it is used to compute the LHA (local hour angle) that is then an input to Equation (16) for τ from Jacchia's original work, which is simply:
LHA = SAT(1) - SUN(1)
In Orekit, this occurs here for the scalar implementation.
Note that LHA as an input is not a JB2008-specific input, but it's part of all Jacchia models. In fact, this shows up in Appendix B of Vallado as well, where he suggests an alternative way to compute it. Separately from LHA (which should be computed using TOD / inertial coordinates as it represents the angle in the celestial equator), the Jacchia models also require the geocentric/geodetic latitude and height inputs. This got me at first when going through the JB2008.for
code as SAT
is an array with right ascension (1) and then geodetic latitude (2) and height (3), which is definitely not a common coordinate system combination.
My concern is that when using the generic Orekit Atmosphere
interface and its getDensity
method, the inputs that are passed into the original interface method are:
- Equivalent of
SUN
:sunInBody.getLongitude()
,sunInBody.getLatitude()
- Equivalent of
SAT
:inBody.getLongitude()
,inBody.getLatitude()
,inBody.getAltitude()
where sunInBody = earth.transform(sunPos, ecef, date)
and inBody = earth.transform(position, frame, date)
are GeodeticPoint
s as can be seen here.
Therefore, to me it looks like the generic Orekit getDensity
implementation is passing in the wrong values to the internal original getDensity
interface. It appears to be passing in the geodetic longitude instead of the inertial right ascension, resulting in a different LHA.
I may have missed something as I was going through the code, so please feel free to close this issue in that case, but I wanted to get another set of eyes on this. Thank you!