I referred to the DSSTPropagation tutorial in 12.2 and wrote a piece of code, passing in the six elements of the satellite orbit and outputting the average elements of the semi-major axis. However, it didn’t work. My PropagationType outputType was indeed propagationtype.mean. The following is part of the code. Also, the program runs a bit slowly.
public static void inputOrbitHeight(List<KeplerOrbit> keplerlist) throws IOException {
for (KeplerOrbit keplerOrbit : keplerlist) {
// check mandatory input parameters
if (keplerOrbit.getOrbitDate() == null) {
throw new IOException("Orbit date is not defined.");
}
// All dates in UTC
final TimeScale utc = TimeScalesFactory.getUTC();
//地球自转速率
final double rotationRate = Constants.WGS84_EARTH_ANGULAR_VELOCITY;
// Propagator data
final int degree = 6;
final int order = 6;
// Potential coefficients providers
final AbsoluteDate orbitDate = new AbsoluteDate(keplerOrbit.getOrbitDate(), utc);
final UnnormalizedSphericalHarmonicsProvider unnormalized =
GravityFieldFactory.getConstantUnnormalizedProvider(degree, order, orbitDate);
final NormalizedSphericalHarmonicsProvider normalized =
GravityFieldFactory.getConstantNormalizedProvider(degree, order, orbitDate);
// Central body attraction coefficient (m³/s²)
final double mu = unnormalized.getMu();
// Earth frame definition
final Frame earthFrame= FramesFactory.getITRF(IERSConventions.IERS_2010, true);
// Orbit definition
final Frame frame = FramesFactory.getEME2000();
final Orbit orbit = new KeplerianOrbit(keplerOrbit.getSemiMajorAxis(),
keplerOrbit.getEccentricity(),
keplerOrbit.getInclination(),
keplerOrbit.getPerigeeArgument(),
keplerOrbit.getRaan(),
keplerOrbit.getMeanAnomaly(),
PositionAngleType.MEAN, frame, orbitDate, mu);
// DSST propagator definition
final PropagationType initialType =PropagationType.OSCULATING;
final PropagationType outputType = PropagationType.MEAN;
List<String> shortPeriodCoefficients = null;
double mass = 31;
double fixedStepSize = -1.;
double minStep = 6000.0;
double maxStep = 86400.0;
double dP = 1.0;
final DSSTPropagator dsstProp = createDSSTProp(orbit, mass,
initialType, outputType,
fixedStepSize, minStep, maxStep, dP,
shortPeriodCoefficients);
dsstProp.setInterpolationGridToFixedNumberOfPoints(3);
final double ae = unnormalized.getAe();
// Central Body Force Model with un-normalized coefficients
dsstProp.addForceModel(new DSSTZonal(unnormalized));
dsstProp.addForceModel(new DSSTTesseral(earthFrame, rotationRate, unnormalized));
// 3rd body
dsstProp.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getSun(), mu));
dsstProp.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getMoon(), mu));
// Drag
final double Cd =2.2;
final double dragArea = 0.13;
final OneAxisEllipsoid earth = new OneAxisEllipsoid(ae, Constants.WGS84_EARTH_FLATTENING, earthFrame);
final Atmosphere atm = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
dsstProp.addForceModel(new DSSTAtmosphericDrag(atm, Cd, dragArea, mu));
// Solar Radiation Pressure
final double Cr = 1.0;
final double srpArea = 0.13;
dsstProp.addForceModel(new DSSTSolarRadiationPressure(Cr, srpArea, CelestialBodyFactory.getSun(), earth, mu));
// Simulation properties
final AbsoluteDate start = orbit.getDate();
Date currentTime = new Date();
AbsoluteDate currentDate = new AbsoluteDate(currentTime, utc);
// DSST Propagation
final EphemerisGenerator dsstGenerator = dsstProp.getEphemerisGenerator();
final double dsstOn = System.currentTimeMillis();
dsstProp.propagate(start, currentDate);
final double dsstOff = System.currentTimeMillis();
System.out.println("DSST execution time (without large file write) : " + (dsstOff - dsstOn) / 1000.);
System.out.println("writing file...");
final BoundedPropagator dsstEphem = dsstGenerator.getGeneratedEphemeris();
// dsstEphem.getMultiplexer().add(outStep, new DSSTPropagation.OutputHandler(output,
// displayKeplerian, displayEquinoctial, displayCartesian,
// null));
SpacecraftState intermediateState = dsstEphem.propagate(start, currentDate);
double a =intermediateState.getOrbit().getA();
String satHieght = Double.toString((a-Constants.WGS84_EARTH_EQUATORIAL_RADIUS)/1000);
System.out.println(intermediateState.getOrbit().getDate().toString(utc)+" satheigh："+satHieght+"Km"+" a："+a+"m");
}
}