DSSTPropagator.java

  1. /* Copyright 2002-2013 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.propagation.semianalytical.dsst;

  19. import java.io.NotSerializableException;
  20. import java.io.Serializable;
  21. import java.util.ArrayList;
  22. import java.util.Arrays;
  23. import java.util.List;

  25. import org.apache.commons.math3.ode.AbstractIntegrator;
  26. import org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator;
  27. import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
  28. import org.orekit.attitudes.Attitude;
  29. import org.orekit.attitudes.AttitudeProvider;
  30. import org.orekit.bodies.CelestialBody;
  31. import org.orekit.bodies.CelestialBodyFactory;
  32. import org.orekit.errors.OrekitException;
  33. import org.orekit.errors.OrekitMessages;
  34. import org.orekit.errors.PropagationException;
  35. import org.orekit.forces.ForceModel;
  36. import org.orekit.forces.SphericalSpacecraft;
  37. import org.orekit.forces.drag.Atmosphere;
  38. import org.orekit.forces.drag.DragForce;
  39. import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
  40. import org.orekit.forces.gravity.ThirdBodyAttraction;
  41. import org.orekit.forces.gravity.potential.GravityFieldFactory;
  42. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  43. import org.orekit.forces.radiation.SolarRadiationPressure;
  44. import org.orekit.frames.Frame;
  45. import org.orekit.frames.FramesFactory;
  46. import org.orekit.orbits.EquinoctialOrbit;
  47. import org.orekit.orbits.Orbit;
  48. import org.orekit.orbits.OrbitType;
  49. import org.orekit.orbits.PositionAngle;
  50. import org.orekit.propagation.Propagator;
  51. import org.orekit.propagation.SpacecraftState;
  52. import org.orekit.propagation.conversion.OsculatingToMeanElementsConverter;
  53. import org.orekit.propagation.events.EventDetector;
  54. import org.orekit.propagation.integration.AbstractIntegratedPropagator;
  55. import org.orekit.propagation.integration.StateMapper;
  56. import org.orekit.propagation.numerical.NumericalPropagator;
  57. import org.orekit.propagation.semianalytical.dsst.forces.DSSTAtmosphericDrag;
  58. import org.orekit.propagation.semianalytical.dsst.forces.DSSTCentralBody;
  59. import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
  60. import org.orekit.propagation.semianalytical.dsst.forces.DSSTSolarRadiationPressure;
  61. import org.orekit.propagation.semianalytical.dsst.forces.DSSTThirdBody;
  62. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  63. import org.orekit.time.AbsoluteDate;
  64. import org.orekit.utils.IERSConventions;

  66. /**
  67.  * This class propagates {@link org.orekit.orbits.Orbit orbits} using the DSST theory.
  68.  * <p>
  69.  * Whereas analytical propagators are configured only thanks to their various
  70.  * constructors and can be used immediately after construction, such a semianalytical
  71.  * propagator configuration involves setting several parameters between construction
  72.  * time and propagation time, just as numerical propagators.
  73.  * </p>
  74.  * <p>
  75.  * The configuration parameters that can be set are:
  76.  * </p>
  77.  * <ul>
  78.  * <li>the initial spacecraft state ({@link #setInitialState(SpacecraftState)})</li>
  79.  * <li>the various force models ({@link #addForceModel(DSSTForceModel)},
  80.  * {@link #removeForceModels()})</li>
  81.  * <li>the discrete events that should be triggered during propagation (
  82.  * {@link #addEventDetector(org.orekit.propagation.events.EventDetector)},
  83.  * {@link #clearEventsDetectors()})</li>
  84.  * <li>the binding logic with the rest of the application ({@link #setSlaveMode()},
  85.  * {@link #setMasterMode(double, org.orekit.propagation.sampling.OrekitFixedStepHandler)},
  86.  * {@link #setMasterMode(org.orekit.propagation.sampling.OrekitStepHandler)},
  87.  * {@link #setEphemerisMode()}, {@link #getGeneratedEphemeris()})</li>
  88.  * </ul>
  89.  * <p>
  90.  * From these configuration parameters, only the initial state is mandatory.
  91.  * The default propagation settings are in {@link OrbitType#EQUINOCTIAL equinoctial}
  92.  * parameters with {@link PositionAngle#TRUE true} longitude argument.
  93.  * The central attraction coefficient used to define the initial orbit will be used.
  94.  * However, specifying only the initial state would mean the propagator would use
  95.  * only keplerian forces. In this case, the simpler
  96.  * {@link org.orekit.propagation.analytical.KeplerianPropagator KeplerianPropagator}
  97.  * class would be more effective.
  98.  * </p>
  99.  * <p>
  100.  * The underlying numerical integrator set up in the constructor may also have
  101.  * its own configuration parameters. Typical configuration parameters for adaptive
  102.  * stepsize integrators are the min, max and perhaps start step size as well as
  103.  * the absolute and/or relative errors thresholds.
  104.  * </p>
  105.  * <p>
  106.  * The state that is seen by the integrator is a simple six elements double array.
  107.  * These six elements are:
  108.  * <ul>
  109.  * <li>the {@link org.orekit.orbits.EquinoctialOrbit equinoctial orbit parameters}
  110.  * (a, e<sub>x</sub>, e<sub>y</sub>, h<sub>x</sub>, h<sub>y</sub>, &lambda;<sub>m</sub>)
  111.  * in meters and radians,</li>
  112.  * </ul>
  113.  * </p>
  114.  * <p>
  115.  * The same propagator can be reused for several orbit extrapolations,
  116.  * by resetting the initial state without modifying the other configuration
  117.  * parameters. However, the same instance cannot be used simultaneously by
  118.  * different threads, the class is <em>not</em> thread-safe.
  119.  * </p>
  120.  *
  121.  * @see SpacecraftState
  122.  * @see DSSTForceModel
  123.  * @author Romain Di Costanzo
  124.  * @author Pascal Parraud
  125.  */
  126. public class DSSTPropagator extends AbstractIntegratedPropagator {

  128.     /** Retrograde factor I.
  129.      *  <p>
  130.      *  DSST model needs equinoctial orbit as internal representation.
  131.      *  Classical equinoctial elements have discontinuities when inclination
  132.      *  is close to zero. In this representation, I = +1. <br>
  133.      *  To avoid this discontinuity, another representation exists and equinoctial
  134.      *  elements can be expressed in a different way, called "retrograde" orbit.
  135.      *  This implies I = -1. <br>
  136.      *  As Orekit doesn't implement the retrograde orbit, I is always set to +1.
  137.      *  But for the sake of consistency with the theory, the retrograde factor
  138.      *  has been kept in the formulas.
  139.      *  </p>
  140.      */
  141.     private static final int I = 1;

  143.     /** State mapper holding the force models. */
  144.     private MeanPlusShortPeriodicMapper mapper;

  146.     /** Create a new instance of DSSTPropagator.
  147.      *  <p>
  148.      *  After creation, there are no perturbing forces at all.
  149.      *  This means that if {@link #addForceModel addForceModel}
  150.      *  is not called after creation, the integrated orbit will
  151.      *  follow a keplerian evolution only.
  152.      *  </p>
  153.      *  @param integrator numerical integrator to use for propagation.
  154.      */
  155.     public DSSTPropagator(final AbstractIntegrator integrator) {
  156.         super(integrator);
  157.         initMapper();
  158.         // DSST uses only equinoctial orbits and mean longitude argument
  159.         setOrbitType(OrbitType.EQUINOCTIAL);
  160.         setPositionAngleType(PositionAngle.TRUE);
  161.         setAttitudeProvider(DEFAULT_LAW);
  162.     }

  164.     /** Set the initial state with osculating orbital elements.
  165.      *  @param initialState initial state (defined with osculating elements)
  166.      *  @throws PropagationException if the initial state cannot be set
  167.      */
  168.     public void setInitialState(final SpacecraftState initialState)
  169.         throws PropagationException {
  170.         setInitialState(initialState, true);
  171.     }

  173.     /** Set the initial state.
  174.      *  @param initialState initial state
  175.      *  @param withOsculatingElements true if the orbital state is defined with osculating elements
  176.      *  @throws PropagationException if the initial state cannot be set
  177.      */
  178.     public void setInitialState(final SpacecraftState initialState,
  179.                                 final boolean withOsculatingElements)
  180.         throws PropagationException {
  181.         mapper.setInitialIsOsculating(withOsculatingElements);
  182.         resetInitialState(initialState);
  183.     }

  185.     /** Reset the initial state.
  186.      *
  187.      *  @param state new initial state
  188.      *  @throws PropagationException if initial state cannot be reset
  189.      */
  190.     public void resetInitialState(final SpacecraftState state) throws PropagationException {
  191.         super.setStartDate(state.getDate());
  192.         super.resetInitialState(state);
  193.     }

  195.     /** Check if the initial state is provided in osculating elements.
  196.      * @return true if initial state is provided in osculating elements
  197.      */
  198.     public boolean initialIsOsculating() {
  199.         return mapper.initialIsOsculating();
  200.     }

  202.     /** Add a force model to the global perturbation model.
  203.      *  <p>
  204.      *  If this method is not called at all,
  205.      *  the integrated orbit will follow a keplerian evolution only.
  206.      *  </p>
  207.      *  @param force perturbing {@link DSSTForceModel force} to add
  208.      *  @see #removeForceModels()
  209.      */
  210.     public void addForceModel(final DSSTForceModel force) {
  211.         mapper.addForceModel(force);
  212.     }

  214.     /** Remove all perturbing force models from the global perturbation model.
  215.      *  <p>
  216.      *  Once all perturbing forces have been removed (and as long as no new force model is added),
  217.      *  the integrated orbit will follow a keplerian evolution only.
  218.      *  </p>
  219.      *  @see #addForceModel(DSSTForceModel)
  220.      */
  221.     public void removeForceModels() {
  222.         mapper.removeForceModels();
  223.     }

  225.     /** Override the default value of the parameter.
  226.      *  <p>
  227.      *  By default, if the initial orbit is defined as osculating,
  228.      *  it will be averaged over 2 satellite revolutions.
  229.      *  This can be changed by using this method.
  230.      *  </p>
  231.      *  @param satelliteRevolution number of satellite revolutions to use for converting osculating to mean
  232.      *                             elements
  233.      */
  234.     public void setSatelliteRevolution(final int satelliteRevolution) {
  235.         mapper.setSatelliteRevolution(satelliteRevolution);
  236.     }

  238.     /** Get the number of satellite revolutions to use for converting osculating to mean elements.
  239.      *  @return number of satellite revolutions to use for converting osculating to mean elements
  240.      */
  241.     public int getSatelliteRevolution() {
  242.         return mapper.getSatelliteRevolution();
  243.     }

  245.     /** Method called just before integration.
  246.      * <p>
  247.      * The default implementation does nothing, it may be specialized in subclasses.
  248.      * </p>
  249.      * @param initialState initial state
  250.      * @param tEnd target date at which state should be propagated
  251.      * @exception OrekitException if hook cannot be run
  252.      */
  253.     protected void beforeIntegration(final SpacecraftState initialState,
  254.                                      final AbsoluteDate tEnd)
  255.         throws OrekitException {

  257.         // compute common auxiliary elements
  258.         final AuxiliaryElements aux = new AuxiliaryElements(initialState.getOrbit(), I);

  260.         // initialize all perturbing forces
  261.         for (final DSSTForceModel force : mapper.getForceModels()) {
  262.             force.initialize(aux);
  263.         }

  265.     }

  267.     /** {@inheritDoc} */
  268.     protected StateMapper createMapper(final AbsoluteDate referenceDate, final double mu,
  269.                                        final OrbitType orbitType, final PositionAngle positionAngleType,
  270.                                        final AttitudeProvider attitudeProvider, final Frame frame) {

  272.         // create a mapper with the common settings provided as arguments
  273.         final MeanPlusShortPeriodicMapper newMapper =
  274.                 new MeanPlusShortPeriodicMapper(referenceDate, mu, attitudeProvider, frame);

  276.         // copy the specific settings from the existing mapper
  277.         if (mapper != null) {
  278.             for (final DSSTForceModel forceModel : mapper.getForceModels()) {
  279.                 newMapper.addForceModel(forceModel);
  280.             }
  281.             newMapper.setSatelliteRevolution(mapper.getSatelliteRevolution());
  282.             newMapper.setInitialIsOsculating(mapper.initialIsOsculating());
  283.         }

  285.         mapper = newMapper;
  286.         return mapper;

  288.     }

  290.     /** Internal mapper using mean parameters plus short periodic terms. */
  291.     private static class MeanPlusShortPeriodicMapper extends StateMapper implements Serializable {

  293.         /** Serializable UID. */
  294.         private static final long serialVersionUID = 20130621L;

  296.         /** Flag specifying whether the initial orbital state is given with osculating elements. */
  297.         private boolean                    initialIsOsculating;

  299.         /** Force models used to compute short periodic terms. */
  300.         private final transient List<DSSTForceModel> forceModels;

  302.         /** Number of satellite revolutions in the averaging interval. */
  303.         private int                        satelliteRevolution;

  305.         /** Simple constructor.
  306.          * <p>
  307.          * The position parameter type is meaningful only if {@link
  308.          * #getOrbitType() propagation orbit type}
  309.          * support it. As an example, it is not meaningful for propagation
  310.          * in {@link OrbitType#CARTESIAN Cartesian} parameters.
  311.          * </p>
  312.          * @param referenceDate reference date
  313.          * @param mu central attraction coefficient (m<sup>3</sup>/s<sup>2</sup>)
  314.          * @param attitudeProvider attitude provider
  315.          * @param frame inertial frame
  316.          */
  317.         public MeanPlusShortPeriodicMapper(final AbsoluteDate referenceDate, final double mu,
  318.                                            final AttitudeProvider attitudeProvider, final Frame frame) {

  320.             super(referenceDate, mu, OrbitType.EQUINOCTIAL, PositionAngle.MEAN, attitudeProvider, frame);

  322.             this.forceModels = new ArrayList<DSSTForceModel>();

  324.             // Default averaging period for conversion from osculating to mean elements
  325.             this.satelliteRevolution = 2;

  327.             this.initialIsOsculating = true;

  329.         }

  331.         /** {@inheritDoc} */
  332.         public SpacecraftState mapArrayToState(final double t, final double[] y)
  333.             throws OrekitException {

  335.             final AbsoluteDate date = mapDoubleToDate(t);

  337.             // add short periodic variations to mean elements to get osculating elements
  338.             // (the loop may not be performed if there are no force models, in the
  339.             //  case we want to remain in mean parameters only)
  340.             final double[] osculatingElements = y.clone();
  341.             for (final DSSTForceModel forceModel : forceModels) {
  342.                 final double[] shortPeriodic = forceModel.getShortPeriodicVariations(date, y);
  343.                 for (int i = 0; i < shortPeriodic.length; i++) {
  344.                     osculatingElements[i] += shortPeriodic[i];
  345.                 }
  346.             }

  348.             final double mass = y[6];
  349.             if (mass <= 0.0) {
  350.                 throw new PropagationException(OrekitMessages.SPACECRAFT_MASS_BECOMES_NEGATIVE, mass);
  351.             }

  353.             final Orbit orbit       = OrbitType.EQUINOCTIAL.mapArrayToOrbit(y, PositionAngle.MEAN, date, getMu(), getFrame());
  354.             final Attitude attitude = getAttitudeProvider().getAttitude(orbit, date, getFrame());

  356.             return new SpacecraftState(orbit, attitude, mass);

  358.         }

  360.         /** {@inheritDoc} */
  361.         public void mapStateToArray(final SpacecraftState state, final double[] y)
  362.             throws OrekitException {

  364.             final Orbit meanOrbit;
  365.             if (!initialIsOsculating) {
  366.                 // the state is considered to be already a mean state
  367.                 meanOrbit = state.getOrbit();
  368.             } else {
  369.                 // the state is considered to be an osculating state
  370.                 final Propagator propagator = createPropagator(state);
  371.                 meanOrbit = new OsculatingToMeanElementsConverter(state, satelliteRevolution, propagator).convert().getOrbit();
  372.             }

  374.             OrbitType.EQUINOCTIAL.mapOrbitToArray(meanOrbit, PositionAngle.MEAN, y);
  375.             y[6] = state.getMass();

  377.         }

  379.         /** Add a force model to the global perturbation model.
  380.          *  <p>
  381.          *  If this method is not called at all,
  382.          *  the integrated orbit will follow a keplerian evolution only.
  383.          *  </p>
  384.          *  @param force perturbing {@link DSSTForceModel force} to add
  385.          *  @see #removeForceModels()
  386.          */
  387.         public void addForceModel(final DSSTForceModel force) {
  388.             forceModels.add(force);
  389.         }

  391.         /** Remove all perturbing force models from the global perturbation model.
  392.          *  <p>
  393.          *  Once all perturbing forces have been removed (and as long as no new force model is added),
  394.          *  the integrated orbit will follow a keplerian evolution only.
  395.          *  </p>
  396.          *  @see #addForceModel(DSSTForceModel)
  397.          */
  398.         public void removeForceModels() {
  399.             forceModels.clear();
  400.         }

  402.         /** Get the force models.
  403.          * @return force models
  404.          */
  405.         public List<DSSTForceModel> getForceModels() {
  406.             return forceModels;
  407.         }

  409.         /** Set the number of satellite revolutions to use for converting osculating to mean elements.
  410.          *  <p>
  411.          *  By default, if the initial orbit is defined as osculating,
  412.          *  it will be averaged over 2 satellite revolutions.
  413.          *  This can be changed by using this method.
  414.          *  </p>
  415.          *  @param satelliteRevolution number of satellite revolutions to use for converting osculating to mean
  416.          *                             elements
  417.          */
  418.         public void setSatelliteRevolution(final int satelliteRevolution) {
  419.             this.satelliteRevolution = satelliteRevolution;
  420.         }

  422.         /** Get the number of satellite revolutions to use for converting osculating to mean elements.
  423.          *  @return number of satellite revolutions to use for converting osculating to mean elements
  424.          */
  425.         public int getSatelliteRevolution() {
  426.             return satelliteRevolution;
  427.         }

  429.         /** Set the osculating parameters flag.
  430.          * @param initialIsOsculating if true, the initial state is provided in osculating elements
  431.          */
  432.         public void setInitialIsOsculating(final boolean initialIsOsculating) {
  433.             this.initialIsOsculating = initialIsOsculating;
  434.         }

  436.         /** Check if the initial state is provided in osculating elements.
  437.          * @return true if initial state is provided in osculating elements
  438.          */
  439.         public boolean initialIsOsculating() {
  440.             return initialIsOsculating;
  441.         }

  443.         /** Create a reference numerical propagator to convert orbit to mean elements.
  444.          *  @param initialState initial state
  445.          *  @return propagator
  446.          *  @throws OrekitException if some numerical force model cannot be built
  447.          */
  448.         private Propagator createPropagator(final SpacecraftState initialState)
  449.             throws OrekitException {
  450.             final Orbit initialOrbit = initialState.getOrbit();
  451.             final double[][] tol = NumericalPropagator.tolerances(1.0, initialOrbit, OrbitType.EQUINOCTIAL);
  452.             final double minStep = 1.;
  453.             final double maxStep = 200.;
  454.             final AdaptiveStepsizeIntegrator integ = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
  455.             integ.setInitialStepSize(100.);

  457.             final NumericalPropagator propagator = new NumericalPropagator(integ);
  458.             propagator.setOrbitType(OrbitType.EQUINOCTIAL);
  459.             propagator.setInitialState(initialState);

  461.             // Define the same force model as the DSST
  462.             for (final DSSTForceModel force : forceModels) {
  463.                 if (force instanceof DSSTCentralBody) {
  464.                     // Central body
  465.                     final UnnormalizedSphericalHarmonicsProvider provider = ((DSSTCentralBody) force).getProvider();
  466.                     final ForceModel holmesFeatherstone =
  467.                             new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true),
  468.                                                                   GravityFieldFactory.getNormalizedProvider(provider));
  469.                     propagator.addForceModel(holmesFeatherstone);
  470.                 } else if (force instanceof DSSTThirdBody) {
  471.                     // Third body
  472.                     final CelestialBody body = ((DSSTThirdBody) force).getBody();
  473.                     final ForceModel third   = new ThirdBodyAttraction(body);
  474.                     propagator.addForceModel(third);
  475.                 } else if (force instanceof DSSTAtmosphericDrag) {
  476.                     // Atmospheric drag
  477.                     final Atmosphere atm = ((DSSTAtmosphericDrag) force).getAtmosphere();
  478.                     final double area    = ((DSSTAtmosphericDrag) force).getArea();
  479.                     final double cd      = ((DSSTAtmosphericDrag) force).getCd();
  480.                     final SphericalSpacecraft scr = new SphericalSpacecraft(area, cd, 0., 0.);
  481.                     final ForceModel drag = new DragForce(atm, scr);
  482.                     propagator.addForceModel(drag);
  483.                 } else if (force instanceof DSSTSolarRadiationPressure) {
  484.                     // Solar radiation pressure
  485.                     final double ae   = ((DSSTSolarRadiationPressure) force).getEquatorialRadius();
  486.                     final double area = ((DSSTSolarRadiationPressure) force).getArea();
  487.                     final double cr   = ((DSSTSolarRadiationPressure) force).getCr();
  488.                     // Convert DSST SRP coefficient convention to numerical's one
  489.                     final double kr   = 3.25 - 2.25 * cr;
  490.                     final SphericalSpacecraft scr = new SphericalSpacecraft(area, 0., 0., kr);
  491.                     final ForceModel pressure = new SolarRadiationPressure(CelestialBodyFactory.getSun(), ae, scr);
  492.                     propagator.addForceModel(pressure);
  493.                 }
  494.             }
  495.             return propagator;
  496.         }

  498.         /** Replace the instance with a data transfer object for serialization.
  499.          * @return data transfer object that will be serialized
  500.          * @exception NotSerializableException if one of the force models cannot be serialized
  501.          */
  502.         private Object writeReplace() throws NotSerializableException {

  504.             // Check the force models can be serialized
  505.             final DSSTForceModel[] serializableorceModels = new DSSTForceModel[forceModels.size()];
  506.             for (int i = 0; i < serializableorceModels.length; ++i) {
  507.                 final DSSTForceModel forceModel = forceModels.get(i);
  508.                 if (forceModel instanceof Serializable) {
  509.                     serializableorceModels[i] = forceModel;
  510.                 } else {
  511.                     throw new NotSerializableException(forceModel.getClass().getName());
  512.                 }
  513.             }
  514.             return new DataTransferObject(getReferenceDate(), getMu(), getAttitudeProvider(), getFrame(),
  515.                                           initialIsOsculating, serializableorceModels, satelliteRevolution);
  516.         }

  518.         /** Internal class used only for serialization. */
  519.         private static class DataTransferObject implements Serializable {

  521.             /** Serializable UID. */
  522.             private static final long serialVersionUID = 20130621L;

  524.             /** Reference date. */
  525.             private final AbsoluteDate referenceDate;

  527.             /** Central attraction coefficient (m<sup>3</sup>/s<sup>2</sup>). */
  528.             private final double mu;

  530.             /** Attitude provider. */
  531.             private final AttitudeProvider attitudeProvider;

  533.             /** Inertial frame. */
  534.             private final Frame frame;

  536.             /** Flag specifying whether the initial orbital state is given with osculating elements. */
  537.             private final boolean initialIsOsculating;

  539.             /** Force models to use for short periodic terms computation. */
  540.             private final DSSTForceModel[] forceModels;

  542.             /** Number of satellite revolutions in the averaging interval. */
  543.             private final int satelliteRevolution;

  545.             /** Simple constructor.
  546.              * @param referenceDate reference date
  547.              * @param mu central attraction coefficient (m<sup>3</sup>/s<sup>2</sup>)
  548.              * @param attitudeProvider attitude provider
  549.              * @param frame inertial frame
  550.              * @param initialIsOsculating if true, initial orbital state is given with osculating elements
  551.              * @param forceModels force models to use for short periodic terms computation
  552.              * @param satelliteRevolution number of satellite revolutions in the averaging interval
  553.              */
  554.             public DataTransferObject(final AbsoluteDate referenceDate, final double mu,
  555.                                       final AttitudeProvider attitudeProvider, final Frame frame,
  556.                                       final boolean initialIsOsculating,
  557.                                       final DSSTForceModel[] forceModels, final int satelliteRevolution) {
  558.                 this.referenceDate       = referenceDate;
  559.                 this.mu                  = mu;
  560.                 this.attitudeProvider    = attitudeProvider;
  561.                 this.frame               = frame;
  562.                 this.initialIsOsculating = initialIsOsculating;
  563.                 this.forceModels         = forceModels;
  564.                 this.satelliteRevolution = satelliteRevolution;
  565.             }

  567.             /** Replace the deserialized data transfer object with a {@link MeanPlusShortPeriodicMapper}.
  568.              * @return replacement {@link MeanPlusShortPeriodicMapper}
  569.              */
  570.             private Object readResolve() {
  571.                 final MeanPlusShortPeriodicMapper mapper =
  572.                         new MeanPlusShortPeriodicMapper(referenceDate, mu, attitudeProvider, frame);
  573.                 for (final DSSTForceModel forceModel : forceModels) {
  574.                     mapper.addForceModel(forceModel);
  575.                 }
  576.                 mapper.setSatelliteRevolution(satelliteRevolution);
  577.                 mapper.setInitialIsOsculating(initialIsOsculating);
  578.                 return mapper;
  579.             }

  581.         }

  583.     }

  585.     /** {@inheritDoc} */
  586.     protected MainStateEquations getMainStateEquations(final AbstractIntegrator integrator) {
  587.         return new Main(integrator);
  588.     }

  590.     /** Internal class for mean parameters integration. */
  591.     private class Main implements MainStateEquations {

  593.         /** Derivatives array. */
  594.         private final double[] yDot;

  596.         /** Simple constructor.
  597.          * @param integrator numerical integrator to use for propagation.
  598.          */
  599.         public Main(final AbstractIntegrator integrator) {
  600.             yDot = new double[7];

  602.             for (final DSSTForceModel forceModel : mapper.getForceModels()) {
  603.                 final EventDetector[] modelDetectors = forceModel.getEventsDetectors();
  604.                 if (modelDetectors != null) {
  605.                     for (final EventDetector detector : modelDetectors) {
  606.                         setUpEventDetector(integrator, detector);
  607.                     }
  608.                 }
  609.             }

  611.         }

  613.         /** {@inheritDoc} */
  614.         public double[] computeDerivatives(final SpacecraftState state) throws OrekitException {

  616.             // compute common auxiliary elements
  617.             final AuxiliaryElements aux = new AuxiliaryElements(state.getOrbit(), I);

  619.             // initialize all perturbing forces
  620.             for (final DSSTForceModel force : mapper.getForceModels()) {
  621.                 force.initializeStep(aux);
  622.             }

  624.             Arrays.fill(yDot, 0.0);

  626.             // compute the contributions of all perturbing forces
  627.             for (final DSSTForceModel forceModel : mapper.getForceModels()) {
  628.                 final double[] daidt = forceModel.getMeanElementRate(state);
  629.                 for (int i = 0; i < daidt.length; i++) {
  630.                     yDot[i] += daidt[i];
  631.                 }
  632.             }

  634.             // finalize derivatives by adding the Kepler contribution
  635.             final EquinoctialOrbit orbit = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(state.getOrbit());
  636.             orbit.addKeplerContribution(PositionAngle.MEAN, getMu(), yDot);

  638.             return yDot.clone();
  639.         }

  641.     }

  643.     /** Estimate tolerance vectors for an AdaptativeStepsizeIntegrator.
  644.      *  <p>
  645.      *  The errors are estimated from partial derivatives properties of orbits,
  646.      *  starting from a scalar position error specified by the user.
  647.      *  Considering the energy conservation equation V = sqrt(mu (2/r - 1/a)),
  648.      *  we get at constant energy (i.e. on a Keplerian trajectory):
  649.      *
  650.      *  <pre>
  651.      *  V<sup>2</sup> r |dV| = mu |dr|
  652.      *  </pre>
  653.      *
  654.      *  So we deduce a scalar velocity error consistent with the position error. From here, we apply
  655.      *  orbits Jacobians matrices to get consistent errors on orbital parameters.
  656.      *  </p>
  657.      *  <p>
  658.      *  The tolerances are only <em>orders of magnitude</em>, and integrator tolerances are only
  659.      *  local estimates, not global ones. So some care must be taken when using these tolerances.
  660.      *  Setting 1mm as a position error does NOT mean the tolerances will guarantee a 1mm error
  661.      *  position after several orbits integration.
  662.      *  </p>
  663.      *
  664.      * @param dP user specified position error (m)
  665.      * @param orbit reference orbit
  666.      * @return a two rows array, row 0 being the absolute tolerance error
  667.      *                       and row 1 being the relative tolerance error
  668.      * @exception PropagationException if Jacobian is singular
  669.      */
  670.     public static double[][] tolerances(final double dP,
  671.                                         final Orbit orbit)
  672.         throws PropagationException {

  674.         return NumericalPropagator.tolerances(dP, orbit, OrbitType.EQUINOCTIAL);

  676.     }

  678. }
Created with JaCoCo 0.6.3.201306030806