# Wrong MEAN orbital parameters computation

Hello!

While computing the mean orbital parameters of a LEO satellite, I noted that the results are not the ones I expected.
In order to compute them, I used the code on source code repository in src->turorials->java/fr/cs/examples->conversion-> PropagatorConversion, modified with LEO orbital parameters.
I was expecting a mean SMA of about 7179 km, while the result was of 7185 km.
I plotted the results and this is what I found out:

The mean of the elements is not correct. My guess on this is that, since the numerical propagator takes cartesian orbit type, the orbital parameters are averaged in cartesian and then re-converted in keplerian and this could create the mistake.

Let me know

Giacomo

Hi Giacomo,

Welcome to the Orekit community!

I’m not sure I understand what you’re trying to do.
How did you compute the mean orbital elements ?
The tutorial you’re pointing to is a conversion from numerical to Keplerian propagator. So there is no reference to mean orbital elements. Unless you consider the Keplerian elements to be your mean elements.

If you want to compute mean elements I suggest you try the DSST propagator.

Do you have a sample code of what you’re aiming to do ?

Maxime

Yes, I consider keplerian elements to be my mean elements (they should be the mean elements, aren’t they?).
The idea is to use a numerical propagator with a full force model, and than to convert the propagator into a keplerian so that I can obtain the keplerian (constant) orbital parameters.
Hereafter I attached the code that produces the results of the previous figure:

public class PropagatorConversion {

public static void main(String[] args) {
//Configuration data//
File orekitData = new File("orekit-data");
DataProvidersManager manager = DataProvidersManager.getInstance();
//////////////////////

// gravity field
NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(2, 0);
double mu =  provider.getMu();

// inertial frame
Frame inertialFrame = FramesFactory.getEME2000();

// Initial date
AbsoluteDate initialDate = new AbsoluteDate(2022,11,01,1,48,9.303,
TimeScalesFactory.getUTC());
// Initial orbit
final double a     = 7182.25536152*1000;                // semi major axis in meters
final double e     = 0.00111776;                        // eccentricity
final double i     = FastMath.toRadians(98.73774250);   // inclination
final double omega = FastMath.toRadians(65.11464822);   // perigee argument
final double raan  = FastMath.toRadians(2.23553023);    // right ascention of ascending node
final double TA    = FastMath.toRadians(151.17294416);  // true anomaly
Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, TA, PositionAngle.TRUE,
inertialFrame, initialDate, mu);
final double period = initialOrbit.getKeplerianPeriod();

// Initial state definition
final SpacecraftState initialState = new SpacecraftState(initialOrbit);

// Adaptive step integrator with a minimum step of 0.001 and a maximum step of 1000
final double minStep    = 0.001;
final double maxStep    = 1000.;
final double dP         = 1.e-2;
final OrbitType orbType = OrbitType.CARTESIAN;
final double[][] tol = NumericalPropagator.tolerances(dP, initialOrbit, orbType);
final AbstractIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep,
tol[0], tol[1]);
// Propagator
NumericalPropagator numProp = new NumericalPropagator(integrator);
numProp.setInitialState(initialState);
numProp.setOrbitType(orbType);

//Propagator Force Model//////////////////////////////////////////////////////////////////////////////////////////////////                                                                            //
//1) Gravity Field & Third Body                                                                                         //
//
ForceModel holmesFeatherstone =                                                                                         //
new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010,true),provider);  //
//
//2) Drag                                                                                                               //
final DragForce drag =                                                                                                  //
new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),                                                     //
Constants.WGS84_EARTH_FLATTENING,                             //
FramesFactory.getITRF(IERSConventions.IERS_2010, true))),     //
new IsotropicDrag(2.5, 1.2));                                                                         //
//
ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();                                                      //
OneAxisEllipsoid earth =                                                                                                //
new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765,                                                                   //
FramesFactory.getITRF(IERSConventions.IERS_2010, true));                                       //                                                                                     //
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// Propagator factory
PropagatorBuilder builder = new KeplerianPropagatorBuilder(initialOrbit, PositionAngle.TRUE, dP);

// Propagator converter
PropagatorConverter fitter = new FiniteDifferencePropagatorConverter(builder, 1.e-6, 5000);

// Resulting propagator
KeplerianPropagator kepProp = (KeplerianPropagator)fitter.convert(numProp, 2*period, 251);

// Step handlers
StatesHandler numStepHandler = new StatesHandler();
StatesHandler kepStepHandler = new StatesHandler();

// Set up operating mode for the propagator as master mode
// with fixed step and specialized step handler
numProp.setMasterMode(60., numStepHandler);
kepProp.setMasterMode(60., kepStepHandler);

// Extrapolate from the initial to the final date
numProp.propagate(initialDate, initialDate.shiftedBy(10.*period));
kepProp.propagate(initialDate, initialDate.shiftedBy(10.*period));
// retrieve the states
List<SpacecraftState> numStates = numStepHandler.getStates();
List<SpacecraftState> kepStates = kepStepHandler.getStates();

// Print the results on the output file
System.out.println("Mean SMA (constant):"+ kepProp.getInitialState().getA());
System.out.println("");
System.out.println("Osculating SMA");
for (SpacecraftState numState : numStates) {
for (SpacecraftState kepState : kepStates) {
if (numState.getDate().compareTo(kepState.getDate()) == 0) {
System.out.println(numState.getA());
break;
}
}
}
}

/** Specialized step handler.
* <p>This class extends the step handler in order to handle states at each step.<p>
*/
private static class StatesHandler implements OrekitFixedStepHandler {

/** Points. */
private final List<SpacecraftState> states;

public StatesHandler() {
// prepare an empty list of SpacecraftState
states = new ArrayList<SpacecraftState>();
}

public void handleStep(SpacecraftState currentState, boolean isLast) {

}

/** Get the list of handled orbital elements.
* @return orbital elements list
*/
public List<SpacecraftState> getStates() {
return states;
}

}


}

Thank you!!

Giacomo

The mean parameters in the sense of some dynamical model is not the arithmetic mean.
This is more obvious if instead of using a low eccentricity orbit you use a very eccentric one (say a geostationary transfer orbit, with e ≈ 0.72. In thi case, instead of having an osculating semi major axis that look almost sinusoidal, you would end up with a curve with a very flat bottom and high peaks repeating at a regular rate. The flat bottom would correspond to apogee crossings and the high peaks to perigee crossings, as the J_2 zonal term effect is much larger when you are close to Earth, so it distorts the orbit a lot. The “mean” meaning here would rather be the a that you can compute using Keplerian equation T=2\pi\sqrt{a^3/\mu} and using the separation of the peaks as the T. You would than see that this mean semi major axis in much closer to the semi major axis near apogee than to the arithmetic mean.

Anyway, using the propagator convertor to get the mean elements as you did is a proper way, and if it gives back a semi major axis that is not the arithmetic mean, it is probably right. One thing you could check is to compare the Cartesian distance between the output of the three models: the numerical one, the converted Keplerian one and the Keplerian that uses the arithmetic mean. The comparison is simpler to understand if done in a local orbital frame (like LofType.VNC). I guess you will see an ever increasing difference along the tangential direction between the two Keplerian models.