Java tutorial
/* Copyright 2002-2015 CS Systmes d'Information * Licensed to CS Systmes d'Information (CS) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. * CS licenses this file to You under the Apache License, Version 2.0 * (the "License"); you may not use this file except in compliance with * the License. You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ package fr.cs.examples.propagation; import java.io.BufferedWriter; import java.io.File; import java.io.FileInputStream; import java.io.FileWriter; import java.io.IOException; import java.text.ParseException; import java.util.ArrayList; import java.util.Formatter; import java.util.List; import java.util.Locale; import java.util.NoSuchElementException; import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; import org.apache.commons.math3.ode.AbstractIntegrator; import org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator; import org.apache.commons.math3.ode.nonstiff.ClassicalRungeKuttaIntegrator; import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator; import org.apache.commons.math3.util.FastMath; import org.apache.commons.math3.util.MathUtils; import org.orekit.Utils; import org.orekit.bodies.CelestialBodyFactory; import org.orekit.bodies.OneAxisEllipsoid; import org.orekit.errors.OrekitException; import org.orekit.forces.SphericalSpacecraft; import org.orekit.forces.drag.Atmosphere; import org.orekit.forces.drag.DragForce; import org.orekit.forces.drag.HarrisPriester; import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel; import org.orekit.forces.gravity.ThirdBodyAttraction; import org.orekit.forces.gravity.potential.GravityFieldFactory; import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider; import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider; import org.orekit.forces.radiation.SolarRadiationPressure; import org.orekit.frames.Frame; import org.orekit.frames.FramesFactory; import org.orekit.orbits.CartesianOrbit; import org.orekit.orbits.CircularOrbit; import org.orekit.orbits.EquinoctialOrbit; import org.orekit.orbits.KeplerianOrbit; import org.orekit.orbits.Orbit; 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.semianalytical.dsst.DSSTPropagator; import org.orekit.propagation.semianalytical.dsst.forces.DSSTAtmosphericDrag; import org.orekit.propagation.semianalytical.dsst.forces.DSSTCentralBody; import org.orekit.propagation.semianalytical.dsst.forces.DSSTSolarRadiationPressure; import org.orekit.propagation.semianalytical.dsst.forces.DSSTThirdBody; import org.orekit.time.AbsoluteDate; import org.orekit.time.TimeScale; import org.orekit.time.TimeScalesFactory; import org.orekit.utils.Constants; import org.orekit.utils.PVCoordinates; import fr.cs.examples.KeyValueFileParser; /** Orekit tutorial for semi-analytical extrapolation using the DSST. * <p> * The parameters are read from the input file dsst-propagation.in located in the user's * home directory (see commented example at src/tutorial/ressources/dsst-propagation.in). * The results are written to the ouput file dsst-propagation.out in the same directory. * </p> * <p> * Comparison between the DSST propagator and the numerical propagator can be optionally * performed. Numerical results are written to the ouput file numerical-propagation.out. * </p> * * @author Romain Di Costanzo * @author Pascal Parraud */ public class DSSTPropagation { /** Program entry point. * @param args program arguments */ public static void main(String[] args) { try { // configure Orekit data acces Utils.setDataRoot("tutorial-orekit-data"); // input/output (in user's home directory) File input = new File(new File(System.getProperty("user.home")), "dsst-propagation.in"); File output = new File(input.getParentFile(), "dsst-propagation.out"); new DSSTPropagation().run(input, output); } catch (IOException ioe) { System.err.println(ioe.getLocalizedMessage()); System.exit(1); } catch (IllegalArgumentException iae) { System.err.println(iae.getLocalizedMessage()); System.exit(1); } catch (OrekitException oe) { System.err.println(oe.getLocalizedMessage()); System.exit(1); } catch (ParseException pe) { System.err.println(pe.getLocalizedMessage()); System.exit(1); } } /** Input parameter keys. */ private static enum ParameterKey { ORBIT_DATE, ORBIT_CIRCULAR_A, ORBIT_CIRCULAR_EX, ORBIT_CIRCULAR_EY, ORBIT_CIRCULAR_I, ORBIT_CIRCULAR_RAAN, ORBIT_CIRCULAR_ALPHA, ORBIT_EQUINOCTIAL_A, ORBIT_EQUINOCTIAL_EX, ORBIT_EQUINOCTIAL_EY, ORBIT_EQUINOCTIAL_HX, ORBIT_EQUINOCTIAL_HY, ORBIT_EQUINOCTIAL_LAMBDA, ORBIT_KEPLERIAN_A, ORBIT_KEPLERIAN_E, ORBIT_KEPLERIAN_I, ORBIT_KEPLERIAN_PA, ORBIT_KEPLERIAN_RAAN, ORBIT_KEPLERIAN_ANOMALY, ORBIT_ANGLE_TYPE, ORBIT_CARTESIAN_PX, ORBIT_CARTESIAN_PY, ORBIT_CARTESIAN_PZ, ORBIT_CARTESIAN_VX, ORBIT_CARTESIAN_VY, ORBIT_CARTESIAN_VZ, ORBIT_IS_OSCULATING, START_DATE, DURATION, DURATION_IN_DAYS, OUTPUT_STEP, FIXED_INTEGRATION_STEP, NUMERICAL_COMPARISON, CENTRAL_BODY_ORDER, CENTRAL_BODY_DEGREE, THIRD_BODY_MOON, THIRD_BODY_SUN, MASS, DRAG, DRAG_CD, DRAG_SF, SOLAR_RADIATION_PRESSURE, SOLAR_RADIATION_PRESSURE_CR, SOLAR_RADIATION_PRESSURE_SF; } private void run(final File input, final File output) throws IOException, IllegalArgumentException, OrekitException, ParseException { // read input parameters KeyValueFileParser<ParameterKey> parser = new KeyValueFileParser<ParameterKey>(ParameterKey.class); parser.parseInput(new FileInputStream(input)); // check mandatory input parameters if (!parser.containsKey(ParameterKey.ORBIT_DATE)) { throw new IOException("Orbit date is not defined."); } if (!parser.containsKey(ParameterKey.DURATION) && !parser.containsKey(ParameterKey.DURATION_IN_DAYS)) { throw new IOException("Propagation duration is not defined."); } // All dates in UTC final TimeScale utc = TimeScalesFactory.getUTC(); final int degree = parser.getInt(ParameterKey.CENTRAL_BODY_DEGREE); final int order = FastMath.min(degree, parser.getInt(ParameterKey.CENTRAL_BODY_ORDER)); // Potential coefficients providers final UnnormalizedSphericalHarmonicsProvider unnormalized = GravityFieldFactory .getConstantUnnormalizedProvider(degree, order); final NormalizedSphericalHarmonicsProvider normalized = GravityFieldFactory .getConstantNormalizedProvider(degree, order); // Central body attraction coefficient (m/s) final double mu = unnormalized.getMu(); // Earth frame definition (for faster computation) final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame(); // Orbit definition (inertial frame is EME2000) final Orbit orbit = createOrbit(parser, FramesFactory.getEME2000(), utc, mu); // DSST propagator definition double mass = 1000.0; if (parser.containsKey(ParameterKey.MASS)) { mass = parser.getDouble(ParameterKey.MASS); } Boolean isOsculating = false; if (parser.containsKey(ParameterKey.ORBIT_IS_OSCULATING)) { isOsculating = parser.getBoolean(ParameterKey.ORBIT_IS_OSCULATING); } double fixedStepSize = -1.; if (parser.containsKey(ParameterKey.FIXED_INTEGRATION_STEP)) { fixedStepSize = parser.getDouble(ParameterKey.FIXED_INTEGRATION_STEP); } final DSSTPropagator dsstProp = createDSSTProp(orbit, mass, isOsculating, fixedStepSize); // Set Force models setForceModel(parser, unnormalized, earthFrame, dsstProp); // Simulation properties AbsoluteDate start; if (parser.containsKey(ParameterKey.START_DATE)) { start = parser.getDate(ParameterKey.START_DATE, utc); } else { start = parser.getDate(ParameterKey.ORBIT_DATE, utc); } double duration = 0.; if (parser.containsKey(ParameterKey.DURATION)) { duration = parser.getDouble(ParameterKey.DURATION); } if (parser.containsKey(ParameterKey.DURATION_IN_DAYS)) { duration = parser.getDouble(ParameterKey.DURATION_IN_DAYS) * Constants.JULIAN_DAY; } double outStep = parser.getDouble(ParameterKey.OUTPUT_STEP); // Add orbit handler OrbitHandler dsstHandler = new OrbitHandler(); dsstProp.setMasterMode(outStep, dsstHandler); // DSST Propagation final double dsstOn = System.currentTimeMillis(); dsstProp.propagate(start, start.shiftedBy(duration)); final double dsstOff = System.currentTimeMillis(); System.out.println("DSST execution time: " + (dsstOff - dsstOn) / 1000.); // Print results printOutput(output, dsstHandler, start); System.out.println("DSST results saved as file " + output); // Check if we want to compare numerical to DSST propagator (default is false) if (parser.containsKey(ParameterKey.NUMERICAL_COMPARISON) && parser.getBoolean(ParameterKey.NUMERICAL_COMPARISON)) { if (!isOsculating) { System.out.println("\nWARNING:"); System.out.println( "The DSST propagator considers a mean orbit while the numerical will consider an osculating one."); System.out.println("The comparison will be meaningless.\n"); } // output (in user's home directory) File output_num = new File(input.getParentFile(), "numerical-propagation.out"); // Numerical propagator definition final NumericalPropagator numProp = createNumProp(orbit, mass); // Set Force models setForceModel(parser, normalized, earthFrame, numProp); // Add orbit handler OrbitHandler numHandler = new OrbitHandler(); numProp.setMasterMode(outStep, numHandler); // Numerical Propagation final double numOn = System.currentTimeMillis(); numProp.propagate(start, start.shiftedBy(duration)); final double numOff = System.currentTimeMillis(); System.out.println("Numerical execution time: " + (numOff - numOn) / 1000.); // Print results printOutput(output_num, numHandler, start); System.out.println("Numerical results saved as file " + output_num); } } /** Create an orbit from input parameters * @param parser input file parser * @param frame inertial frame * @param scale time scale * @param mu central attraction coefficient * @throws NoSuchElementException if input parameters are mising * @throws IOException if input parameters are invalid */ private Orbit createOrbit(final KeyValueFileParser<ParameterKey> parser, final Frame frame, final TimeScale scale, final double mu) throws NoSuchElementException, IOException { // Orbit definition Orbit orbit; PositionAngle angleType = PositionAngle.MEAN; if (parser.containsKey(ParameterKey.ORBIT_ANGLE_TYPE)) { angleType = PositionAngle.valueOf(parser.getString(ParameterKey.ORBIT_ANGLE_TYPE).toUpperCase()); } if (parser.containsKey(ParameterKey.ORBIT_KEPLERIAN_A)) { orbit = new KeplerianOrbit(parser.getDouble(ParameterKey.ORBIT_KEPLERIAN_A) * 1000., parser.getDouble(ParameterKey.ORBIT_KEPLERIAN_E), parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_I), parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_PA), parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_RAAN), parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_ANOMALY), angleType, frame, parser.getDate(ParameterKey.ORBIT_DATE, scale), mu); } else if (parser.containsKey(ParameterKey.ORBIT_EQUINOCTIAL_A)) { orbit = new EquinoctialOrbit(parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_A) * 1000., parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_EX), parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_EY), parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_HX), parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_HY), parser.getAngle(ParameterKey.ORBIT_EQUINOCTIAL_LAMBDA), angleType, frame, parser.getDate(ParameterKey.ORBIT_DATE, scale), mu); } else if (parser.containsKey(ParameterKey.ORBIT_CIRCULAR_A)) { orbit = new CircularOrbit(parser.getDouble(ParameterKey.ORBIT_CIRCULAR_A) * 1000., parser.getDouble(ParameterKey.ORBIT_CIRCULAR_EX), parser.getDouble(ParameterKey.ORBIT_CIRCULAR_EY), parser.getAngle(ParameterKey.ORBIT_CIRCULAR_I), parser.getAngle(ParameterKey.ORBIT_CIRCULAR_RAAN), parser.getAngle(ParameterKey.ORBIT_CIRCULAR_ALPHA), angleType, frame, parser.getDate(ParameterKey.ORBIT_DATE, scale), mu); } else if (parser.containsKey(ParameterKey.ORBIT_CARTESIAN_PX)) { final double[] pos = { parser.getDouble(ParameterKey.ORBIT_CARTESIAN_PX) * 1000., parser.getDouble(ParameterKey.ORBIT_CARTESIAN_PY) * 1000., parser.getDouble(ParameterKey.ORBIT_CARTESIAN_PZ) * 1000. }; final double[] vel = { parser.getDouble(ParameterKey.ORBIT_CARTESIAN_VX) * 1000., parser.getDouble(ParameterKey.ORBIT_CARTESIAN_VY) * 1000., parser.getDouble(ParameterKey.ORBIT_CARTESIAN_VZ) * 1000. }; orbit = new CartesianOrbit(new PVCoordinates(new Vector3D(pos), new Vector3D(vel)), frame, parser.getDate(ParameterKey.ORBIT_DATE, scale), mu); } else { throw new IOException("Orbit definition is incomplete."); } return orbit; } /** Set up the DSST Propagator * * @param orbit initial orbit * @param mass S/C mass (kg) * @param isOsculating if orbital elements are osculating * @param fixedStepSize step size for fixed step integrator (s) * @throws OrekitException */ private DSSTPropagator createDSSTProp(final Orbit orbit, final double mass, final boolean isOsculating, final double fixedStepSize) throws OrekitException { AbstractIntegrator integrator; if (fixedStepSize > 0.) { integrator = new ClassicalRungeKuttaIntegrator(fixedStepSize); } else { final double minStep = orbit.getKeplerianPeriod(); final double maxStep = minStep * 100.; final double[][] tol = DSSTPropagator.tolerances(1.0, orbit); integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]); ((AdaptiveStepsizeIntegrator) integrator).setInitialStepSize(10. * minStep); } DSSTPropagator dsstProp = new DSSTPropagator(integrator); dsstProp.setInitialState(new SpacecraftState(orbit, mass), isOsculating); return dsstProp; } /** Create the numerical propagator * * @param orbit initial orbit * @param mass S/C mass (kg) * @throws OrekitException */ private NumericalPropagator createNumProp(final Orbit orbit, final double mass) throws OrekitException { final double[][] tol = NumericalPropagator.tolerances(1.0, orbit, orbit.getType()); final double minStep = 1.e-3; final double maxStep = 1.e+3; AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]); integrator.setInitialStepSize(100.); NumericalPropagator numProp = new NumericalPropagator(integrator); numProp.setInitialState(new SpacecraftState(orbit, mass)); return numProp; } /** Set DSST propagator force models * * @param parser input file parser * @param unnormalized spherical harmonics provider * @param earthFrame Earth rotating frame * @param dsstProp DSST propagator * @throws IOException * @throws OrekitException */ private void setForceModel(final KeyValueFileParser<ParameterKey> parser, final UnnormalizedSphericalHarmonicsProvider unnormalized, final Frame earthFrame, final DSSTPropagator dsstProp) throws IOException, OrekitException { final double ae = unnormalized.getAe(); final int degree = parser.getInt(ParameterKey.CENTRAL_BODY_DEGREE); final int order = parser.getInt(ParameterKey.CENTRAL_BODY_ORDER); if (order > degree) { throw new IOException("Potential order cannot be higher than potential degree"); } // Central Body Force Model with un-normalized coefficients dsstProp.addForceModel( new DSSTCentralBody(earthFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, unnormalized)); // 3rd body (SUN) if (parser.containsKey(ParameterKey.THIRD_BODY_SUN) && parser.getBoolean(ParameterKey.THIRD_BODY_SUN)) { dsstProp.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getSun())); } // 3rd body (MOON) if (parser.containsKey(ParameterKey.THIRD_BODY_MOON) && parser.getBoolean(ParameterKey.THIRD_BODY_MOON)) { dsstProp.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getMoon())); } // Drag if (parser.containsKey(ParameterKey.DRAG) && parser.getBoolean(ParameterKey.DRAG)) { 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, parser.getDouble(ParameterKey.DRAG_CD), parser.getDouble(ParameterKey.DRAG_SF))); } // Solar Radiation Pressure if (parser.containsKey(ParameterKey.SOLAR_RADIATION_PRESSURE) && parser.getBoolean(ParameterKey.SOLAR_RADIATION_PRESSURE)) { dsstProp.addForceModel(new DSSTSolarRadiationPressure( parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_CR), parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_SF), CelestialBodyFactory.getSun(), ae)); } } /** Set numerical propagator force models * * @param parser input file parser * @param normalized spherical harmonics provider * @param earthFrame Earth rotating frame * @param numProp numerical propagator * @throws IOException * @throws OrekitException */ private void setForceModel(final KeyValueFileParser<ParameterKey> parser, final NormalizedSphericalHarmonicsProvider normalized, final Frame earthFrame, final NumericalPropagator numProp) throws IOException, OrekitException { final double ae = normalized.getAe(); final int degree = parser.getInt(ParameterKey.CENTRAL_BODY_DEGREE); final int order = parser.getInt(ParameterKey.CENTRAL_BODY_ORDER); if (order > degree) { throw new IOException("Potential order cannot be higher than potential degree"); } // Central Body (normalized coefficients) numProp.addForceModel(new HolmesFeatherstoneAttractionModel(earthFrame, normalized)); // 3rd body (SUN) if (parser.containsKey(ParameterKey.THIRD_BODY_SUN) && parser.getBoolean(ParameterKey.THIRD_BODY_SUN)) { numProp.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun())); } // 3rd body (MOON) if (parser.containsKey(ParameterKey.THIRD_BODY_MOON) && parser.getBoolean(ParameterKey.THIRD_BODY_MOON)) { numProp.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon())); } // Drag if (parser.containsKey(ParameterKey.DRAG) && parser.getBoolean(ParameterKey.DRAG)) { final OneAxisEllipsoid earth = new OneAxisEllipsoid(ae, Constants.WGS84_EARTH_FLATTENING, earthFrame); final Atmosphere atm = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6); final SphericalSpacecraft ssc = new SphericalSpacecraft(parser.getDouble(ParameterKey.DRAG_SF), parser.getDouble(ParameterKey.DRAG_CD), 0., 0.); numProp.addForceModel(new DragForce(atm, ssc)); } // Solar Radiation Pressure if (parser.containsKey(ParameterKey.SOLAR_RADIATION_PRESSURE) && parser.getBoolean(ParameterKey.SOLAR_RADIATION_PRESSURE)) { final double cR = parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_CR); // cR being the DSST SRP coef and assuming a spherical spacecraft, the conversion is: // cR = 1 + (1 - kA) * (1 - kR) * 4 / 9 // with kA arbitrary sets to 0 final double kR = 3.25 - 2.25 * cR; final SphericalSpacecraft ssc = new SphericalSpacecraft( parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_SF), 0., 0., kR); numProp.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(), ae, ssc)); } } /** Print the results in the output file * * @param handler orbit handler * @param output output file * @param sart start date of propagation * @throws OrekitException * @throws IOException */ private void printOutput(final File output, final OrbitHandler handler, final AbsoluteDate start) throws OrekitException, IOException { // Output format: // time_from_start, a, e, i, raan, pa, aM, h, k, p, q, lM, px, py, pz, vx, vy, vz final String format = new String( " %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e"); final BufferedWriter buffer = new BufferedWriter(new FileWriter(output)); buffer.write( "## time_from_start(s) a(km) e i(deg) "); buffer.write( " raan(deg) pa(deg) mean_anomaly(deg) ey/h "); buffer.write( " ex/k hy/p hx/q mean_longitude_arg(deg)"); buffer.write( " Xposition(km) Yposition(km) Zposition(km) Xvelocity(km/s) "); buffer.write(" Yvelocity(km/s) Zvelocity(km/s)"); buffer.newLine(); for (Orbit o : handler.getOrbits()) { final Formatter f = new Formatter(new StringBuilder(), Locale.ENGLISH); // Time from start (s) final double time = o.getDate().durationFrom(start); // Semi-major axis (km) final double a = o.getA() / 1000.; // Keplerian elements // Eccentricity final double e = o.getE(); // Inclination (degrees) final double i = Math.toDegrees(MathUtils.normalizeAngle(o.getI(), FastMath.PI)); // Right Ascension of Ascending Node (degrees) KeplerianOrbit ko = new KeplerianOrbit(o); final double ra = Math .toDegrees(MathUtils.normalizeAngle(ko.getRightAscensionOfAscendingNode(), FastMath.PI)); // Perigee Argument (degrees) final double pa = Math.toDegrees(MathUtils.normalizeAngle(ko.getPerigeeArgument(), FastMath.PI)); // Mean Anomaly (degrees) final double am = Math .toDegrees(MathUtils.normalizeAngle(ko.getAnomaly(PositionAngle.MEAN), FastMath.PI)); // Equinoctial elements // ey/h component of eccentricity vector final double h = o.getEquinoctialEy(); // ex/k component of eccentricity vector final double k = o.getEquinoctialEx(); // hy/p component of inclination vector final double p = o.getHy(); // hx/q component of inclination vector final double q = o.getHx(); // Mean Longitude Argument (degrees) final double lm = Math.toDegrees(MathUtils.normalizeAngle(o.getLM(), FastMath.PI)); // Cartesian elements // Position along X in inertial frame (km) final double px = o.getPVCoordinates().getPosition().getX() / 1000.; // Position along Y in inertial frame (km) final double py = o.getPVCoordinates().getPosition().getY() / 1000.; // Position along Z in inertial frame (km) final double pz = o.getPVCoordinates().getPosition().getZ() / 1000.; // Velocity along X in inertial frame (km/s) final double vx = o.getPVCoordinates().getVelocity().getX() / 1000.; // Velocity along Y in inertial frame (km/s) final double vy = o.getPVCoordinates().getVelocity().getY() / 1000.; // Velocity along Z in inertial frame (km/s) final double vz = o.getPVCoordinates().getVelocity().getZ() / 1000.; buffer.write( f.format(format, time, a, e, i, ra, pa, am, h, k, p, q, lm, px, py, pz, vx, vy, vz).toString()); buffer.newLine(); f.close(); } buffer.close(); } /** Specialized step handler catching the orbit at each step. */ private static class OrbitHandler implements OrekitFixedStepHandler { /** List of orbits. */ private final List<Orbit> orbits; private OrbitHandler() { // initialise an empty list of orbit orbits = new ArrayList<Orbit>(); } /** {@inheritDoc} */ public void init(final SpacecraftState s0, final AbsoluteDate t) { } /** {@inheritDoc} */ public void handleStep(SpacecraftState currentState, boolean isLast) { // fill in the list with the orbit from the current step orbits.add(currentState.getOrbit()); } /** Get the list of propagated orbits. * @return orbits */ public List<Orbit> getOrbits() { return orbits; } } }