[Orekit Developers] STM/Jacobians w.r.t. Keplerian Parameters

["Ward, Evan" <Evan.Ward@nrl.navy.mil>]

Hi,

I'm trying to compute the State Transition Matrix (STM) for Keplerian orbit w.r.t. Keplerian parameters. I'm not sure if I'm doing it right so I thought I'd post my question here. So I'm trying to compute dx(t)/dx(0) where x = {a, e, i,  arg. of perigee, RAAN, M (mean anomaly)}. According to Montenbruck & Gill (section 7.1.1 of Satellite Orbits) the result should be the identity matrix except for the dM(t)/da(0).  term. In Orekit I've set the orbitType to Keplerian and the positionAngleType to Mean, but the result looks very different from what I would expect. I've attached my code and included the output below. The code is in groovy and it should be able to run if you want to try it. Am I doing something wrong? Is this a bug?

Initial State

Keplerian parameters: {a: 6678137.0; e: 0.1; i: 50.0; pa: 0.0; raan: 0.0; v: 0.0;}

Final state

Keplerian parameters: {a: 6678137.0; e: 0.1; i: 50.0; pa: 0.0; raan: 0.0; v: 359.9999999999999;}

STM

         1.084033316          1048746.583     -1.881271601e-07          338443.1765          217547.0804          689070.0095 

     8.986475977e-08        -0.4460115259      5.972999872e-14        -0.2344750310        -0.1507176447        -0.5273227631 

    -5.987960454e-21      1.482147738e-14          1.008352980     -2.359223927e-15       0.003748479827      5.162537065e-15 

    -1.049794901e-06          19.17181914      -0.009257359011          1.103560184        0.05508206389         -1.305994282 

     3.388131789e-21      9.503509091e-14        0.01440189399      9.742207041e-15          1.017867705      2.688127498e-14 

     7.973387020e-09         -23.67118846      1.041684976e-12        -0.8494761263        -0.5460327287         0.9198999327 

Expected dM(t)/da(t0)

-1.411288501683835E-6

Best Regards,

Evan

 

/**
 * @author Evan Ward
 */
@Grab(group = 'org.orekit', module = 'orekit', version = '9.0')
// for easy development only
@Grab(group = 'org.orekit', module = 'orekit', version = '9.0', classifier = 'sources')

import org.orekit.data.DataProvidersManager
import org.orekit.forces.gravity.NewtonianAttraction
import org.orekit.frames.FramesFactory
import org.orekit.orbits.KeplerianOrbit
import org.orekit.orbits.OrbitType
import org.orekit.orbits.PositionAngle
import org.orekit.propagation.SpacecraftState
import org.orekit.propagation.numerical.NumericalPropagator
import org.orekit.propagation.numerical.PartialDerivativesEquations
import org.orekit.time.AbsoluteDate
import org.orekit.utils.Constants;
import org.hipparchus.linear.MatrixUtils
import org.hipparchus.linear.RealMatrix
import org.hipparchus.ode.nonstiff.DormandPrince853Integrator


// set supplemental data
System.setProperty(DataProvidersManager.OREKIT_DATA_PATH, "/home/ward/eop/");
def eci = FramesFactory.getGCRF()
def epoch = new AbsoluteDate()

def positionAngle = PositionAngle.MEAN
def gm = Constants.EIGEN5C_EARTH_MU
def orbit = new KeplerianOrbit(6378137 + 300e3, 0.1, Math.toRadians(50), 0, 0, 0,
        positionAngle, eci, epoch, gm)
def orbitType = OrbitType.KEPLERIAN
def tol = NumericalPropagator.tolerances(1, orbit, orbitType)
def integrator = new DormandPrince853Integrator(1e-1, 3600, tol[0], tol[1])
def np = new NumericalPropagator(integrator)
np.setOrbitType(orbitType)
np.setPositionAngleType(positionAngle)
np.addForceModel(new NewtonianAttraction(gm))
def ic = new SpacecraftState(orbit)
def pde = new PartialDerivativesEquations("stm", np)
ic = pde.setInitialJacobians(ic)
np.setInitialState(ic)
def mapper = pde.getMapper()

def targetDate = epoch.shiftedBy(orbit.getKeplerianPeriod())
def state = np.propagate(targetDate)
def stm = new double[6][6]
mapper.getStateJacobian(state, stm)

println "Initial State"
println orbit
println "Final state"
println state.getOrbit()
println "STM"
for (def row  : stm) {
    for (def value  : row) {
        printf "%20.10g ", value
    }
    println ''
}
println "Expected dM(t)/da(t0)"
println(-3.0 / 2.0 * orbit.getKeplerianMeanMotion() / orbit.getA() * (targetDate.durationFrom(epoch)))

Attachment: smime.p7s
Description: S/MIME cryptographic signature