Different result for the numerical propagation when using different orbital type

Orekit usage
1

Hello orekit community,

When I tried to test the Numerical Propagator from Matlab, I found an odd phenomenon: when the propagation type was changed from KEPLERIAN to CARTESIAN, the propagation results are different. Why?

Here is the code with Keplerian propagation type (propagationType = OrbitType.KEPLERIAN), which is transfered from the orekit tutorial.

clear
clc
import java.io.File;
import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
import org.hipparchus.util.FastMath;
import org.orekit.data.DataContext;
import org.orekit.data.DataProvidersManager;
import org.orekit.data.ZipJarCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.forces.ForceModel;
import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
import org.orekit.forces.gravity.potential.GravityFieldFactory;
import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
import org.orekit.forces.gravity.ThirdBodyAttraction
import org.orekit.forces.gravity.Relativity
import org.orekit.forces.radiation.IsotropicRadiationSingleCoefficient
import org.orekit.forces.radiation.SolarRadiationPressure
import org.orekit.forces.drag.IsotropicDrag
import org.orekit.forces.drag.DragForce
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.orbits.KeplerianOrbit;
import org.orekit.orbits.Orbit;
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.sampling.OrekitFixedStepHandler;
import org.orekit.propagation.sampling.OrekitStepNormalizer
import org.orekit.propagation.EphemerisGenerator;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.IERSConventions;
import org.orekit.bodies.CelestialBodyFactory
import org.orekit.models.earth.ReferenceEllipsoid
import org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation
import org.orekit.models.earth.atmosphere.NRLMSISE00
import org.orekit.models.earth.atmosphere.DTM2000
import org.orekit.models.earth.atmosphere.JB2008
import org.orekit.models.earth.atmosphere.SimpleExponentialAtmosphere
import org.orekit.models.earth.atmosphere.HarrisPriester
import org.orekit.propagation.numerical.PartialDerivativesEquations

%% Configure Orekit. The file orekit-data-master.zip
DM = org.orekit.data.DataContext.getDefault().getDataProvidersManager();
crawler = org.orekit.data.DirectoryCrawler(File('/home/tao/文档/Codes_Tao/orekit/orekit_11.3.2/orekit-data'));
DM.clearProviders();
DM.addProvider(crawler);

%% Initial orbit parameters
a = 24396159;    % semi major axis in meters
e = 0.72831215;  % eccentricity
i = (7.0)/180*pi; % inclination
omega = (180)/180*pi; % perigee argument
raan = (261)/180*pi; %right ascension of ascending node
lM = 0.0; % mean anomaly

%% Defines frames
eci = FramesFactory.getEME2000();
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
wgs84ellipsoid = ReferenceEllipsoid.getWgs84(itrf);

%% Initial date in UTC time scale
utc = TimeScalesFactory.getUTC();
initialDate = AbsoluteDate(2004, 01, 01, 23, 30, 00.000, utc);

%% Setup orbit propagator
mu = 3.986004415e+14;
% Orbit construction as Keplerian
initialOrbit = KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.MEAN, eci, initialDate, mu);
initialState = SpacecraftState(initialOrbit);

% Adaptive step integrator with a minimum step of 0.001 and a maximum step of 1000
minStep = 0.001;
maxstep = 1000.0;
positionTolerance = 1.0;
propagationType = OrbitType.KEPLERIAN;
tolerances = NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType);
integrator = DormandPrince853Integrator(minStep, maxstep, tolerances(1), tolerances(2));
integrator.setInitialStepSize(1);
% Propagator
propagator = NumericalPropagator(integrator);
propagator.setOrbitType(propagationType);

%% Force Model
% Earth gravity field with degree 64 and order 64
provider = GravityFieldFactory.getNormalizedProvider(64, 64);
holmesFeatherstone = HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010,true),provider);
propagator.addForceModel(holmesFeatherstone);

% Moon and Sun perturbations
moon = CelestialBodyFactory.getMoon();
sun = CelestialBodyFactory.getSun();
moon_3dbodyattraction = ThirdBodyAttraction(moon);
propagator.addForceModel(moon_3dbodyattraction);
sun_3dbodyattraction = ThirdBodyAttraction(sun);
propagator.addForceModel(sun_3dbodyattraction);

