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.