AbstractDragForceModel.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.drag;

  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.analysis.differentiation.DSFactory;
  21. import org.hipparchus.analysis.differentiation.DerivativeStructure;
  22. import org.hipparchus.analysis.differentiation.Gradient;
  23. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  24. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  25. import org.orekit.forces.ForceModel;
  26. import org.orekit.frames.Frame;
  27. import org.orekit.frames.StaticTransform;
  28. import org.orekit.models.earth.atmosphere.Atmosphere;
  29. import org.orekit.propagation.FieldSpacecraftState;
  30. import org.orekit.time.AbsoluteDate;
  31. import org.orekit.time.FieldAbsoluteDate;
  32. import org.orekit.utils.FieldPVCoordinates;

  34. import java.util.Arrays;

  36. /**
  37.  * Base class for drag force models.
  38.  * @see DragForce
  39.  * @author Bryan Cazabonne
  40.  * @since 10.2
  41.  */
  42. public abstract class AbstractDragForceModel implements ForceModel {

  44.     /** Atmospheric model. */
  45.     private final Atmosphere atmosphere;

  47.     /**
  48.      * Flag to use (first-order) finite differences instead of automatic differentiation when computing density derivatives w.r.t. position.
  49.      */
  50.     private final boolean useFiniteDifferencesOnDensityWrtPosition;

  52.     /**
  53.      * Constructor with default value for finite differences flag.
  54.      * @param atmosphere atmospheric model
  55.      */
  56.     protected AbstractDragForceModel(final Atmosphere atmosphere) {
  57.         this(atmosphere, true);
  58.     }

  60.     /**
  61.      * Constructor.
  62.      * @param atmosphere atmospheric model
  63.      * @param useFiniteDifferencesOnDensityWrtPosition flag to use finite differences to compute density derivatives w.r.t.
  64.      *                                                 position (is less accurate but can be faster depending on model)
  65.      * @since 12.1
  66.      */
  67.     protected AbstractDragForceModel(final Atmosphere atmosphere, final boolean useFiniteDifferencesOnDensityWrtPosition) {
  68.         this.atmosphere = atmosphere;
  69.         this.useFiniteDifferencesOnDensityWrtPosition = useFiniteDifferencesOnDensityWrtPosition;
  70.     }

  72.     /** Get the atmospheric model.
  73.      * @return atmosphere model
  74.      * @since 12.1
  75.      */
  76.     public Atmosphere getAtmosphere() {
  77.         return atmosphere;
  78.     }

  80.     /** {@inheritDoc} */
  81.     @Override
  82.     public boolean dependsOnPositionOnly() {
  83.         return false;
  84.     }

  86.     /** Check if a field state corresponds to derivatives with respect to state.
  87.      * @param state state to check
  88.      * @param <T> type of the field elements
  89.      * @return true if state corresponds to derivatives with respect to state
  90.      */
  91.     protected <T extends CalculusFieldElement<T>> boolean isDSStateDerivative(final FieldSpacecraftState<T> state) {
  92.         try {
  93.             final DerivativeStructure dsMass = (DerivativeStructure) state.getMass();
  94.             final int o = dsMass.getOrder();
  95.             final int p = dsMass.getFreeParameters();

  97.             // To be in the desired case:
  98.             // Order must be 1 (first order derivatives only)
  99.             // Number of parameters must be 6 (PV), 7 (PV + drag coefficient) or 8 (PV + drag coefficient + lift ratio)
  100.             if (o != 1 || p != 6 && p != 7 && p != 8) {
  101.                 return false;
  102.             }

  104.             // Check that the first 6 parameters are position and velocity
  105.             @SuppressWarnings("unchecked")
  106.             final FieldPVCoordinates<DerivativeStructure> pv = (FieldPVCoordinates<DerivativeStructure>) state.getPVCoordinates();
  107.             return isVariable(pv.getPosition().getX(), 0) &&
  108.                    isVariable(pv.getPosition().getY(), 1) &&
  109.                    isVariable(pv.getPosition().getZ(), 2) &&
  110.                    isVariable(pv.getVelocity().getX(), 3) &&
  111.                    isVariable(pv.getVelocity().getY(), 4) &&
  112.                    isVariable(pv.getVelocity().getZ(), 5);
  113.         } catch (ClassCastException cce) {
  114.             return false;
  115.         }
  116.     }

  118.     /** Check if a field state corresponds to derivatives with respect to state.
  119.      * @param state state to check
  120.      * @param <T> type of the field elements
  121.      * @return true if state corresponds to derivatives with respect to state
  122.      */
  123.     protected <T extends CalculusFieldElement<T>> boolean isGradientStateDerivative(final FieldSpacecraftState<T> state) {
  124.         try {
  125.             final Gradient gMass = (Gradient) state.getMass();
  126.             final int p = gMass.getFreeParameters();

  128.             // To be in the desired case:
  129.             // Number of parameters must be 6 (PV), 7 (PV + drag coefficient) or 8 (PV + drag coefficient + lift ratio)
  130.             if (p != 6 && p != 7 && p != 8) {
  131.                 return false;
  132.             }

  134.             // Check that the first 6 parameters are position and velocity
  135.             @SuppressWarnings("unchecked")
  136.             final FieldPVCoordinates<Gradient> pv = (FieldPVCoordinates<Gradient>) state.getPVCoordinates();
  137.             return isVariable(pv.getPosition().getX(), 0) &&
  138.                    isVariable(pv.getPosition().getY(), 1) &&
  139.                    isVariable(pv.getPosition().getZ(), 2) &&
  140.                    isVariable(pv.getVelocity().getX(), 3) &&
  141.                    isVariable(pv.getVelocity().getY(), 4) &&
  142.                    isVariable(pv.getVelocity().getZ(), 5);
  143.         } catch (ClassCastException cce) {
  144.             return false;
  145.         }
  146.     }

  148.     /**
  149.      * Evaluate the Field density.
  150.      * @param s spacecraft state
  151.      * @return atmospheric density
  152.      * @param <T> field type
  153.      * @since 12.1
  154.      */
  155.     @SuppressWarnings("unchecked")
  156.     protected <T extends CalculusFieldElement<T>> T getFieldDensity(final FieldSpacecraftState<T> s) {
  157.         final FieldAbsoluteDate<T> date     = s.getDate();
  158.         final Frame                frame    = s.getFrame();
  159.         final FieldVector3D<T>     position = s.getPosition();
  160.         if (isGradientStateDerivative(s)) {
  161.             if (useFiniteDifferencesOnDensityWrtPosition) {
  162.                 return (T) this.getGradientDensityWrtStateUsingFiniteDifferences(date.toAbsoluteDate(), frame, (FieldVector3D<Gradient>) position);
  163.             } else {
  164.                 return (T) this.getGradientDensityWrtState(date.toAbsoluteDate(), frame, (FieldVector3D<Gradient>) position);
  165.             }
  166.         } else if (isDSStateDerivative(s)) {
  167.             if (useFiniteDifferencesOnDensityWrtPosition) {
  168.                 return (T) this.getDSDensityWrtStateUsingFiniteDifferences(date.toAbsoluteDate(), frame, (FieldVector3D<DerivativeStructure>) position);
  169.             } else {
  170.                 return (T) this.getDSDensityWrtState(date.toAbsoluteDate(), frame, (FieldVector3D<DerivativeStructure>) position);
  171.             }
  172.         } else {
  173.             return atmosphere.getDensity(date, position, frame);
  174.         }
  175.     }

  177.     /** Check if a derivative represents a specified variable.
  178.      * @param ds derivative to check
  179.      * @param index index of the variable
  180.      * @return true if the derivative represents a specified variable
  181.      */
  182.     protected boolean isVariable(final DerivativeStructure ds, final int index) {
  183.         final double[] derivatives = ds.getAllDerivatives();
  184.         boolean check = true;
  185.         for (int i = 1; i < derivatives.length; ++i) {
  186.             check &= derivatives[i] == ((index + 1 == i) ? 1.0 : 0.0);
  187.         }
  188.         return check;
  189.     }

  191.     /** Check if a derivative represents a specified variable.
  192.      * @param g derivative to check
  193.      * @param index index of the variable
  194.      * @return true if the derivative represents a specified variable
  195.      */
  196.     protected boolean isVariable(final Gradient g, final int index) {
  197.         final double[] derivatives = g.getGradient();
  198.         boolean check = true;
  199.         for (int i = 0; i < derivatives.length; ++i) {
  200.             check &= derivatives[i] == ((index == i) ? 1.0 : 0.0);
  201.         }
  202.         return check;
  203.     }

  205.     /** Compute density and its derivatives.
  206.      * Using finite differences for the derivatives.
  207.      * And doing the actual computation only for the derivatives with respect to position (others are set to 0.).
  208.      * <p>
  209.      * From a theoretical point of view, this method computes the same values
  210.      * as {@link Atmosphere#getDensity(FieldAbsoluteDate, FieldVector3D, Frame)} in the
  211.      * specific case of {@link DerivativeStructure} with respect to state, so
  212.      * it is less general. However, it can be faster depending the Field implementation.
  213.      * </p>
  214.      * <p>
  215.      * The derivatives should be computed with respect to position. The input
  216.      * parameters already take into account the free parameters (6, 7 or 8 depending
  217.      * on derivation with respect to drag coefficient and lift ratio being considered or not)
  218.      * and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
  219.      * with respect to position. Free parameters at indices 3, 4 and 5 correspond
  220.      * to derivatives with respect to velocity (these derivatives will remain zero
  221.      * as the atmospheric density does not depend on velocity). Free parameter
  222.      * at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
  223.      * and/or lift ratio (one of these or both).
  224.      * This 2 last derivatives will remain zero as atmospheric density does not depend on them.
  225.      * </p>
  226.      * @param date current date
  227.      * @param frame inertial reference frame for state (both orbit and attitude)
  228.      * @param position position of spacecraft in inertial frame
  229.      * @return the density and its derivatives
  230.      */
  231.     protected DerivativeStructure getDSDensityWrtStateUsingFiniteDifferences(final AbsoluteDate date,
  232.                                                                              final Frame frame,
  233.                                                                              final FieldVector3D<DerivativeStructure> position) {

  235.         // Retrieve derivation properties for parameter T
  236.         // It is implied here that T is a DerivativeStructure
  237.         // With order 1 and 6, 7 or 8 free parameters
  238.         // This is all checked before in method isStateDerivatives
  239.         final DSFactory factory = position.getX().getFactory();

  241.         // Build a DerivativeStructure using only derivatives with respect to position
  242.         final DSFactory factory3 = new DSFactory(3, 1);
  243.         final FieldVector3D<DerivativeStructure> position3 =
  244.                         new FieldVector3D<>(factory3.variable(0, position.getX().getReal()),
  245.                                             factory3.variable(1,  position.getY().getReal()),
  246.                                             factory3.variable(2,  position.getZ().getReal()));

  248.         // Get atmosphere properties in atmosphere own frame
  249.         final Frame      atmFrame  = atmosphere.getFrame();
  250.         final StaticTransform toBody = frame.getStaticTransformTo(atmFrame, date);
  251.         final FieldVector3D<DerivativeStructure> posBodyDS = toBody.transformPosition(position3);
  252.         final Vector3D   posBody   = posBodyDS.toVector3D();

  254.         // Estimate density model by finite differences and composition
  255.         // Using a delta of 1m
  256.         final double delta  = 1.0;
  257.         final double x      = posBody.getX();
  258.         final double y      = posBody.getY();
  259.         final double z      = posBody.getZ();
  260.         final double rho0   = atmosphere.getDensity(date, posBody, atmFrame);
  261.         final double dRhodX = (atmosphere.getDensity(date, new Vector3D(x + delta, y,         z),         atmFrame) - rho0) / delta;
  262.         final double dRhodY = (atmosphere.getDensity(date, new Vector3D(x,         y + delta, z),         atmFrame) - rho0) / delta;
  263.         final double dRhodZ = (atmosphere.getDensity(date, new Vector3D(x,         y,         z + delta), atmFrame) - rho0) / delta;
  264.         final double[] dXdQ = posBodyDS.getX().getAllDerivatives();
  265.         final double[] dYdQ = posBodyDS.getY().getAllDerivatives();
  266.         final double[] dZdQ = posBodyDS.getZ().getAllDerivatives();

  268.         // Density with derivatives:
  269.         // - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
  270.         // - Others are set to 0.
  271.         final int p = factory.getCompiler().getFreeParameters();
  272.         final double[] rhoAll = new double[p + 1];
  273.         rhoAll[0] = rho0;
  274.         for (int i = 1; i < 4; ++i) {
  275.             rhoAll[i] = dRhodX * dXdQ[i] + dRhodY * dYdQ[i] + dRhodZ * dZdQ[i];
  276.         }

  278.         return factory.build(rhoAll);
  279.     }

  281.     /** Compute density and its derivatives.
  282.      * And doing the actual computation only for the derivatives with respect to position (others are set to 0.).
  283.      * <p>
  284.      * The derivatives should be computed with respect to position. The input
  285.      * parameters already take into account the free parameters (6, 7 or 8 depending
  286.      * on derivation with respect to drag coefficient and lift ratio being considered or not)
  287.      * and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
  288.      * with respect to position. Free parameters at indices 3, 4 and 5 correspond
  289.      * to derivatives with respect to velocity (these derivatives will remain zero
  290.      * as the atmospheric density does not depend on velocity). Free parameter
  291.      * at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
  292.      * and/or lift ratio (one of these or both).
  293.      * This 2 last derivatives will remain zero as atmospheric density does not depend on them.
  294.      * </p>
  295.      * @param date current date
  296.      * @param frame inertial reference frame for state (both orbit and attitude)
  297.      * @param position position of spacecraft in inertial frame
  298.      * @return the density and its derivatives
  299.      */
  300.     protected DerivativeStructure getDSDensityWrtState(final AbsoluteDate date, final Frame frame,
  301.                                                        final FieldVector3D<DerivativeStructure> position) {

  303.         // Retrieve derivation properties for parameter T
  304.         // It is implied here that T is a DerivativeStructure
  305.         // With order 1 and 6, 7 or 8 free parameters
  306.         // This is all checked before in method isStateDerivatives
  307.         final DSFactory factory = position.getX().getFactory();

  309.         // Build a DerivativeStructure using only derivatives with respect to position
  310.         final DSFactory factory3 = new DSFactory(3, 1);
  311.         final FieldVector3D<DerivativeStructure> position3 =
  312.                 new FieldVector3D<>(factory3.variable(0, position.getX().getReal()),
  313.                         factory3.variable(1,  position.getY().getReal()),
  314.                         factory3.variable(2,  position.getZ().getReal()));

  316.         // Get atmosphere properties in atmosphere own frame
  317.         final Frame      atmFrame  = atmosphere.getFrame();
  318.         final StaticTransform toBody = frame.getStaticTransformTo(atmFrame, date);
  319.         final FieldVector3D<DerivativeStructure> posBodyDS = toBody.transformPosition(position3);
  320.         final FieldAbsoluteDate<DerivativeStructure> fieldDate = new FieldAbsoluteDate<>(position3.getX().getField(), date);
  321.         final DerivativeStructure density = atmosphere.getDensity(fieldDate, posBodyDS, atmFrame);

  323.         // Density with derivatives:
  324.         // - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
  325.         // - Others are set to 0.
  326.         final double[] derivatives = Arrays.copyOf(density.getAllDerivatives(), factory.getCompiler().getSize());
  327.         return factory.build(derivatives);
  328.     }

  330.     /** Compute density and its derivatives.
  331.      * Using finite differences for the derivatives.
  332.      * And doing the actual computation only for the derivatives with respect to position (others are set to 0.).
  333.      * <p>
  334.      * From a theoretical point of view, this method computes the same values
  335.      * as {@link Atmosphere#getDensity(FieldAbsoluteDate, FieldVector3D, Frame)} in the
  336.      * specific case of {@link Gradient} with respect to state, so
  337.      * it is less general. However, it can be faster depending the Field implementation.
  338.      * </p>
  339.      * <p>
  340.      * The derivatives should be computed with respect to position. The input
  341.      * parameters already take into account the free parameters (6, 7 or 8 depending
  342.      * on derivation with respect to drag coefficient and lift ratio being considered or not)
  343.      * and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
  344.      * with respect to position. Free parameters at indices 3, 4 and 5 correspond
  345.      * to derivatives with respect to velocity (these derivatives will remain zero
  346.      * as the atmospheric density does not depend on velocity). Free parameter
  347.      * at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
  348.      * and/or lift ratio (one of these or both).
  349.      * This 2 last derivatives will remain zero as atmospheric density does not depend on them.
  350.      * </p>
  351.      * @param date current date
  352.      * @param frame inertial reference frame for state (both orbit and attitude)
  353.      * @param position position of spacecraft in inertial frame
  354.      * @return the density and its derivatives
  355.      */
  356.     protected Gradient getGradientDensityWrtStateUsingFiniteDifferences(final AbsoluteDate date,
  357.                                                                         final Frame frame,
  358.                                                                         final FieldVector3D<Gradient> position) {

  360.         // Build a Gradient using only derivatives with respect to position
  361.         final FieldVector3D<Gradient> position3 =
  362.                         new FieldVector3D<>(Gradient.variable(3, 0, position.getX().getReal()),
  363.                                             Gradient.variable(3, 1,  position.getY().getReal()),
  364.                                             Gradient.variable(3, 2,  position.getZ().getReal()));

  366.         // Get atmosphere properties in atmosphere own frame
  367.         final Frame      atmFrame  = atmosphere.getFrame();
  368.         final StaticTransform toBody = frame.getStaticTransformTo(atmFrame, date);
  369.         final FieldVector3D<Gradient> posBodyDS = toBody.transformPosition(position3);
  370.         final Vector3D   posBody   = posBodyDS.toVector3D();

  372.         // Estimate density model by finite differences and composition
  373.         // Using a delta of 1m
  374.         final double delta  = 1.0;
  375.         final double x      = posBody.getX();
  376.         final double y      = posBody.getY();
  377.         final double z      = posBody.getZ();
  378.         final double rho0   = atmosphere.getDensity(date, posBody, atmFrame);
  379.         final double dRhodX = (atmosphere.getDensity(date, new Vector3D(x + delta, y,         z),         atmFrame) - rho0) / delta;
  380.         final double dRhodY = (atmosphere.getDensity(date, new Vector3D(x,         y + delta, z),         atmFrame) - rho0) / delta;
  381.         final double dRhodZ = (atmosphere.getDensity(date, new Vector3D(x,         y,         z + delta), atmFrame) - rho0) / delta;
  382.         final double[] dXdQ = posBodyDS.getX().getGradient();
  383.         final double[] dYdQ = posBodyDS.getY().getGradient();
  384.         final double[] dZdQ = posBodyDS.getZ().getGradient();

  386.         // Density with derivatives:
  387.         // - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
  388.         // - Others are set to 0.
  389.         final int p = position.getX().getFreeParameters();
  390.         final double[] rhoAll = new double[p];
  391.         for (int i = 0; i < 3; ++i) {
  392.             rhoAll[i] = dRhodX * dXdQ[i] + dRhodY * dYdQ[i] + dRhodZ * dZdQ[i];
  393.         }

  395.         return new Gradient(rho0, rhoAll);
  396.     }

  398.     /** Compute density and its derivatives.
  399.      * <p>
  400.      * The derivatives should be computed with respect to position. The input
  401.      * parameters already take into account the free parameters (6, 7 or 8 depending
  402.      * on derivation with respect to drag coefficient and lift ratio being considered or not)
  403.      * and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
  404.      * with respect to position. Free parameters at indices 3, 4 and 5 correspond
  405.      * to derivatives with respect to velocity (these derivatives will remain zero
  406.      * as the atmospheric density does not depend on velocity). Free parameter
  407.      * at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
  408.      * and/or lift ratio (one of these or both).
  409.      * This 2 last derivatives will remain zero as atmospheric density does not depend on them.
  410.      * </p>
  411.      * @param date current date
  412.      * @param frame inertial reference frame for state (both orbit and attitude)
  413.      * @param position position of spacecraft in inertial frame
  414.      * @return the density and its derivatives
  415.      * @since 12.1
  416.      */
  417.     protected Gradient getGradientDensityWrtState(final AbsoluteDate date, final Frame frame,
  418.                                                   final FieldVector3D<Gradient> position) {

  420.         // Build a Gradient using only derivatives with respect to position
  421.         final int positionDimension = 3;
  422.         final FieldVector3D<Gradient> position3 =
  423.                 new FieldVector3D<>(Gradient.variable(positionDimension, 0, position.getX().getReal()),
  424.                         Gradient.variable(positionDimension, 1,  position.getY().getReal()),
  425.                         Gradient.variable(positionDimension, 2,  position.getZ().getReal()));

  427.         // Get atmosphere properties in atmosphere own frame
  428.         final Frame      atmFrame  = atmosphere.getFrame();
  429.         final StaticTransform toBody = frame.getStaticTransformTo(atmFrame, date);
  430.         final FieldVector3D<Gradient> posBodyGradient = toBody.transformPosition(position3);
  431.         final FieldAbsoluteDate<Gradient> fieldDate = new FieldAbsoluteDate<>(position3.getX().getField(), date);
  432.         final Gradient density = atmosphere.getDensity(fieldDate, posBodyGradient, atmFrame);

  434.         // Density with derivatives:
  435.         // - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
  436.         // - Others are set to 0.
  437.         final double[] derivatives = Arrays.copyOf(density.getGradient(), position.getX().getFreeParameters());
  438.         return new Gradient(density.getValue(), derivatives);
  439.     }
  440. }
Created with JaCoCo 0.8.12.202403310830