I use Orekit to write a piece of code to output the semi-major axis of every second in a customized period of time, and I want to verify whether it is consistent with the output data of STK. I find that the semi-major axis data at the same time point has a gap of 6-7km. I don’t know where I went wrong and how to modify it to make the output semi-major axis data of Orekit consistent with the data of STK.
public void test1() throws OrekitException, ParseException {
//加载固定数据
Dataloading dataloading = new Dataloading();
dataloading.dataLoade();
//时间设置
String startTime = "2025-03-03 08:00:00";
String endTime = "2025-03-04 15:00:00";
SimpleDateFormat format = new SimpleDateFormat("yyyy-MM-dd HH:mm:ss");
Date StartTime = format.parse(startTime);
Date EndTime = format.parse(endTime);
// Initial date in UTC time scale
final TimeScale utc = TimeScalesFactory.getUTC();
AbsoluteDate start = new AbsoluteDate(StartTime, utc);
AbsoluteDate end = new AbsoluteDate(EndTime, utc);
final AbsoluteDate initialDate = new AbsoluteDate(2025, 3, 3, 7, 23, 56.000, utc);
final double a = 6919145.67847071; // semi-major axis (m)
final double e = 0.00190973; // eccentricity
final double i = FastMath.toRadians(97.50841189); // inclination (°)
final double omega = FastMath.toRadians(169.64388288); // perigee argument (°)
final double raan = FastMath.toRadians(228.43982325); // right ascension of ascending node (°)
final double meanAnomaly = FastMath.toRadians(243.11338095); // mean anomaly (°)
// Inertial frame
final Frame inertialFrame = FramesFactory.getEME2000();
// gravitation coefficient
final double mu = 3.986004415e+14;
// Orbit construction as Keplerian
final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, meanAnomaly, PositionAngleType.MEAN, inertialFrame, initialDate, mu);
// Initialize state
double mass = 31; // kg
final SpacecraftState initialState = new SpacecraftState(initialOrbit, mass);
// steps limits
final double minStep = 0.001;
final double maxStep = 1000;
final double initStep = 60;
// error control parameters (absolute and relative)
final double positionError = 10.0;
final double[][] tolerances = NumericalPropagator.tolerances(positionError, initialOrbit, initialOrbit.getType());
// set up mathematical integrator
AdaptiveStepsizeIntegrator integrator =
new DormandPrince853Integrator(minStep, maxStep, tolerances[0], tolerances[1]);
integrator.setInitialStepSize(initStep);
//final double stepSize = 10;
//final AbstractIntegrator integrator = new ClassicalRungeKuttaIntegrator(stepSize);
final NumericalPropagator numericalPropagator = new NumericalPropagator(integrator);
// Earth and frame：定义地球框架
final Frame earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
final BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING, earthFrame);
int degree = 21;
int order = 21;
NormalizedSphericalHarmonicsProvider gravProvider = GravityFieldFactory.getNormalizedProvider(degree, order); // WGS84 potential coefficients
ForceModel gravModel = new HolmesFeatherstoneAttractionModel(earthFrame, gravProvider);
numericalPropagator.addForceModel(gravModel);
// Solid Earth Tides
SolidTides tides = new SolidTides(earthFrame,
Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
Constants.EIGEN5C_EARTH_MU,
TideSystem.ZERO_TIDE,
IERSConventions.IERS_2010,
TimeScalesFactory.getUT1(IERSConventions.IERS_2010, false),
CelestialBodyFactory.getSun(),
CelestialBodyFactory.getMoon());
numericalPropagator.addForceModel(tides);
// Relativity
Relativity rel = new Relativity(Constants.EIGEN5C_EARTH_MU);
numericalPropagator.addForceModel(rel);
// Third Body Attraction
numericalPropagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
numericalPropagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
// Solar Radiation Pressure Model
double area = 0.62396452; // m2
double cr = 1.0; // reflectance coefficient
SolarRadiationPressure srp =
new SolarRadiationPressure(CelestialBodyFactory.getSun(),
new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING, earthFrame),
new IsotropicRadiationSingleCoefficient(area, cr));
numericalPropagator.addForceModel(srp);
//Drag Model
double cd = 2.2; // drag coefficient
CssiSpaceWeatherData inputParameters = new CssiSpaceWeatherData("SpaceWeather-All-v1.2.txt");
Atmosphere atmosphere = new NRLMSISE00(inputParameters, CelestialBodyFactory.getSun(), earth);
ForceModel dragModel = new DragForce(atmosphere, new IsotropicDrag(area, cd));
numericalPropagator.addForceModel(dragModel);
// Initialize propagation
numericalPropagator.setInitialState(initialState);
numericalPropagator.setPositionAngleType(PositionAngleType.MEAN);
// Print the result
double step =1.0; //时间步长
for (double dt = 0; dt <= end.durationFrom(start); dt += step) {
AbsoluteDate date = start.shiftedBy(dt);
// Propagate from the initial date to the first raising or for the fixed duration
final SpacecraftState state = numericalPropagator.propagate(start, end);
System.out.println("时间："+date+"轨道高度："+state.getOrbit().getA());
}
}
