SmallManeuverAnalyticalModel.java

  1. /* Copyright 2002-2025 CS GROUP
  2.  * Licensed to CS GROUP (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.forces.maneuvers;

  19. import java.util.Arrays;

  21. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  22. import org.hipparchus.util.FastMath;
  23. import org.orekit.frames.Frame;
  24. import org.orekit.orbits.Orbit;
  25. import org.orekit.orbits.OrbitType;
  26. import org.orekit.orbits.PositionAngleType;
  27. import org.orekit.propagation.SpacecraftState;
  28. import org.orekit.propagation.analytical.AdapterPropagator;
  29. import org.orekit.time.AbsoluteDate;
  30. import org.orekit.utils.Constants;

  32. /** Analytical model for small maneuvers.
  33.  * <p>The aim of this model is to compute quickly the effect at date t₁
  34.  * of a small maneuver performed at an earlier date t₀. Both the
  35.  * direct effect of the maneuver and the Jacobian of this effect with respect to
  36.  * maneuver parameters are available.
  37.  * </p>
  38.  * <p>These effect are computed analytically using two Jacobian matrices:
  39.  * <ol>
  40.  *   <li>J₀: Jacobian of Keplerian or equinoctial elements with respect
  41.  *   to Cartesian parameters at date t₀ allows to compute
  42.  *   maneuver effect as a change in orbital elements at maneuver date t₀,</li>
  43.  *   <li>J<sub>1/0</sub>: Jacobian of Keplerian or equinoctial elements
  44.  *   at date t₁ with respect to Keplerian or equinoctial elements
  45.  *   at date t₀ allows to propagate the change in orbital elements
  46.  *   to final date t₁.</li>
  47.  * </ol>
  48.  *
  49.  * <p>
  50.  * The second Jacobian, J<sub>1/0</sub>, is computed using a simple Keplerian
  51.  * model, i.e. it is the identity except for the mean motion row which also includes
  52.  * an off-diagonal element due to semi-major axis change.
  53.  * </p>
  54.  * <p>
  55.  * The orbital elements change at date t₁ can be added to orbital elements
  56.  * extracted from state, and the final elements taking account the changes are then
  57.  * converted back to appropriate type, which may be different from Keplerian or
  58.  * equinoctial elements.
  59.  * </p>
  60.  * <p>
  61.  * Note that this model takes <em>only</em> Keplerian effects into account. This means
  62.  * that using only this class to compute an inclination maneuver in Low Earth Orbit will
  63.  * <em>not</em> change ascending node drift rate despite inclination has changed (the
  64.  * same would be true for a semi-major axis change of course). In order to take this
  65.  * drift into account, an instance of {@link
  66.  * org.orekit.propagation.analytical.J2DifferentialEffect J2DifferentialEffect}
  67.  * must be used together with an instance of this class.
  68.  * </p>
  69.  * @author Luc Maisonobe
  70.  */
  71. public class SmallManeuverAnalyticalModel implements AdapterPropagator.DifferentialEffect {

  73.     /** State at maneuver date (before maneuver occurrence). */
  74.     private final SpacecraftState state0;

  76.     /** Inertial velocity increment. */
  77.     private final Vector3D inertialDV;

  79.     /** Mass change ratio. */
  80.     private final double massRatio;

  82.     /** Type of orbit used for internal Jacobians. */
  83.     private final OrbitType type;

  85.     /** Initial Keplerian (or equinoctial) Jacobian with respect to maneuver. */
  86.     private final double[][] j0;

  88.     /** Time derivative of the initial Keplerian (or equinoctial) Jacobian with respect to maneuver. */
  89.     private double[][] j0Dot;

  91.     /** Mean anomaly change factor. */
  92.     private final double ksi;

  94.     /** Build a maneuver defined in spacecraft frame with default orbit type.
  95.      * @param state0 state at maneuver date, <em>before</em> the maneuver
  96.      * is performed
  97.      * @param dV velocity increment in spacecraft frame
  98.      * @param isp engine specific impulse (s)
  99.      */
  100.     public SmallManeuverAnalyticalModel(final SpacecraftState state0,
  101.                                         final Vector3D dV, final double isp) {
  102.         this(state0, state0.getFrame(),
  103.              state0.getAttitude().getRotation().applyInverseTo(dV),
  104.              isp);
  105.     }

  107.     /** Build a maneuver defined in spacecraft frame.
  108.      * @param state0 state at maneuver date, <em>before</em> the maneuver
  109.      * is performed
  110.      * @param orbitType orbit type to be used later on in Jacobian conversions
  111.      * @param dV velocity increment in spacecraft frame
  112.      * @param isp engine specific impulse (s)
  113.      * @since 12.1 orbit type added as input
  114.      */
  115.     public SmallManeuverAnalyticalModel(final SpacecraftState state0, final OrbitType orbitType,
  116.                                         final Vector3D dV, final double isp) {
  117.         this(state0, orbitType, state0.getFrame(),
  118.              state0.getAttitude().getRotation().applyInverseTo(dV),
  119.              isp);
  120.     }

  122.     /** Build a maneuver defined in user-specified frame.
  123.      * @param state0 state at maneuver date, <em>before</em> the maneuver
  124.      * is performed
  125.      * @param frame frame in which velocity increment is defined
  126.      * @param dV velocity increment in specified frame
  127.      * @param isp engine specific impulse (s)
  128.      */
  129.     public SmallManeuverAnalyticalModel(final SpacecraftState state0, final Frame frame,
  130.                                         final Vector3D dV, final double isp) {
  131.         // No orbit type specified, use equinoctial orbit type if possible, Keplerian if nearly hyperbolic orbits
  132.         this(state0, (state0.getOrbit().getE() < 0.9) ? OrbitType.EQUINOCTIAL : OrbitType.KEPLERIAN, frame, dV, isp);
  133.     }

  135.     /** Build a maneuver defined in user-specified frame.
  136.      * @param state0 state at maneuver date, <em>before</em> the maneuver
  137.      * is performed
  138.      * @param orbitType orbit type to be used later on in Jacobian conversions
  139.      * @param frame frame in which velocity increment is defined
  140.      * @param dV velocity increment in specified frame
  141.      * @param isp engine specific impulse (s)
  142.      * @since 12.1 orbit type added as input
  143.      */
  144.     public SmallManeuverAnalyticalModel(final SpacecraftState state0, final OrbitType orbitType,
  145.                                         final Frame frame, final Vector3D dV, final double isp) {

  147.         this.state0    = state0;
  148.         this.massRatio = FastMath.exp(-dV.getNorm() / (Constants.G0_STANDARD_GRAVITY * isp));
  149.         this.type = orbitType;

  151.         // compute initial Jacobian
  152.         final double[][] fullJacobian = new double[6][6];
  153.         j0 = new double[6][3];
  154.         final Orbit orbit0 = orbitType.convertType(state0.getOrbit());
  155.         orbit0.getJacobianWrtCartesian(PositionAngleType.MEAN, fullJacobian);
  156.         for (int i = 0; i < j0.length; ++i) {
  157.             System.arraycopy(fullJacobian[i], 3, j0[i], 0, 3);
  158.         }

  160.         // use lazy evaluation for j0Dot, as it is used only when Jacobians are evaluated
  161.         j0Dot = null;

  163.         // compute maneuver effect on Keplerian (or equinoctial) elements
  164.         inertialDV = frame.getStaticTransformTo(state0.getFrame(), state0.getDate())
  165.                         .transformVector(dV);

  167.         // compute mean anomaly change: dM(t1) = dM(t0) + ksi * da * (t1 - t0)
  168.         final Orbit orbit = state0.getOrbit();
  169.         final double mu = orbit.getMu();
  170.         final double a  = orbit.getA();
  171.         ksi = -1.5 * FastMath.sqrt(mu / a) / (a * a);

  173.     }

  175.     /** Get the date of the maneuver.
  176.      * @return date of the maneuver
  177.      */
  178.     public AbsoluteDate getDate() {
  179.         return state0.getDate();
  180.     }

  182.     /** Get the inertial velocity increment of the maneuver.
  183.      * @return velocity increment in a state-dependent inertial frame
  184.      * @see #getInertialFrame()
  185.      */
  186.     public Vector3D getInertialDV() {
  187.         return inertialDV;
  188.     }

  190.     /** Get the inertial frame in which the velocity increment is defined.
  191.      * @return inertial frame in which the velocity increment is defined
  192.      * @see #getInertialDV()
  193.      */
  194.     public Frame getInertialFrame() {
  195.         return state0.getFrame();
  196.     }

  198.     /** Compute the effect of the maneuver on an orbit.
  199.      * @param orbit1 original orbit at t₁, without maneuver
  200.      * @return orbit at t₁, taking the maneuver
  201.      * into account if t₁ &gt; t₀
  202.      * @see #apply(SpacecraftState)
  203.      * @see #getJacobian(Orbit, PositionAngleType, double[][])
  204.      */
  205.     public Orbit apply(final Orbit orbit1) {

  207.         if (orbit1.getDate().compareTo(state0.getDate()) <= 0) {
  208.             // the maneuver has not occurred yet, don't change anything
  209.             return orbit1;
  210.         }

  212.         return orbit1.getType().convertType(updateOrbit(orbit1));

  214.     }

  216.     /** Compute the effect of the maneuver on a spacecraft state.
  217.      * @param state1 original spacecraft state at t₁,
  218.      * without maneuver
  219.      * @return spacecraft state at t₁, taking the maneuver
  220.      * into account if t₁ &gt; t₀
  221.      * @see #apply(Orbit)
  222.      * @see #getJacobian(Orbit, PositionAngleType, double[][])
  223.      */
  224.     public SpacecraftState apply(final SpacecraftState state1) {

  226.         if (state1.getDate().compareTo(state0.getDate()) <= 0) {
  227.             // the maneuver has not occurred yet, don't change anything
  228.             return state1;
  229.         }

  231.         return new SpacecraftState(state1.getOrbit().getType().convertType(updateOrbit(state1.getOrbit())),
  232.                                    state1.getAttitude()).withMass(updateMass(state1.getMass()));

  234.     }

  236.     /** Compute the effect of the maneuver on an orbit.
  237.      * @param orbit1 original orbit at t₁, without maneuver
  238.      * @return orbit at t₁, always taking the maneuver into account, always in the internal type
  239.      */
  240.     private Orbit updateOrbit(final Orbit orbit1) {

  242.         // compute maneuver effect
  243.         final double dt = orbit1.getDate().durationFrom(state0.getDate());
  244.         final double x  = inertialDV.getX();
  245.         final double y  = inertialDV.getY();
  246.         final double z  = inertialDV.getZ();
  247.         final double[] delta = new double[6];
  248.         for (int i = 0; i < delta.length; ++i) {
  249.             delta[i] = j0[i][0] * x + j0[i][1] * y + j0[i][2] * z;
  250.         }
  251.         delta[5] += ksi * delta[0] * dt;

  253.         // convert current orbital state to Keplerian or equinoctial elements
  254.         final double[] parameters    = new double[6];
  255.         type.mapOrbitToArray(type.convertType(orbit1), PositionAngleType.MEAN, parameters, null);
  256.         for (int i = 0; i < delta.length; ++i) {
  257.             parameters[i] += delta[i];
  258.         }

  260.         // build updated orbit as Keplerian or equinoctial elements
  261.         return type.mapArrayToOrbit(parameters, null, PositionAngleType.MEAN,
  262.                                     orbit1.getDate(), orbit1.getMu(), orbit1.getFrame());

  264.     }

  266.     /** Compute the Jacobian of the orbit with respect to maneuver parameters.
  267.      * <p>
  268.      * The Jacobian matrix is a 6x4 matrix. Element jacobian[i][j] corresponds to
  269.      * the partial derivative of orbital parameter i with respect to maneuver
  270.      * parameter j. The rows order is the same order as used in {@link
  271.      * Orbit#getJacobianWrtCartesian(PositionAngleType, double[][]) Orbit.getJacobianWrtCartesian}
  272.      * method. Columns (0, 1, 2) correspond to the velocity increment coordinates
  273.      * (ΔV<sub>x</sub>, ΔV<sub>y</sub>, ΔV<sub>z</sub>) in the
  274.      * inertial frame returned by {@link #getInertialFrame()}, and column 3
  275.      * corresponds to the maneuver date t₀.
  276.      * </p>
  277.      * @param orbit1 original orbit at t₁, without maneuver
  278.      * @param positionAngleType type of the position angle to use
  279.      * @param jacobian placeholder 6x4 (or larger) matrix to be filled with the Jacobian, if matrix
  280.      * is larger than 6x4, only the 6x4 upper left corner will be modified
  281.      * @see #apply(Orbit)
  282.      */
  283.     public void getJacobian(final Orbit orbit1, final PositionAngleType positionAngleType,
  284.                             final double[][] jacobian) {

  286.         final double dt = orbit1.getDate().durationFrom(state0.getDate());
  287.         if (dt < 0) {
  288.             // the maneuver has not occurred yet, Jacobian is null
  289.             for (int i = 0; i < 6; ++i) {
  290.                 Arrays.fill(jacobian[i], 0, 4, 0.0);
  291.             }
  292.             return;
  293.         }

  295.         // derivatives of Keplerian/equinoctial elements with respect to velocity increment
  296.         final double x  = inertialDV.getX();
  297.         final double y  = inertialDV.getY();
  298.         final double z  = inertialDV.getZ();
  299.         for (int i = 0; i < 6; ++i) {
  300.             System.arraycopy(j0[i], 0, jacobian[i], 0, 3);
  301.         }
  302.         for (int j = 0; j < 3; ++j) {
  303.             jacobian[5][j] += ksi * dt * j0[0][j];
  304.         }

  306.         // derivatives of Keplerian/equinoctial elements with respect to date
  307.         evaluateJ0Dot();
  308.         for (int i = 0; i < 6; ++i) {
  309.             jacobian[i][3] = j0Dot[i][0] * x + j0Dot[i][1] * y + j0Dot[i][2] * z;
  310.         }
  311.         final double da = j0[0][0] * x + j0[0][1] * y + j0[0][2] * z;
  312.         jacobian[5][3] += ksi * (jacobian[0][3] * dt - da);

  314.         if (orbit1.getType() != type || positionAngleType != PositionAngleType.MEAN) {

  316.             // convert to derivatives of Cartesian parameters
  317.             final double[][] j2         = new double[6][6];
  318.             final double[][] pvJacobian = new double[6][4];
  319.             final Orbit updated         = updateOrbit(orbit1);
  320.             updated.getJacobianWrtParameters(PositionAngleType.MEAN, j2);
  321.             for (int i = 0; i < 6; ++i) {
  322.                 for (int j = 0; j < 4; ++j) {
  323.                     pvJacobian[i][j] = j2[i][0] * jacobian[0][j] + j2[i][1] * jacobian[1][j] +
  324.                                     j2[i][2] * jacobian[2][j] + j2[i][3] * jacobian[3][j] +
  325.                                     j2[i][4] * jacobian[4][j] + j2[i][5] * jacobian[5][j];
  326.                 }
  327.             }

  329.             // convert to derivatives of specified parameters
  330.             final double[][] j3 = new double[6][6];
  331.             orbit1.getType().convertType(updated).getJacobianWrtCartesian(positionAngleType, j3);
  332.             for (int j = 0; j < 4; ++j) {
  333.                 for (int i = 0; i < 6; ++i) {
  334.                     jacobian[i][j] = j3[i][0] * pvJacobian[0][j] + j3[i][1] * pvJacobian[1][j] +
  335.                                     j3[i][2] * pvJacobian[2][j] + j3[i][3] * pvJacobian[3][j] +
  336.                                     j3[i][4] * pvJacobian[4][j] + j3[i][5] * pvJacobian[5][j];
  337.                 }
  338.             }

  340.         }

  342.     }

  344.     /** Lazy evaluation of the initial Jacobian time derivative.
  345.      */
  346.     private void evaluateJ0Dot() {

  348.         if (j0Dot == null) {

  350.             j0Dot = new double[6][3];
  351.             final double dt = 1.0e-5 / state0.getOrbit().getKeplerianMeanMotion();
  352.             final Orbit orbit = type.convertType(state0.getOrbit());

  354.             // compute shifted Jacobians
  355.             final double[][] j0m1 = new double[6][6];
  356.             orbit.shiftedBy(-1 * dt).getJacobianWrtCartesian(PositionAngleType.MEAN, j0m1);
  357.             final double[][] j0p1 = new double[6][6];
  358.             orbit.shiftedBy(+1 * dt).getJacobianWrtCartesian(PositionAngleType.MEAN, j0p1);

  360.             // evaluate derivative by finite differences
  361.             for (int i = 0; i < j0Dot.length; ++i) {
  362.                 final double[] m1Row    = j0m1[i];
  363.                 final double[] p1Row    = j0p1[i];
  364.                 final double[] j0DotRow = j0Dot[i];
  365.                 for (int j = 0; j < 3; ++j) {
  366.                     j0DotRow[j] = (p1Row[j + 3] - m1Row[j + 3]) / (2 * dt);
  367.                 }
  368.             }

  370.         }

  372.     }

  374.     /** Update a spacecraft mass due to maneuver.
  375.      * @param mass masse before maneuver
  376.      * @return mass after maneuver
  377.      */
  378.     public double updateMass(final double mass) {
  379.         return massRatio * mass;
  380.     }

  382. }
Created with JaCoCo 0.8.12.202403310830