% Solar radiation pressure
isotropicRadiationSingleCoeff = IsotropicRadiationSingleCoefficient(2, 1.8);
solarRadiationPressure = SolarRadiationPressure(sun, wgs84ellipsoid.getEquatorialRadius(), isotropicRadiationSingleCoeff);
propagator.addForceModel(solarRadiationPressure);

% Atmospheric drag
atmosphere = HarrisPriester(sun, wgs84ellipsoid, 4);
isotropicDrag = IsotropicDrag(2, 2.2);
dragForce = DragForce(atmosphere, isotropicDrag);
propagator.addForceModel(dragForce);

% Relativity
relativityForce = Relativity(mu);
propagator.addForceModel(relativityForce);

% Set up initial state in the propagator
propagator.setInitialState(initialState);

% Set the propagator to ephemeris mode
generator = propagator.getEphemerisGenerator();

% Propagation with storage of the results in an integrated ephemeris
finalState = propagator.propagate(initialDate.shiftedBy(86400));
fprintf('Numerical propagation:\n  Final date: %s\n  %s\n',...
    finalState.getDate(), finalState.getOrbit());

The result is as follows:

Final date: 2004-01-02T23:30:00.000Z
  Keplerian parameters: {a: 2.4306954943147082E7; e: 0.7271675845162197; i: 6.992085368448246; pa: 180.89242772298732; raan: 260.5528338210272; v: 881.0653368532536;}
Final date: 2004-01-02T23:30:00.000Z
  Cartesian parameters: {P(-1.6852489192813065E7, -3.2562045930076305E7, -1383336.7413289594), V(981.4835023261265, -2083.233491686334, 160.67757106067373)}

However, when I replace the propagation type with CARTEISAN (propagationType = OrbitType.KARTESIAN), the results are different. With one-day propagation, the results are:

Final date: 2004-01-02T23:30:00.000Z
  Keplerian parameters: {a: 2.4306757298081696E7; e: 0.7271648865581377; i: 6.9921400647043965; pa: -179.10744495020242; raan: -99.44720347751696; v: 161.07023772189126;}
Final date: 2004-01-02T23:30:00.000Z
  Cartesian parameters: {P(-1.6850629930141047E7, -3.256536287103041E7, -1383053.1755013545), V(981.7338678315112, -2082.786221158252, 160.70027751420278)}

I want to know what leads to the different, when the force model is identical?

Thank you in advance

Xuefeng Tao.

2

Hi @taoxuefeng11,

Welcome to the community !

I think that either your integrator maxStep is too big or your positionTolerance is too small.

What must happen is that, since your integrator is too tolerant, it accepts steps that are too big.
Then, since Keplerian formalism is “more linear” than Cartesian formalism when propagating, the steps in Cartesian have a greater error that builds up during propagation, leading to the huge discrepancy you get after one day.

Try reducing the max step size to 100s for example, or lower the tolerance (but I had to get down to 1.e-9 to get similar results).

Best,
Maxime

2 Likes
3

Thanks very much! I try to reduce the max step size to 100s for both orbital type, and the results are as follows:

for propagationType = OrbitType.KEPLERIAN,

Final date: 2004-01-02T23:30:00.000Z
  Keplerian parameters: {a: 2.4306406960611366E7; e: 0.7271614884135825; i: 6.992123520864885; pa: 180.89260423807363; raan: 260.5527467554603; v: 881.071073244691;}
Final date: 2004-01-02T23:30:00.000Z
  Cartesian parameters: {P(-1.68501059535136E7, -3.2565504670960248E7, -1382979.9075461684), V(981.7949904140435, -2082.7169254899954, 160.7060949118338)}

for propagationType = OrbitType.KARTESIAN,

Final date: 2004-01-02T23:30:00.000Z
  Keplerian parameters: {a: 2.4306406962926593E7; e: 0.7271614884392142; i: 6.992123520621731; pa: -179.10739575947082; raan: -99.4472532447076; v: 161.07107325046488;}
Final date: 2004-01-02T23:30:00.000Z
  Cartesian parameters: {P(-1.6850105951840095E7, -3.256550467887109E7, -1382979.9071233398), V(981.7949906140076, -2082.7169248836517, 160.7060949188555)}

Now the results are nearly the same. Maybe in this instance, the max step size 1000 is too large for both KEPLERIAN and KARTESIAN type. Thanks again for your valuable suggestion!

Xuefeng Tao.

1 Like