CustomAtmosphere interface with Derivatives for STM Harvester

Orekit Python Wrapper usage
1

I am attempting to set up a custom atmosphere density model provider in conjunction with State-Transition-Matrix based covariance propagation. This has been mostly successful so far, but I have hit a snag with the getVelocity() method in the PythonAtmosphere interface. The implementation of this method is as follows:

def getVelocity(
        self, date: AbsoluteDate, position: Vector3D, frame: Frame
    ) -> Vector3D:
        # fixedToFrame: FieldKinematicTransform = self.getFrame().getKinematicTransformTo(
        #     frame, date
        # )
        # posInFixed: FieldVector3D = fixedToFrame.getStaticInverse().transformPosition(
        #     position
        # )
        
        # pvBody: FieldPVCoordinates = FieldPVCoordinates(posInFixed, Vector3D.ZERO)
        # pvFrame: FieldPVCoordinates = fixedToFrame.transformOnlyPV(pvBody)
        # velocity = pvFrame.getVelocity()
        # return velocity
        FieldFlag = False
        if isinstance(date,FieldAbsoluteDate):
            date = FieldAbsoluteDate.toAbsoluteDate(date)
            FieldFlag = True
        if isinstance(position,FieldVector3D):
            position = FieldVector3D.toVector3D(position)
            FieldFlag = True
            
        fixedToFrame: KinematicTransform = self.getFrame().getKinematicTransformTo(
            frame, date
        )
        posInFixed: Vector3D = fixedToFrame.getStaticInverse().transformPosition(
            position
        )
        pvBody: PVCoordinates = PVCoordinates(posInFixed, Vector3D.ZERO)
        pvFrame: PVCoordinates = fixedToFrame.transformOnlyPV(pvBody)
        if not FieldFlag:
            return pvFrame.getVelocity()
        else:
            FieldpvFrame: FieldPVCoordinates = pvFrame.toDerivativeStructurePV(1)
            velocity = FieldpvFrame.getVelocity()
            return velocity

Note the default (non-derivative) interface works perfectly without the STM call (i.e. FieldFlag is false). My issues arise when I try to start providing (zero) derivatives to support the MatricesHarvester interface. When calling the STM propagator through this interface, I get the following error:

"<super: <class 'JavaError'>, <JavaError object>>\n    Java stacktrace:\njava.lang.ClassCastException: org.hipparchus.analysis.differentiation.Gradient cannot be cast to org.hipparchus.analysis.differentiation.DerivativeStructure\r\n\tat org.hipparchus.analysis.differentiation.DerivativeStructure.subtract(DerivativeStructure.java:65)\r\n\tat org.hipparchus.geometry.euclidean.threed.FieldVector3D.subtract(FieldVector3D.java:563)\r\n\tat org.orekit.forces.drag.DragForce.acceleration(DragForce.java:110)\r\n\tat org.orekit.propagation.numerical.StateTransitionMatrixGenerator.computePartials(StateTransitionMatrixGenerator.java:255)\r\n\tat org.orekit.propagation.numerical.StateTransitionMatrixGenerator.combinedDerivatives(StateTransitionMatrixGenerator.java:184)\r\n\tat org.orekit.propagation.integration.AbstractIntegratedPropagator$ConvertedSecondaryStateEquations.computeDerivatives(AbstractIntegratedPropagator.java:881)\r\n\tat org.hipparchus.ode.ExpandableODE.computeDerivatives(ExpandableODE.java:139)\r\n\tat org.hipparchus.ode.AbstractIntegrator.computeDerivatives(AbstractIntegrator.java:259)\r\n\tat org.hipparchus.ode.AbstractIntegrator.initIntegration(AbstractIntegrator.java:204)\r\n\tat org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:195)\r\n\tat org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:509)\r\n\tat org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:440)\r\n"

If I am reading this correctly, the method I am using to embed zeroed derivatives must return a gradient to mesh with the interface the MatricesHarvester requests. I’ve done some digging and experimenting, and I can’t get a working conversion to a Gradient-type derivative. Am I going about this the correct way, or is there an easier solution for this issue?

2

Hi,

Yes, you need to use Gradient, not DerivativeStructure. However you should be able to retrieve the correct Field from the input. So to make it run (without properly computing the derivatives), you can do something along these lines to create the vector:
FieldVector3D(date.getField(), velocity)

Cheers,
Romain

3

Thank you for the quick reply! I believe I understand what the correction you recommended does, though I am not clear how a velocity Field is being pulled from an absolutedate. Is this just to give it a functional field, or is this actually analytically correct assuming the analytic answer is zero?