Added an analytical model for small maneuvers.

The model all allows to compute at date t1 the effect of a
maneuver performed at date t0, without relying on any propagation (the
model applies directly to SpacecraftState instances).

src/main/java/org/orekit/forces/maneuvers/SmallManeuverAnalyticalModel.java

+/* Copyright 2002-2011 CS Communication & Systèmes
+ * Licensed to CS Communication & Systèmes (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 org.orekit.forces.maneuvers;
+
+import java.io.Serializable;
+
+import org.apache.commons.math.geometry.euclidean.threed.Vector3D;
+import org.apache.commons.math.util.FastMath;
+import org.orekit.errors.OrekitException;
+import org.orekit.frames.Frame;
+import org.orekit.orbits.Orbit;
+import org.orekit.orbits.OrbitType;
+import org.orekit.orbits.PositionAngle;
+import org.orekit.propagation.SpacecraftState;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.utils.Constants;
+
+/** Analytical model for small maneuvers.
+ * <p>The aim of this model is to compute quickly the effect at date t<sub>1</sub>
+ * of a small maneuver performed at an earlier date t<sub>0</sub>. This effect is
+ * computed analytically using two Jacobian matrices:
+ * <ol>
+ *   <li>J<sub>0</sub>: Jacobian of Keplerian or equinoctial elements with respect
+ *   to cartesian parameters at date t<sub>0</sub></li> allows to compute
+ *   maneuver effect as an orbital elements change,
+ *   <li>J<sub>1/0</sub></li>: Jacobian of Keplerian or equinoctial elements
+ *   at date t<sub>1</sub> with respect to Keplerian or equinoctial elements
+ *   at date t<sub>0</sub> allows to propagate the orbital elements change to
+ *   final date.
+ * </ol>
+ * </p>
+ * <p>
+ * The second Jacobian, J<sub>1/0</sub></li>, is computed using a simple Keplerian
+ * model, i.e. it is the identity except for the mean motion row which also includes
+ * an off-diagonal element due to semi-major axis change.
+ * </p>
+ * <p>
+ * The orbital elements change at date t<sub>1</sub> can be added to orbital elements
+ * extracted from state, and the final elements taking account the changes are then
+ * converted back to appropriate type, which may be different from Keplerian or
+ * equinoctial elements.
+ * </p>
+ * @author Luc Maisonobe
+ */
+public class SmallManeuverAnalyticalModel implements Serializable {
+
+    /** Serializable UID. */
+    private static final long serialVersionUID = -3804833976337652097L;
+
+    /** Maneuver date. */
+    final AbsoluteDate t0;
+
+    /** Mass change ratio. */
+    private final double massRatio;
+
+    /** Type of orbit used for internal Jacobians. */
+    final OrbitType type;
+
+    /** Keplerian (or equinoctial) differential change due to maneuver. */
+    final double[] delta;
+
+    /** Mean motion change. */
+    final double deltaN;
+
+    /** Build a maneuver defined in spacecraft frame.
+     * @param state0 state at maneuver date, <em>before</em> the maneuver
+     * is performed
+     * @param dV velocity increment in spacecraft frame
+     * @param isp engine specific impulse (s)
+     * @exception OrekitException if spacecraft frame cannot be transformed
+     */
+    
+    public SmallManeuverAnalyticalModel(final SpacecraftState state0,
+                                        final Vector3D dV, final double isp)
+        throws OrekitException {
+        this(state0, state0.getFrame(),
+             state0.getAttitude().getRotation().applyInverseTo(dV),
+             isp);
+    }
+
+    /** Build a maneuver defined in user-specified frame.
+     * @param state0 state at maneuver date, <em>before</em> the maneuver
+     * is performed
+     * @param frame frame in which velocity increment is defined
+     * @param dV velocity increment in specified frame
+     * @param isp engine specific impulse (s)
+     * @exception OrekitException if velocity increment frame cannot be transformed
+     */
+    
+    public SmallManeuverAnalyticalModel(final SpacecraftState state0, final Frame frame,
+                                        final Vector3D dV, final double isp)
+        throws OrekitException {
+
+        this.t0        = state0.getDate();
+        this.massRatio = FastMath.exp(-dV.getNorm() / (Constants.G0_STANDARD_GRAVITY * isp));
+
+        // use equinoctial orbit type if possible, Keplerian if nearly hyperbolic orbits
+        type = (state0.getE() < 0.9) ? OrbitType.EQUINOCTIAL : OrbitType.KEPLERIAN;
+
+        // compute initial Jacobian
+        final double[][] jacobian = new double[6][6];
+        final Orbit orbit0 = type.convertType(state0.getOrbit());
+        orbit0.getJacobianWrtCartesian(PositionAngle.MEAN, jacobian);
+
+        // compute maneuver effect on Keplerian (or equinoctial) elements
+        final Vector3D inertDV = frame.getTransformTo(state0.getFrame(), t0).transformVector(dV);
+        delta = new double[6];
+        for (int i = 0; i < delta.length; ++i) {
+            delta[i] = jacobian[i][3] * inertDV.getX() +
+                       jacobian[i][4] * inertDV.getY() +
+                       jacobian[i][5] * inertDV.getZ();
+        }
+
+        // compute mean motion change: dM(t1) = dM(t0) + deltaN * (t1 - t0)
+        final double mu = state0.getMu();
+        final double a  = state0.getA();
+        deltaN = -1.5 * delta[0] * FastMath.sqrt(mu / a) / (a * a);
+
+    }
+
+    /** Compute the effect of the maneuver on a spacecraft state.
+     * @param state1 original spacecraft state at t<sub>1</sub>,
+     * without maneuver
+     * @return spacecraft state at t<sub>1</sub>, taking the maneuver
+     * into account if t<sub>1</sub> &gt; t<sub>0</sub>
+     */
+    public SpacecraftState applyManeuver(final SpacecraftState state1) {
+
+        if (state1.getDate().compareTo(t0) <= 0) {
+            // the maneuver has not occurred yet, don't change anything
+            return state1;
+        }
+
+        // convert current orbital state to Keplerian or equinoctial elements
+        final double[] parameters = new double[6];
+        type.mapOrbitToArray(type.convertType(state1.getOrbit()),
+                             PositionAngle.MEAN, parameters);
+
+        // add maneuver effect
+        for (int i = 0; i < parameters.length; ++i) {
+            parameters[i] += delta[i];
+        }
+        parameters[5] += deltaN * state1.getDate().durationFrom(t0);
+
+        // build updated orbit as Keplerian or equinoctial elements
+        final Orbit updated = type.mapArrayToOrbit(parameters, PositionAngle.MEAN,
+                                                   state1.getDate(), state1.getMu(),
+                                                   state1.getFrame());
+
+        // return a new state with the same orbit type as provided
+        return new SpacecraftState(state1.getOrbit().getType().convertType(updated),
+                                   state1.getAttitude(), massRatio * state1.getMass());
+
+    }
+
+}

