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
5

Hello Romain, thanks for the suggestion. I gave it a naive try but I’m not quite there yet

fieldOrbit = FieldOrbit(FramesFactory.getGCRF(),date, Constants.WGS84_EARTH_MU)

date is an instance of <class 'org.orekit.time.AbsoluteDate'>. However, the following error pops up

NotImplementedError: ('instantiating java class', <class 'org.orekit.orbits.FieldOrbit'>)

I’m using Orekit through its Python bindings. Could you confirm that a FieldOrbit instance is necessary in order to create a FieldSpacecraftState object that can be passed to the force model so as to retrieve the correct acceleration method ?

That is, the method that return a FieldVector3D object for which a proper gradient method is available.

public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration​(FieldSpacecraftState<T> s, T[] parameters)
6

Hi,

FieldOrbit is an abstract class, so you can’t instantiate it as is. You can use FieldCartesianOrbit for example.
You can look at this class in Java for inspiration, in particular method buildBasicGradientSpacecraftState:

Cheers,
Romain.

7

Thanks Romain, your script was indeed of great help. I got to the point where I have a FieldSpacecraftState object instantiated

(Pdb) fieldSpacecraftState
<FieldSpacecraftState: FieldSpacecraftState{orbit=Cartesian parameters: {P(0.0, 4509927.0504078, 4509927.0504078), V(0.0, 0.0, 0.0)}, attitude=org.orekit.attitudes.FieldAttitude@2de50ee4, mass=org.hipparchus.analysis.differentiation.Gradient@c3e3eef0, additional={}, additionalDot={}}>

The velocity was zeroed out at initialization which was expected.

Going back to your initial comment

One way of getting the force’s Jacobian is to use the Field version of acceleration with Gradient.

I thought calling acceleration with fieldSpacecraftState as opposed to a plain SpacecraftState object would do the trick and call the Field version of the acceleration method. Unfortunately, it doesn’t return a differentiable acceleration, but instead throws an invalid argument error

(Pdb)         nonCentralAcceleration = forceModel.acceleration(fieldSpacecraftState,[Constants.WGS84_EARTH_MU]) 
*** orekit.InvalidArgsError: (<class 'org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel'>, 'acceleration', (<FieldSpacecraftState: FieldSpacecraftState{orbit=Cartesian parameters: {P(0.0, 4509927.0504078, 4509927.0504078), V(0.0, 0.0, 0.0)}, attitude=org.orekit.attitudes.FieldAttitude@2de50ee4, mass=org.hipparchus.analysis.differentiation.Gradient@c3e3eef0, additional={}, additionalDot={}}>, [398600441800000.0]))

This may be a complete misunderstanding on my end of how Orekit exposes force field partial derivatives

8

i think the problem is that your array of parameters is of float instead of Gradient too

9

Thanks, that worked, forceModel.acceleration does return a <FieldVector3D> object. I’m seeing its getX method returns an object of type <CalculusFieldElement: org.hipparchus.analysis.differentiation.Gradient@d4dc9f2b> which is promising, but I’m having a hard time relating it to what I am after, namely the partial derivative of the acceleration vector relative to the cartesian coordinates.

I might be close, but I would be very grateful if you could either complete the above or point me in the right direction

10

If you do getGradient after getX you’ll get what you want for the first component

11

Hello Romain,

Unfortunately, the getX() method of the <FieldVector3D> class returns an object of type <CalculusFieldElement: org.hipparchus.analysis.differentiation.Gradient@a70e97da> that doesn’t appear to have a getGradient method.


print(type(nonCentralAcceleration))
# <class 'org.hipparchus.geometry.euclidean.threed.Vector3D'>

Xcomp = nonCentralAccelerationGradient.getX()
print(dir(X_comp))

'''
['__class__', '__delattr__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', 
'__hash__', '__init__', '__init_subclass__', 
'__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', 
'__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '_jobject', 'abs', 'acos', 'acosh', 'add', 'asin', 'asinh', 'atan', 'atan2', 'atanh', 'boxfn_', 'cast_', 'cbrt', 'ceil',
 'class', 'class_', 'copySign', 'cos', 'cosh', 'divide', 'equals', 'exp', 'expm1', 'exponent', 
'field', 'finite', 'floor', 'getClass', 'getExponent', 'getField', 'getPi', 'getReal', 'hashCode', 'hypot', 'infinite', 'instance_', 'isFinite', 'isInfinite', 'isNaN', 'isZero', 'linearCombination', 
'log', 'log10', 'log1p', 'multiply', 'naN', 'negate',
 'newInstance', 'norm', 'notify', 'notifyAll', 'of_', 'parameters_', 'pi', 'pow', 'real', 'reciprocal', 'remainder', 'rint', 'rootN', 'round', 
'scalb', 'sign', 'sin', 'sinCos', 'sinh', 'sinhCosh',
 'sqrt', 'square', 'subtract', 'tan', 'tanh', 'toDegrees', 'toRadians', 'toString', 'ulp', 'wait', 'wrapfn_', 'zero']
'''

Your answer is in line with the documentation of Hipparcus ( a getGradient method should indeed be present), so something is amiss here …

I’m running Orekit 12.2. I can move the conversation to the dedicated subforum instead of continuing the off-topic conversation here.

12

Which Python wrapper are you using? JCC or jpype?

13

My guess would be JCC

import orekit
print(orekit.JCC_VERSION)
# '3.14'
14

Ok so try doing Gradient.cast_(x) before calling getGradient

15

The cast did the trick :partying_face: Thank you !

16

FYI I’ve added on the develop branch a DerivativeStateUtils with functions to build Gradient FieldOrbit and FieldAbsolutePVCoordinates from their standard counterpart

1 Like