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.conversion; import java.io.File; import java.io.FileNotFoundException; import java.io.PrintStream; import java.util.ArrayList; import java.util.List; import org.apache.commons.math3.ode.AbstractIntegrator; 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.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.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.analytical.KeplerianPropagator; import org.orekit.propagation.conversion.FiniteDifferencePropagatorConverter; import org.orekit.propagation.conversion.KeplerianPropagatorBuilder; import org.orekit.propagation.conversion.PropagatorBuilder; import org.orekit.propagation.conversion.PropagatorConverter; import org.orekit.propagation.numerical.NumericalPropagator; import org.orekit.propagation.sampling.OrekitFixedStepHandler; import org.orekit.time.AbsoluteDate; import org.orekit.time.TimeScalesFactory; import org.orekit.utils.IERSConventions; import fr.cs.examples.Autoconfiguration; /** Orekit tutorial for propagator conversion. * @author Pascal Parraud */ public class PropagatorConversion { /** Program entry point. * @param args program arguments (unused here) */ public static void main(String[] args) { try { // configure Orekit Autoconfiguration.configureOrekit(); // gravity field NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(2, 0); double mu = provider.getMu(); // inertial frame Frame inertialFrame = FramesFactory.getEME2000(); // Initial date AbsoluteDate initialDate = new AbsoluteDate(2004, 01, 01, 23, 30, 00.000, TimeScalesFactory.getUTC()); // Initial orbit (GTO) final double a = 24396159; // semi major axis in meters final double e = 0.72831215; // eccentricity final double i = FastMath.toRadians(7); // inclination final double omega = FastMath.toRadians(180); // perigee argument final double raan = FastMath.toRadians(261); // right ascention of ascending node final double lM = 0; // mean anomaly Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.MEAN, inertialFrame, initialDate, mu); final double period = initialOrbit.getKeplerianPeriod(); // Initial state definition final 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 = 1000.; final double dP = 1.e-2; final OrbitType orbType = OrbitType.CARTESIAN; final double[][] tol = NumericalPropagator.tolerances(dP, initialOrbit, orbType); final AbstractIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]); // Propagator NumericalPropagator numProp = new NumericalPropagator(integrator); numProp.setInitialState(initialState); numProp.setOrbitType(orbType); // Force Models: // 1 - Perturbing gravity field (only J2 is considered here) ForceModel gravity = new HolmesFeatherstoneAttractionModel( FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider); // Add force models to the propagator numProp.addForceModel(gravity); // Propagator factory PropagatorBuilder builder = new KeplerianPropagatorBuilder(mu, inertialFrame, OrbitType.KEPLERIAN, PositionAngle.TRUE); // Propagator converter PropagatorConverter fitter = new FiniteDifferencePropagatorConverter(builder, 1.e-6, 5000); // Resulting propagator KeplerianPropagator kepProp = (KeplerianPropagator) fitter.convert(numProp, 2 * period, 251); // Step handlers StatesHandler numStepHandler = new StatesHandler(); StatesHandler kepStepHandler = new StatesHandler(); // Set up operating mode for the propagator as master mode // with fixed step and specialized step handler numProp.setMasterMode(60., numStepHandler); kepProp.setMasterMode(60., kepStepHandler); // Extrapolate from the initial to the final date numProp.propagate(initialDate.shiftedBy(10. * period)); kepProp.propagate(initialDate.shiftedBy(10. * period)); // retrieve the states List<SpacecraftState> numStates = numStepHandler.getStates(); List<SpacecraftState> kepStates = kepStepHandler.getStates(); // Print the results on the output file File output = new File(new File(System.getProperty("user.home")), "elements.dat"); PrintStream stream = new PrintStream(output); stream.println("# date Anum Akep Enum Ekep Inum Ikep LMnum LMkep"); for (SpacecraftState numState : numStates) { for (SpacecraftState kepState : kepStates) { if (numState.getDate().compareTo(kepState.getDate()) == 0) { stream.println(numState.getDate() + " " + numState.getA() + " " + kepState.getA() + " " + numState.getE() + " " + kepState.getE() + " " + FastMath.toDegrees(numState.getI()) + " " + FastMath.toDegrees(kepState.getI()) + " " + FastMath.toDegrees(MathUtils.normalizeAngle(numState.getLM(), FastMath.PI)) + " " + FastMath.toDegrees(MathUtils.normalizeAngle(kepState.getLM(), FastMath.PI))); break; } } } stream.close(); System.out.println("Results saved as file " + output); File output1 = new File(new File(System.getProperty("user.home")), "elts_pv.dat"); PrintStream stream1 = new PrintStream(output1); stream.println("# date pxn pyn pzn vxn vyn vzn pxk pyk pzk vxk vyk vzk"); for (SpacecraftState numState : numStates) { for (SpacecraftState kepState : kepStates) { if (numState.getDate().compareTo(kepState.getDate()) == 0) { final double pxn = numState.getPVCoordinates().getPosition().getX(); final double pyn = numState.getPVCoordinates().getPosition().getY(); final double pzn = numState.getPVCoordinates().getPosition().getZ(); final double vxn = numState.getPVCoordinates().getVelocity().getX(); final double vyn = numState.getPVCoordinates().getVelocity().getY(); final double vzn = numState.getPVCoordinates().getVelocity().getZ(); final double pxk = kepState.getPVCoordinates().getPosition().getX(); final double pyk = kepState.getPVCoordinates().getPosition().getY(); final double pzk = kepState.getPVCoordinates().getPosition().getZ(); final double vxk = kepState.getPVCoordinates().getVelocity().getX(); final double vyk = kepState.getPVCoordinates().getVelocity().getY(); final double vzk = kepState.getPVCoordinates().getVelocity().getZ(); stream1.println(numState.getDate() + " " + pxn + " " + pyn + " " + pzn + " " + vxn + " " + vyn + " " + vzn + " " + pxk + " " + pyk + " " + pzk + " " + vxk + " " + vyk + " " + vzk); break; } } } stream1.close(); System.out.println("Results saved as file " + output1); } catch (OrekitException oe) { System.err.println(oe.getLocalizedMessage()); System.exit(1); } catch (FileNotFoundException fnfe) { System.err.println(fnfe.getLocalizedMessage()); System.exit(1); } } /** Specialized step handler. * <p>This class extends the step handler in order to handle states at each step.<p> */ private static class StatesHandler implements OrekitFixedStepHandler { /** Points. */ private final List<SpacecraftState> states; public StatesHandler() { // prepare an empty list of SpacecraftState states = new ArrayList<SpacecraftState>(); } public void init(final SpacecraftState s0, final AbsoluteDate t) { } public void handleStep(SpacecraftState currentState, boolean isLast) { // add the current state states.add(currentState); } /** Get the list of handled orbital elements. * @return orbital elements list */ public List<SpacecraftState> getStates() { return states; } } }