TriggerDate.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.jacobians;

  19. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  20. import org.hipparchus.linear.MatrixUtils;
  21. import org.hipparchus.linear.QRDecomposition;
  22. import org.hipparchus.linear.RealMatrix;
  23. import org.hipparchus.linear.RealVector;
  24. import org.orekit.forces.maneuvers.Maneuver;
  25. import org.orekit.forces.maneuvers.trigger.ManeuverTriggersResetter;
  26. import org.orekit.propagation.AdditionalDataProvider;
  27. import org.orekit.propagation.SpacecraftState;
  28. import org.orekit.propagation.integration.AdditionalDerivativesProvider;
  29. import org.orekit.time.AbsoluteDate;
  30. import org.orekit.utils.TimeSpanMap;

  32. /** Generator for one column of a Jacobian matrix for special case of trigger dates.
  33.  * <p>
  34.  * Typical use cases for this are estimation of maneuver start and stop date during
  35.  * either orbit determination or maneuver optimization.
  36.  * </p>
  37.  * <p>
  38.  * Let \((t_0, y_0)\) be the state at propagation start, \((t_1, y_1)\) be the state at
  39.  * maneuver trigger time, \((t_t, y_t)\) be the state at any arbitrary time \(t\) during
  40.  * propagation, and \(f_m(t, y)\) be the contribution of the maneuver to the global
  41.  * ODE \(\frac{dy}{dt} = f(t, y)\). We are interested in the Jacobian column
  42.  * \(\frac{\partial y_t}{\partial t_1}\).
  43.  * </p>
  44.  * <p>
  45.  * There are two parts in this Jacobian: the primary part corresponds to the full contribution
  46.  * of the acceleration due to the maneuver as it is delayed by a small amount \(dt_1\), whereas
  47.  * the secondary part corresponds to change of acceleration after maneuver start as the mass
  48.  * depletion is delayed and therefore the spacecraft mass is different from the mass for nominal
  49.  * start time.
  50.  * </p>
  51.  * <p>
  52.  * The primary part is computed as follows. After trigger time \(t_1\) (according to propagation direction),
  53.  * \[\frac{\partial y_t}{\partial t_1} = \pm \frac{\partial y_t}{\partial y_1} f_m(t_1, y_1)\]
  54.  * where the sign depends on \(t_1\) being a start or stop trigger and propagation being forward
  55.  * or backward.
  56.  * </p>
  57.  * <p>
  58.  * We don't have \(\frac{\partial y_t}{\partial y_1}\) available if \(t_1 \neq t_0\), but we
  59.  * have \(\frac{\partial y_t}{\partial y_0}\) at any time since it can be computed by integrating
  60.  * variational equations for numerical propagation or by other closed form expressions for analytical
  61.  * propagators. We use the classical composition rule to recover the state transition matrix with
  62.  * respect to intermediate time \(t_1\):
  63.  * \[\frac{\partial y_t}{\partial y_0} = \frac{\partial y_t}{\partial y_1} \frac{\partial y_1}{\partial y_0}\]
  64.  * We deduce
  65.  * \[\frac{\partial y_t}{\partial y_1} = \frac{\partial y_t}{\partial y_0} \left(\frac{\partial y_1}{\partial y_0}\right)^{-1}\]
  66.  * </p>
  67.  * <p>
  68.  * The contribution of the primary part to the Jacobian column can therefore be computed using the following
  69.  * closed-form expression:
  70.  * \[\frac{\partial y_t}{\partial t_1}
  71.  * = \pm \frac{\partial y_t}{\partial y_0} \left(\frac{\partial y_1}{\partial y_0}\right)^{-1} f_m(t_1, y_1)
  72.  * = \frac{\partial y_t}{\partial y_0} c_1\]
  73.  * where \(c_1\) is the signed contribution of maneuver at \(t_1\) and is computed at trigger time
  74.  * by solving \(\frac{\partial y_1}{\partial y_0} c_1 = \pm f_m(t_1, y_1)\).
  75.  * </p>
  76.  * <p>
  77.  * As the primary part of the column is generated using a closed-form expression, this generator
  78.  * implements the {@link AdditionalDataProvider} interface and stores the column directly
  79.  * in the primary state during propagation.
  80.  * </p>
  81.  * <p>
  82.  * As the closed-form expression requires picking \(c_1\) at trigger time \(t_1\), it works only
  83.  * if propagation starts outside of the maneuver and passes over \(t_1\) during integration.
  84.  * </p>
  85.  * <p>
  86.  * The secondary part is computed as follows. We have acceleration \(\vec{\Gamma} = \frac{\vec{F}}{m}\) and
  87.  * \(m = m_0 - q (t - t_s)\), where \(m\) is current mass, \(m_0\) is initial mass and \(t_s\) is
  88.  * maneuver trigger time. A delay \(dt_s\) on trigger time induces delaying mass depletion.
  89.  * We get:
  90.  * \[d\vec{\Gamma} = \frac{-\vec{F}}{m^2} dm = \frac{-\vec{F}}{m^2} q dt_s = -\vec{\Gamma}\frac{q}{m} dt_s\]
  91.  * From this total differential, we extract the partial derivative of the acceleration
  92.  * \[\frac{\partial\vec{\Gamma}}{\partial t_s} = -\vec{\Gamma}\frac{q}{m}\]
  93.  * </p>
  94.  * <p>
  95.  * The contribution of the secondary part to the Jacobian column can therefore be computed by integrating
  96.  * the partial derivative of the acceleration, to get the partial derivative of the position.
  97.  * </p>
  98.  * <p>
  99.  * As the secondary part of the column is generated using a differential equation, a separate
  100.  * underlying generator implementing the {@link AdditionalDerivativesProvider} interface is set up to
  101.  * perform the integration during propagation.
  102.  * </p>
  103.  * <p>
  104.  * This generator takes care to sum up the primary and secondary parts so the full column of the Jacobian
  105.  * is computed.
  106.  * </p>
  107.  * <p>
  108.  * The implementation takes care to <em>not</em> resetting \(c_1\) at propagation start.
  109.  * This allows to get proper Jacobian if we interrupt propagation in the middle of a maneuver
  110.  * and restart propagation where it left.
  111.  * </p>
  112.  * @author Luc Maisonobe
  113.  * @since 11.1
  114.  * @see MedianDate
  115.  * @see Duration
  116.  */
  117. public class TriggerDate
  118.     implements ManeuverTriggersResetter, AdditionalDataProvider<double[]> {

  120.     /** Dimension of the state. */
  121.     private static final int STATE_DIMENSION = 6;

  123.     /** Threshold for decomposing state transition matrix at trigger time. */
  124.     private static final double DECOMPOSITION_THRESHOLD = 1.0e-10;

  126.     /** Name of the state for State Transition Matrix. */
  127.     private final String stmName;

  129.     /** Name of the parameter corresponding to the column. */
  130.     private final String triggerName;

  132.     /** Mass depletion effect. */
  133.     private final MassDepletionDelay massDepletionDelay;

  135.     /** Start/stop management flag. */
  136.     private final boolean manageStart;

  138.     /** Maneuver force model. */
  139.     private final Maneuver maneuver;

  141.     /** Event detector threshold. */
  142.     private final double threshold;

  144.     /** Signed contribution of maneuver at trigger time ±(∂y₁/∂y₀)⁻¹ fₘ(t₁, y₁). */
  145.     private TimeSpanMap<double[]> contribution;

  147.     /** Trigger date. */
  148.     private AbsoluteDate trigger;

  150.     /** Indicator for forward propagation. */
  151.     private boolean forward;

  153.     /** Simple constructor.
  154.      * @param stmName name of State Transition Matrix state
  155.      * @param triggerName name of the parameter corresponding to the trigger date column
  156.      * @param manageStart if true, we compute derivatives with respect to maneuver start
  157.      * @param maneuver maneuver force model
  158.      * @param threshold event detector threshold
  159.      */
  160.     public TriggerDate(final String stmName, final String triggerName, final boolean manageStart,
  161.                        final Maneuver maneuver, final double threshold) {
  162.         this.stmName            = stmName;
  163.         this.triggerName        = triggerName;
  164.         this.massDepletionDelay = new MassDepletionDelay(triggerName, manageStart, maneuver);
  165.         this.manageStart        = manageStart;
  166.         this.maneuver           = maneuver;
  167.         this.threshold          = threshold;
  168.         this.contribution       = null;
  169.         this.trigger            = null;
  170.         this.forward            = true;
  171.     }

  173.     /** {@inheritDoc} */
  174.     @Override
  175.     public String getName() {
  176.         return triggerName;
  177.     }

  179.     /** {@inheritDoc}
  180.      * <p>
  181.      * The column state can be computed only if the State Transition Matrix state is available.
  182.      * </p>
  183.      */
  184.     @Override
  185.     public boolean yields(final SpacecraftState state) {
  186.         return !(state.hasAdditionalData(stmName) && state.hasAdditionalData(massDepletionDelay.getName()));
  187.     }

  189.     /** Get the mass depletion effect processor.
  190.      * @return mass depletion effect processor
  191.      */
  192.     public MassDepletionDelay getMassDepletionDelay() {
  193.         return massDepletionDelay;
  194.     }

  196.     /** {@inheritDoc} */
  197.     @Override
  198.     public void init(final SpacecraftState initialState, final AbsoluteDate target) {

  200.         // note that we reset contribution or triggered ONLY at start or if we change
  201.         // propagation direction
  202.         // this allows to get proper Jacobian if we interrupt propagation
  203.         // in the middle of a maneuver and restart propagation where it left
  204.         final boolean newForward = target.isAfterOrEqualTo(initialState);
  205.         if (contribution == null || (forward ^ newForward)) {
  206.             contribution = new TimeSpanMap<>(null);
  207.             trigger      = null;
  208.         }

  210.         forward = newForward;

  212.     }

  214.     /** {@inheritDoc} */
  215.     @Override
  216.     public double[] getAdditionalData(final SpacecraftState state) {
  217.         // we check contribution rather than triggered because this method
  218.         // is called after maneuverTriggered and before resetState,
  219.         // when preparing the old state to be reset
  220.         final double[] c = contribution == null ? null : contribution.get(state.getDate());
  221.         if (c == null) {
  222.             // no thrust, no effect
  223.             return new double[STATE_DIMENSION];
  224.         } else {

  226.             // primary effect: full maneuver contribution at (delayed) trigger date
  227.             final double[] effect = getStm(state).operate(c);

  229.             // secondary effect: maneuver change throughout thrust as mass depletion is delayed
  230.             final double[] secondary = state.getAdditionalState(massDepletionDelay.getName());

  232.             // sum up both effects
  233.             for (int i = 0; i < effect.length; ++i) {
  234.                 effect[i] += secondary[i];
  235.             }

  237.             return effect;

  239.         }
  240.     }

  242.     /** {@inheritDoc}*/
  243.     @Override
  244.     public void maneuverTriggered(final SpacecraftState state, final boolean start) {
  245.         trigger = (start == manageStart) ? state.getDate() : null;
  246.     }

  248.     /** {@inheritDoc}*/
  249.     @Override
  250.     public SpacecraftState resetState(final SpacecraftState state) {

  252.         if (trigger == null) {
  253.             // this is not the maneuver trigger we expected (start vs. stop)
  254.             return state;
  255.         }

  257.         // get the acceleration near trigger time
  258.         final SpacecraftState stateWhenFiring = state.shiftedBy((manageStart ? 2 : -2) * threshold);
  259.         final Vector3D        acceleration    = maneuver.acceleration(stateWhenFiring, maneuver.getParameters(state.getDate()));

  261.         // initialize derivatives computation
  262.         final double     sign = (forward == manageStart) ? -1 : +1;
  263.         final RealVector rhs  = MatrixUtils.createRealVector(STATE_DIMENSION);
  264.         rhs.setEntry(3, sign * acceleration.getX());
  265.         rhs.setEntry(4, sign * acceleration.getY());
  266.         rhs.setEntry(5, sign * acceleration.getZ());

  268.         // get State Transition Matrix with respect to Cartesian parameters at trigger time
  269.         final RealMatrix dY1dY0 = getStm(state);

  271.         // store contribution factor for derivatives scm = ±(∂y₁/∂y₀)⁻¹ fₘ(t₁, y₁)
  272.         final double[] c = new QRDecomposition(dY1dY0, DECOMPOSITION_THRESHOLD).getSolver().solve(rhs).toArray();
  273.         if (forward) {
  274.             contribution.addValidAfter(c, state.getDate(), false);
  275.         } else {
  276.             contribution.addValidBefore(c, state.getDate(), false);
  277.         }

  279.         // return unchanged state
  280.         return state;

  282.     }

  284.     /** Extract State Transition Matrix with respect to Cartesian parameters.
  285.      * @param state state containing the State Transition Matrix
  286.      * @return State Transition Matrix
  287.      */
  288.     private RealMatrix getStm(final SpacecraftState state) {
  289.         final double[] p = state.getAdditionalState(stmName);
  290.         final RealMatrix dYdY0 = MatrixUtils.createRealMatrix(STATE_DIMENSION, STATE_DIMENSION);
  291.         int index = 0;
  292.         for (int i = 0; i < STATE_DIMENSION; ++i) {
  293.             for (int j = 0; j < STATE_DIMENSION; ++j) {
  294.                 dYdY0.setEntry(i, j, p[index++]);
  295.             }
  296.         }
  297.         return dYdY0;
  298.     }

  300. }
Created with JaCoCo 0.8.12.202403310830