src/main/java/org/orekit/overview.html

@@ -62,9 +62,17 @@ discrete events. Here is a short list of the features offered by the library:</p
       <li>transparent conversion between all parameters</li>
       <li>automatic binding with frames</li>
       <li>attitude state and derivative</li>
+      <li>Jacobians</li>
       <li>mass management</li>
     </ul>
   </li>
+  <li>Maneuvers
+    <ul>
+      <li>analytical models for small maneuvers without propagation</li>
+      <li>impulse maneuvers for any propagator type</li>
+      <li>continuous maneuvers for numerical propagator type</li>
+    </ul>
+  </li>
   <li>Propagation
     <ul>
       <li>analytical propagation models (Kepler, Eckstein-Heschler, SDP4/SGP4 with 2006 corrections)</li>
@@ -78,7 +86,7 @@ discrete events. Here is a short list of the features offered by the library:</p
               <li>atmospheric drag (DTM2000, Jacchia-Bowman 2006 and simple exponential models)</li>
               <li>third body attraction (with data for Sun, Moon and all solar systems planets)</li>
               <li>radiation pressure with eclipses</li>
-              <li>multiple maneuvers
+              <li>multiple maneuvers</li>
             </ul>
           </li>
           <li>state of the art ODE integrators (adaptive stepsize with error control,

src/site/apt/index.apt

@@ -68,8 +68,18 @@ Overview
 
     * attitude state and derivative
 
+    * Jacobians
+
     * mass management
 
+  ** Maneuvers
+
+    * analytical models for small maneuvers without propagation
+
+    * impulse maneuvers for any propagator type
+
+    *continuous maneuvers for numerical propagator type
+
   ** Propagation
 
     * analytical propagation models

src/site/xdoc/changes.xml

@@ -22,6 +22,9 @@
   <body>
     <release version="6.0" date="TBD"
              description="TBD.">
+      <action dev="luc" type="add">
+        Added an analytical model for the effect at date t1 of a small maneuver performed at date t0.
+      </action>
       <action dev="luc" type="fix">
         Fixed a missing reinitialization of start date when state was reset in numerical propagator.
       </action>

src/test/java/org/orekit/forces/maneuvers/SmallManeuverAnalyticalModelTest.java

+/* Copyright 2002-2011 CS Communication & Systèmes
+ * Licensed to CS Communication & Systèmes (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 org.orekit.forces.maneuvers;
+
+
+import org.apache.commons.math.geometry.euclidean.threed.Vector3D;
+import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator;
+import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;
+import org.apache.commons.math.util.FastMath;
+import org.junit.Assert;
+import org.junit.Before;
+import org.junit.Test;
+import org.orekit.Utils;
+import org.orekit.attitudes.AttitudeProvider;
+import org.orekit.attitudes.LofOffset;
+import org.orekit.errors.OrekitException;
+import org.orekit.frames.FramesFactory;
+import org.orekit.frames.LOFType;
+import org.orekit.orbits.CircularOrbit;
+import org.orekit.orbits.KeplerianOrbit;
+import org.orekit.orbits.Orbit;
+import org.orekit.orbits.PositionAngle;
+import org.orekit.propagation.BoundedPropagator;
+import org.orekit.propagation.SpacecraftState;
+import org.orekit.propagation.numerical.NumericalPropagator;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.time.DateComponents;
+import org.orekit.time.TimeComponents;
+import org.orekit.time.TimeScalesFactory;
+import org.orekit.utils.Constants;
+import org.orekit.utils.PVCoordinates;
+
+public class SmallManeuverAnalyticalModelTest {
+
+    @Test
+    public void testLowEarthOrbit() throws OrekitException {
+
+        Orbit leo = new CircularOrbit(7200000.0, -1.0e-5, 2.0e-4,
+                                      FastMath.toRadians(98.0),
+                                      FastMath.toRadians(123.456),
+                                      0.0, PositionAngle.MEAN,
+                                      FramesFactory.getEME2000(),
+                                      new AbsoluteDate(new DateComponents(2004, 01, 01),
+                                                       new TimeComponents(23, 30, 00.000),
+                                                       TimeScalesFactory.getUTC()),
+                                      Constants.EIGEN5C_EARTH_MU);
+        double mass     = 5600.0;
+        AbsoluteDate t0 = leo.getDate().shiftedBy(1000.0);
+        Vector3D dV     = new Vector3D(-0.01, 0.02, 0.03);
+        double f        = 20.0;
+        double isp      = 315.0;
+        BoundedPropagator withoutManeuver = getEphemeris(leo, mass, t0, Vector3D.ZERO, f, isp);
+        BoundedPropagator withManeuver    = getEphemeris(leo, mass, t0, dV, f, isp);
+        SmallManeuverAnalyticalModel model =
+                new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0), dV, isp);
+
+        for (AbsoluteDate t = withoutManeuver.getMinDate();
+             t.compareTo(withoutManeuver.getMaxDate()) < 0;
+             t = t.shiftedBy(60.0)) {
+            PVCoordinates pvWithout = withoutManeuver.getPVCoordinates(t, leo.getFrame());
+            PVCoordinates pvWith    = withManeuver.getPVCoordinates(t, leo.getFrame());
+            PVCoordinates pvModel   = model.applyManeuver(withoutManeuver.propagate(t)).getPVCoordinates(leo.getFrame());
+            double nominalDeltaP    = new PVCoordinates(pvWith, pvWithout).getPosition().getNorm();
+            double modelError       = new PVCoordinates(pvWith, pvModel).getPosition().getNorm();
+            if (t.compareTo(t0) < 0) {
+                // before maneuver, all positions should be equal
+                Assert.assertEquals(0, nominalDeltaP, 1.0e-10);
+                Assert.assertEquals(0, modelError,    1.0e-10);
+            } else {
+                // after maneuver, model error should be less than 0.8m,
+                // despite nominal deltaP exceeds 1 kilometer after less than 3 orbits
+                if (t.durationFrom(t0) > 0.1 * leo.getKeplerianPeriod()) {
+                    Assert.assertTrue(modelError < 0.009 * nominalDeltaP);
+                }
+                Assert.assertTrue(modelError < 0.8);
+            }
+        }
+
+    }
+
+    @Test
+    public void testEccentricOrbit() throws OrekitException {
+
+        Orbit heo = new KeplerianOrbit(90000000.0, 0.92, FastMath.toRadians(98.0),
+                                       FastMath.toRadians(12.3456),
+                                       FastMath.toRadians(123.456),
+                                       FastMath.toRadians(1.23456), PositionAngle.MEAN,
+                                       FramesFactory.getEME2000(),
+                                       new AbsoluteDate(new DateComponents(2004, 01, 01),
+                                                        new TimeComponents(23, 30, 00.000),
+                                                        TimeScalesFactory.getUTC()),
+                                                        Constants.EIGEN5C_EARTH_MU);
+        double mass     = 5600.0;
+        AbsoluteDate t0 = heo.getDate().shiftedBy(1000.0);
+        Vector3D dV     = new Vector3D(-0.01, 0.02, 0.03);
+        double f        = 20.0;
+        double isp      = 315.0;
+        BoundedPropagator withoutManeuver = getEphemeris(heo, mass, t0, Vector3D.ZERO, f, isp);
+        BoundedPropagator withManeuver    = getEphemeris(heo, mass, t0, dV, f, isp);
+        SmallManeuverAnalyticalModel model =
+                new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0), dV, isp);
+
+        for (AbsoluteDate t = withoutManeuver.getMinDate();
+             t.compareTo(withoutManeuver.getMaxDate()) < 0;
+             t = t.shiftedBy(600.0)) {
+            PVCoordinates pvWithout = withoutManeuver.getPVCoordinates(t, heo.getFrame());
+            PVCoordinates pvWith    = withManeuver.getPVCoordinates(t, heo.getFrame());
+            PVCoordinates pvModel   = model.applyManeuver(withoutManeuver.propagate(t)).getPVCoordinates(heo.getFrame());
+            double nominalDeltaP    = new PVCoordinates(pvWith, pvWithout).getPosition().getNorm();
+            double modelError       = new PVCoordinates(pvWith, pvModel).getPosition().getNorm();
+            if (t.compareTo(t0) < 0) {
+                // before maneuver, all positions should be equal
+                Assert.assertEquals(0, nominalDeltaP, 1.0e-10);
+                Assert.assertEquals(0, modelError,    1.0e-10);
+            } else {
+                // after maneuver, model error should be less than 1700m,
+                // despite nominal deltaP exceeds 300 kilometers at perigee, after 3 orbits
+                if (t.durationFrom(t0) > 0.01 * heo.getKeplerianPeriod()) {
+                    Assert.assertTrue(modelError < 0.005 * nominalDeltaP);
+                }
+                Assert.assertTrue(modelError < 1700);
+            }
+        }
+
+    }
+
+    private BoundedPropagator getEphemeris(final Orbit orbit, final double mass,
+                                           final AbsoluteDate t0, final Vector3D dV,
+                                           final double f, final double isp)
+        throws OrekitException {
+
+        final AttitudeProvider law = new LofOffset(orbit.getFrame(), LOFType.LVLH);
+        final SpacecraftState initialState =
+            new SpacecraftState(orbit, law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
+
+
+        // set up numerical propagator
+        final double dP = 1.0;
+        double[][] tolerances = NumericalPropagator.tolerances(dP, orbit, orbit.getType());
+        AdaptiveStepsizeIntegrator integrator =
+            new DormandPrince853Integrator(0.001, 1000, tolerances[0], tolerances[1]);
+        integrator.setInitialStepSize(orbit.getKeplerianPeriod() / 100.0);
+        final NumericalPropagator propagator = new NumericalPropagator(integrator);
+        propagator.setInitialState(initialState);
+        propagator.setAttitudeProvider(law);
+
+        if (dV.getNorm() > 1.0e-6) {
+            // set up maneuver
+            final double vExhaust = Constants.G0_STANDARD_GRAVITY * isp;
+            final double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust);
+            final ConstantThrustManeuver maneuver =
+                    new ConstantThrustManeuver(t0, dt , f, isp, dV.normalize());
+            propagator.addForceModel(maneuver);
+        }
+
+        propagator.setEphemerisMode();
+        propagator.propagate(t0.shiftedBy(5 * orbit.getKeplerianPeriod()));
+        return propagator.getGeneratedEphemeris();
+
+    }
+
+    @Before
+    public void setUp() {
+        Utils.setDataRoot("regular-data");
+    }
+
+}