Hi, all,
I am trying to use NumericalPropagator to propagate Envisat’s trajectory for 24 hours, but compared with its CPF ephemeris, the errors of the result are really large. There must be something wrong in my code, but I don’t know where.
The initial Date is 2018-10-25 00:00:00.000. I added OceanTides, solid tides, third ocean attraction, solar radiation pressure, and Attitude modes.
Since Envisat is about 780km high, 8211kg, with dimensions 26m10m5m, I figured several parameters used in the code: cd = 2.2, area = 32.0 ( in fact the area = 2610, but as the default mass is 1000kg, I used area = 2610/8.211), cr =1.88
Here is my code:
public class MasterMode {
/** Program entry point.
* @param args program arguments (unused here)
*/
public static void main(String[] args) {
try {
// configure Orekit
File home = new File(System.getProperty("user.home"));
File orekitData = new File(home, "orekit-data");
if (!orekitData.exists()) {
System.err.format(Locale.US, "Failed to find %s folder%n",
orekitData.getAbsolutePath());
System.err.format(Locale.US, "You need to download %s from %s, unzip it in %s and rename it 'orekit-data' for this tutorial to work%n",
"orekit-data-master.zip", "https://gitlab.orekit.org/orekit/orekit-data/-/archive/master/orekit-data-master.zip",
home.getAbsolutePath());
System.exit(1);
}
DataProvidersManager manager = DataProvidersManager.getInstance();
manager.addProvider(new DirectoryCrawler(orekitData));
// gravitation coefficient
double mu = 3.986004415e+14;
// inertial frame
Frame inertialFrame = FramesFactory.getEME2000();
// Initial date
AbsoluteDate initialDate = new AbsoluteDate(2016, 10, 25, 0, 0, 0.000,
TimeScalesFactory.getUTC());
Vector3D p= new Vector3D(-3963665.96200000, 3057637.69600000, 5089564.99900000);
Vector3D v= new Vector3D(-2974.79114922258, 4691.50122498370, -5122.50267471032);
final double r2 = p.getNormSq();
final double r3 = r2 * FastMath.sqrt(r2);
PVCoordinates pv = new PVCoordinates(p,v);
System.out.println("current dir:" + System.getProperty("user.dir"));
Frame ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
Frame j2000 = FramesFactory.getEME2000();
PVCoordinates pv1 = ITRF.getTransformTo(j2000, initialDate).transformPVCoordinates(pv);
TimeStampedPVCoordinates pv2 = new TimeStampedPVCoordinates(initialDate, pv1.getPosition(), pv1.getVelocity(), new Vector3D(-mu / r3, pv1.getPosition()));
System.out.println(" orbit0 j2000: " + pv1);
Orbit initialOrbit = new CartesianOrbit(pv1, j2000, initialDate, mu);
// // Initial state definition
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 = 100.0;
final double positionTolerance = 5.0;
final OrbitType propagationType = OrbitType.CARTESIAN;
final double[][] tolerances =
NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType);
AdaptiveStepsizeIntegrator integrator =
new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
// Propagator
NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(propagationType);
// Force Model (reduced to perturbing gravity field)
final NormalizedSphericalHarmonicsProvider provider =
GravityFieldFactory.getNormalizedProvider(200, 200);
ForceModel holmesFeatherstone =
new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010,
true),
provider);
// Add force model to the propagator
propagator.addForceModel(holmesFeatherstone);
// ocean tides force model
final Frame bodyFrame;
bodyFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
final double equatorialRadius;
equatorialRadius = Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
final double flattening;
flattening = Constants.WGS84_EARTH_FLATTENING;
propagator.addForceModel(new OceanTides((new OneAxisEllipsoid(equatorialRadius, flattening, bodyFrame)).getBodyFrame(),
provider.getAe(), provider.getMu(),
1, 1, IERSConventions.IERS_2010,
TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true)));
// solid tides force model
List<CelestialBody> solidTidesBodies = new ArrayList<CelestialBody>();
solidTidesBodies.add(CelestialBodyFactory.getSun());
solidTidesBodies.add(CelestialBodyFactory.getMoon());
solidTidesBodies.add( CelestialBodyFactory.getJupiter());
//
// solidTidesBodies.add( CelestialBodyFactory.getMars());
//
//
// solidTidesBodies.add( CelestialBodyFactory.getMercury());
//
//
// solidTidesBodies.add( CelestialBodyFactory.getNeptune());
//
//
// solidTidesBodies.add( CelestialBodyFactory.getPluto());
//
//
// solidTidesBodies.add( CelestialBodyFactory.getSaturn());
solidTidesBodies.add( CelestialBodyFactory.getVenus());
propagator.addForceModel(new SolidTides((new OneAxisEllipsoid(equatorialRadius, flattening, bodyFrame)).getBodyFrame(),
provider.getAe(), provider.getMu(),
provider.getTideSystem(), IERSConventions.IERS_2010,
TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true),
solidTidesBodies.toArray(new CelestialBody[solidTidesBodies.size()])));
// third body attraction
propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
// drag
final double cd = 2.2;
final double area = 32.0;
final boolean cdEstimated = true;
MarshallSolarActivityFutureEstimation msafe =
new MarshallSolarActivityFutureEstimation("(?:Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec)\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit}F10\\.(?:txt|TXT)",
MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
// DataProvidersManager manager = DataProvidersManager.getInstance();
manager.feed(msafe.getSupportedNames(), msafe);
Atmosphere atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), new OneAxisEllipsoid(equatorialRadius, flattening, bodyFrame));
propagator.addForceModel(new DragForce(atmosphere, new IsotropicDrag(area, cd)));
// solar radiation pressure
final double cr = 1.88;
propagator.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(),
(new OneAxisEllipsoid(equatorialRadius, flattening, bodyFrame)).getEquatorialRadius(),
new IsotropicRadiationSingleCoefficient(area, cr)));
// post-Newtonian correction force due to general relativity
propagator.addForceModel(new Relativity(provider.getMu()));
final AttitudeMode mode;
mode = AttitudeMode.NADIR_POINTING_WITH_YAW_COMPENSATION;
// Set up initial state in the propagator
propagator.setInitialState(initialState);
// Set up operating mode for the propagator as master mode
// with fixed step and specialized step handler
propagator.setMasterMode(60., new TutorialStepHandler());
// Extrapolate from the initial to the final date
SpacecraftState finalState = propagator.propagate(initialDate.shiftedBy(86400.));
KeplerianOrbit o = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(finalState.getOrbit());
System.out.format(Locale.US, "Final state:%n%s %12.3f %10.8f %10.6f %10.6f %10.6f %10.6f%n",
finalState.getDate(),
o.getA(), o.getE(),
FastMath.toDegrees(o.getI()),
FastMath.toDegrees(o.getPerigeeArgument()),
FastMath.toDegrees(o.getRightAscensionOfAscendingNode()),
FastMath.toDegrees(o.getTrueAnomaly()));
} catch (OrekitException oe) {
System.err.println(oe.getMessage());
}
}
/** Specialized step handler.
* <p>This class extends the step handler in order to print on the output stream at the given step.<p>
* @author Pascal Parraud
*/
private static class TutorialStepHandler implements OrekitFixedStepHandler {
private PrintStream outputStream;
private final File outputFile = new File(new File(System.getProperty("user.home")), "result047.txt");
private TutorialStepHandler() {
//private constructor
}
public void init(final SpacecraftState s0, final AbsoluteDate t, final double step) {
try{
outputStream = new PrintStream(outputFile, "UTF-8");
}catch (IOException ioe) {
throw new OrekitException(ioe, LocalizedCoreFormats.SIMPLE_MESSAGE, ioe.getLocalizedMessage());
}
System.out.println(" date a e" +
" i \u03c9 \u03a9" +
" \u03bd");
}
public void handleStep(SpacecraftState currentState, boolean isLast) {
Frame ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
AbsoluteDate initialDate = currentState.getDate();
outputStream.format(Locale.US, "%s %12.3f %10.8f %10.6f %10.6f %10.6f %10.6f%n",
currentState.getDate(),
o.getPVCoordinates(initialDate, ITRF).getPosition().getX(),
o.getPVCoordinates(initialDate, ITRF).getPosition().getY(),
o.getPVCoordinates(initialDate, ITRF).getPosition().getZ(),
o.getPVCoordinates(initialDate, ITRF).getVelocity().getX(),
o.getPVCoordinates(initialDate, ITRF).getVelocity().getY(),
o.getPVCoordinates(initialDate, ITRF).getVelocity().getZ());
if (isLast) {
System.out.println("this was the last step ");
System.out.println();
}
}
}
/** Attitude modes. */
private static enum AttitudeMode {
NADIR_POINTING_WITH_YAW_COMPENSATION() {
public AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
{
return new YawCompensation(inertialFrame, new NadirPointing(inertialFrame, body));
}
},
CENTER_POINTING_WITH_YAW_STEERING {
public AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
{
return new YawSteering(inertialFrame,
new BodyCenterPointing(inertialFrame, body),
CelestialBodyFactory.getSun(),
Vector3D.PLUS_I);
}
},
LOF_ALIGNED_LVLH {
public AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
{
return new LofOffset(inertialFrame, LOFType.LVLH);
}
},
LOF_ALIGNED_QSW {
public AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
{
return new LofOffset(inertialFrame, LOFType.QSW);
}
},
LOF_ALIGNED_TNW {
public AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
{
return new LofOffset(inertialFrame, LOFType.TNW);
}
},
LOF_ALIGNED_VNC {
public AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
{
return new LofOffset(inertialFrame, LOFType.VNC);
}
},
LOF_ALIGNED_VVLH {
public AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
{
return new LofOffset(inertialFrame, LOFType.VVLH);
}
};
public abstract AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
;
}
}
The results written in the “result_047.txt” are in the ITRF frame. As Envisat is a SLR satellite, I compared the results with its CPF Ephemeris.The errors of x,y,z(ITRF) can reach up to almost 160m after 24 hrs. I am now trying to decrease the errors, but it seems there is nothing more I can do. Is there any perturbation I missed? Shall I change my propagator? Can anyone give me some advices?
By the way, I tried another satellite named Ajisai(altitude 1488km, weight 685 kg ), the errors of x,y,z(ITRF) can reach up to almost 20km after 24 hrs. I used parameters cd = 2.2, area = 3.2*0.685 ( it is a sphere with a diameter of 2.15m ), cr =1.88. The errors are really large……… I am clearly getting something wrong. I would be very grateful for any hints provided!