Holmes Featherstone model restricted to J2 returns asymetrical gravity acceleration

Orekit usage
1

Hello there,

This might be a misunderstanding on my end but I’ve run into strange results upon using the Holmes Featherstone attraction model. Although parameterized to only retain terms up to degree 2, querying the model at a location positioned on the GCRF YZ plane returns an acceleration that features a small, yet unexpected non-zero X component

Below is a code excerpt demonstrating how the model is instantiated in Python and queried. With position_N = np.array([ 0. , 4509927.0504078, 4509927.0504078]) and gravityShDegree = 2 ,

  • nonCentralAcceleration amounts to array([-3.75983971e-05, 1.68947930e-02, -5.54440375e-03])
  • I was expecting it to lie within the YZ plane, and evaluate closer to array([-0. , 0.01687864, -0.00562621])

The coefficient files that is in use appears to be eigen-6s.gfc .

Can anybody shed some light on this ?

    def getEarthGravityAcceleration(self,date,position_N,gravityShDegree):
        '''
        Return Earth gravity acceleration using Orekit's spherical harmonics model
        @args date (orekit.AbsoluteData) date
        @args position_N (np.array) position [m] in GCRF
        @args gravityShDegree (int) maximum degree and order of spherical harmonics expansion
        @return acceleration_N gravity acceleration in GCRF
        '''
        from org.orekit.frames import FramesFactory
        from org.orekit.forces.gravity.potential import GravityFieldFactory
        from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
        from org.hipparchus.geometry.euclidean.threed import Vector3D
        from org.orekit.utils import IERSConventions,PVCoordinates,AbsolutePVCoordinates
        from org.orekit.propagation import SpacecraftState
        from org.orekit.utils import Constants

        gravityProvider = GravityFieldFactory.getNormalizedProvider(gravityShDegree, gravityShDegree)
        forceModel = HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, True), gravityProvider)


        position_N_orekit  = Vector3D(float(position_N[0]),float(position_N[1]),float(position_N[2]))
        velocity_N_orekit  = Vector3D(float(0.0),float(0.0),float(0.0))
        pvCoordinates = PVCoordinates(position_N_orekit, velocity_N_orekit)

        absolutePvCoordinates = AbsolutePVCoordinates(FramesFactory.getGCRF(), date, pvCoordinates)
        spacecraftState = SpacecraftState(absolutePvCoordinates)

        nonCentralAcceleration = forceModel.acceleration(spacecraftState,[Constants.WGS84_EARTH_MU]) 

        nonCentralAcceleration = np.array([
            nonCentralAcceleration.getX(),
            nonCentralAcceleration.getY(),
            nonCentralAcceleration.getZ()])

        return nonCentralAcceleration - Constants.WGS84_EARTH_MU * position_N / np.linalg.norm(position_N) ** 3
2

Maybe this is due to ITRF and GCRF not having the same Z axis? They are close, of course, but ITRF is affected by precession and nutation, so the Z axis wanders a little.

3

Thanks Luc for the quick reply. I had a side question related to the Holmes Featherstone model : I can’t seem to find the gravity gradient matrix (corresponding to the selected SH degree & order) as part of the exposed methods. There is a gradient method that appears to return a vector instead of the expected 3x3 matrix.

4

Hi there,

the method you’re referring to is the gradient of the potential, not the force.
One way of getting the force’s Jacobian is to use the Field version of acceleration with Gradient.

On another note: if you’re only using J2, you may be interested in the class J2OnlyPerturbation.

Cheers,
Romain.

1 Like