AbstractDragForceModel.java

  1. /* Copyright 2002-2020 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.RealFieldElement;
  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.AbstractForceModel;
  26. import org.orekit.frames.Frame;
  27. import org.orekit.frames.Transform;
  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. /**
  35.  * Base class for drag force models.
  36.  * @see DragForce
  37.  * @see TimeSpanDragForce
  38.  * @author Bryan Cazabonne
  39.  * @since 10.2
  40.  */
  41. public abstract class AbstractDragForceModel extends AbstractForceModel {

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

  46.     /**
  47.      * Constructor.
  48.      * @param atmosphere atmospheric model
  49.      */
  50.     protected AbstractDragForceModel(final Atmosphere atmosphere) {
  51.         this.atmosphere = atmosphere;
  52.     }

  54.     /** {@inheritDoc} */
  55.     @Override
  56.     public boolean dependsOnPositionOnly() {
  57.         return false;
  58.     }

  60.     /** Check if a field state corresponds to derivatives with respect to state.
  61.      * @param state state to check
  62.      * @param <T> type of the field elements
  63.      * @return true if state corresponds to derivatives with respect to state
  64.      */
  65.     protected <T extends RealFieldElement<T>> boolean isDSStateDerivative(final FieldSpacecraftState<T> state) {
  66.         try {
  67.             final DerivativeStructure dsMass = (DerivativeStructure) state.getMass();
  68.             final int o = dsMass.getOrder();
  69.             final int p = dsMass.getFreeParameters();

  71.             // To be in the desired case:
  72.             // Order must be 1 (first order derivatives only)
  73.             // Number of parameters must be 6 (PV), 7 (PV + drag coefficient) or 8 (PV + drag coefficient + lift ratio)
  74.             if (o != 1 || (p != 6 && p != 7 && p != 8)) {
  75.                 return false;
  76.             }

  78.             // Check that the first 6 parameters are position and velocity
  79.             @SuppressWarnings("unchecked")
  80.             final FieldPVCoordinates<DerivativeStructure> pv = (FieldPVCoordinates<DerivativeStructure>) state.getPVCoordinates();
  81.             return isVariable(pv.getPosition().getX(), 0) &&
  82.                    isVariable(pv.getPosition().getY(), 1) &&
  83.                    isVariable(pv.getPosition().getZ(), 2) &&
  84.                    isVariable(pv.getVelocity().getX(), 3) &&
  85.                    isVariable(pv.getVelocity().getY(), 4) &&
  86.                    isVariable(pv.getVelocity().getZ(), 5);
  87.         } catch (ClassCastException cce) {
  88.             return false;
  89.         }
  90.     }

  92.     /** Check if a field state corresponds to derivatives with respect to state.
  93.      * @param state state to check
  94.      * @param <T> type of the field elements
  95.      * @return true if state corresponds to derivatives with respect to state
  96.      */
  97.     protected <T extends RealFieldElement<T>> boolean isGradientStateDerivative(final FieldSpacecraftState<T> state) {
  98.         try {
  99.             final Gradient gMass = (Gradient) state.getMass();
  100.             final int p = gMass.getFreeParameters();

  102.             // To be in the desired case:
  103.             // Number of parameters must be 6 (PV), 7 (PV + drag coefficient) or 8 (PV + drag coefficient + lift ratio)
  104.             if (p != 6 && p != 7 && p != 8) {
  105.                 return false;
  106.             }

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

  122.     /** Check if a derivative represents a specified variable.
  123.      * @param ds derivative to check
  124.      * @param index index of the variable
  125.      * @return true if the derivative represents a specified variable
  126.      */
  127.     protected boolean isVariable(final DerivativeStructure ds, final int index) {
  128.         final double[] derivatives = ds.getAllDerivatives();
  129.         boolean check = true;
  130.         for (int i = 1; i < derivatives.length; ++i) {
  131.             check &= derivatives[i] == ((index + 1 == i) ? 1.0 : 0.0);
  132.         }
  133.         return check;
  134.     }

  136.     /** Check if a derivative represents a specified variable.
  137.      * @param g derivative to check
  138.      * @param index index of the variable
  139.      * @return true if the derivative represents a specified variable
  140.      */
  141.     protected boolean isVariable(final Gradient g, final int index) {
  142.         final double[] derivatives = g.getGradient();
  143.         boolean check = true;
  144.         for (int i = 0; i < derivatives.length; ++i) {
  145.             check &= derivatives[i] == ((index == i) ? 1.0 : 0.0);
  146.         }
  147.         return check;
  148.     }

  150.     /** Compute density and its derivatives.
  151.      * Using finite differences for the derivatives.
  152.      * And doing the actual computation only for the derivatives with respect to position (others are set to 0.).
  153.      * <p>
  154.      * From a theoretical point of view, this method computes the same values
  155.      * as {@link Atmosphere#getDensity(FieldAbsoluteDate, FieldVector3D, Frame)} in the
  156.      * specific case of {@link DerivativeStructure} with respect to state, so
  157.      * it is less general. However, it is *much* faster in this important case.
  158.      * </p>
  159.      * <p>
  160.      * The derivatives should be computed with respect to position. The input
  161.      * parameters already take into account the free parameters (6, 7 or 8 depending
  162.      * on derivation with respect to drag coefficient and lift ratio being considered or not)
  163.      * and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
  164.      * with respect to position. Free parameters at indices 3, 4 and 5 correspond
  165.      * to derivatives with respect to velocity (these derivatives will remain zero
  166.      * as the atmospheric density does not depend on velocity). Free parameter
  167.      * at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
  168.      * and/or lift ratio (one of these or both).
  169.      * This 2 last derivatives will remain zero as atmospheric density does not depend on them.
  170.      * </p>
  171.      * @param date current date
  172.      * @param frame inertial reference frame for state (both orbit and attitude)
  173.      * @param position position of spacecraft in inertial frame
  174.      * @return the density and its derivatives
  175.      */
  176.     protected DerivativeStructure getDSDensityWrtStateUsingFiniteDifferences(final AbsoluteDate date,
  177.                                                                              final Frame frame,
  178.                                                                              final FieldVector3D<DerivativeStructure> position) {

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

  186.         // Build a DerivativeStructure using only derivatives with respect to position
  187.         final DSFactory factory3 = new DSFactory(3, 1);
  188.         final FieldVector3D<DerivativeStructure> position3 =
  189.                         new FieldVector3D<>(factory3.variable(0, position.getX().getReal()),
  190.                                             factory3.variable(1,  position.getY().getReal()),
  191.                                             factory3.variable(2,  position.getZ().getReal()));

  193.         // Get atmosphere properties in atmosphere own frame
  194.         final Frame      atmFrame  = atmosphere.getFrame();
  195.         final Transform  toBody    = frame.getTransformTo(atmFrame, date);
  196.         final FieldVector3D<DerivativeStructure> posBodyDS = toBody.transformPosition(position3);
  197.         final Vector3D   posBody   = posBodyDS.toVector3D();

  199.         // Estimate density model by finite differences and composition
  200.         // Using a delta of 1m
  201.         final double delta  = 1.0;
  202.         final double x      = posBody.getX();
  203.         final double y      = posBody.getY();
  204.         final double z      = posBody.getZ();
  205.         final double rho0   = atmosphere.getDensity(date, posBody, atmFrame);
  206.         final double dRhodX = (atmosphere.getDensity(date, new Vector3D(x + delta, y,         z),         atmFrame) - rho0) / delta;
  207.         final double dRhodY = (atmosphere.getDensity(date, new Vector3D(x,         y + delta, z),         atmFrame) - rho0) / delta;
  208.         final double dRhodZ = (atmosphere.getDensity(date, new Vector3D(x,         y,         z + delta), atmFrame) - rho0) / delta;
  209.         final double[] dXdQ = posBodyDS.getX().getAllDerivatives();
  210.         final double[] dYdQ = posBodyDS.getY().getAllDerivatives();
  211.         final double[] dZdQ = posBodyDS.getZ().getAllDerivatives();

  213.         // Density with derivatives:
  214.         // - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
  215.         // - Others are set to 0.
  216.         final int p = factory.getCompiler().getFreeParameters();
  217.         final double[] rhoAll = new double[p + 1];
  218.         rhoAll[0] = rho0;
  219.         for (int i = 1; i < 4; ++i) {
  220.             rhoAll[i] = dRhodX * dXdQ[i] + dRhodY * dYdQ[i] + dRhodZ * dZdQ[i];
  221.         }

  223.         return factory.build(rhoAll);
  224.     }

  226.     /** Compute density and its derivatives.
  227.      * Using finite differences for the derivatives.
  228.      * And doing the actual computation only for the derivatives with respect to position (others are set to 0.).
  229.      * <p>
  230.      * From a theoretical point of view, this method computes the same values
  231.      * as {@link Atmosphere#getDensity(FieldAbsoluteDate, FieldVector3D, Frame)} in the
  232.      * specific case of {@link Gradient} with respect to state, so
  233.      * it is less general. However, it is *much* faster in this important case.
  234.      * </p>
  235.      * <p>
  236.      * The derivatives should be computed with respect to position. The input
  237.      * parameters already take into account the free parameters (6, 7 or 8 depending
  238.      * on derivation with respect to drag coefficient and lift ratio being considered or not)
  239.      * and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
  240.      * with respect to position. Free parameters at indices 3, 4 and 5 correspond
  241.      * to derivatives with respect to velocity (these derivatives will remain zero
  242.      * as the atmospheric density does not depend on velocity). Free parameter
  243.      * at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
  244.      * and/or lift ratio (one of these or both).
  245.      * This 2 last derivatives will remain zero as atmospheric density does not depend on them.
  246.      * </p>
  247.      * @param date current date
  248.      * @param frame inertial reference frame for state (both orbit and attitude)
  249.      * @param position position of spacecraft in inertial frame
  250.      * @return the density and its derivatives
  251.      */
  252.     protected Gradient getGradientDensityWrtStateUsingFiniteDifferences(final AbsoluteDate date,
  253.                                                                         final Frame frame,
  254.                                                                         final FieldVector3D<Gradient> position) {

  256.         // Build a Gradient using only derivatives with respect to position
  257.         final FieldVector3D<Gradient> position3 =
  258.                         new FieldVector3D<>(Gradient.variable(3, 0, position.getX().getReal()),
  259.                                             Gradient.variable(3, 1,  position.getY().getReal()),
  260.                                             Gradient.variable(3, 2,  position.getZ().getReal()));

  262.         // Get atmosphere properties in atmosphere own frame
  263.         final Frame      atmFrame  = atmosphere.getFrame();
  264.         final Transform  toBody    = frame.getTransformTo(atmFrame, date);
  265.         final FieldVector3D<Gradient> posBodyDS = toBody.transformPosition(position3);
  266.         final Vector3D   posBody   = posBodyDS.toVector3D();

  268.         // Estimate density model by finite differences and composition
  269.         // Using a delta of 1m
  270.         final double delta  = 1.0;
  271.         final double x      = posBody.getX();
  272.         final double y      = posBody.getY();
  273.         final double z      = posBody.getZ();
  274.         final double rho0   = atmosphere.getDensity(date, posBody, atmFrame);
  275.         final double dRhodX = (atmosphere.getDensity(date, new Vector3D(x + delta, y,         z),         atmFrame) - rho0) / delta;
  276.         final double dRhodY = (atmosphere.getDensity(date, new Vector3D(x,         y + delta, z),         atmFrame) - rho0) / delta;
  277.         final double dRhodZ = (atmosphere.getDensity(date, new Vector3D(x,         y,         z + delta), atmFrame) - rho0) / delta;
  278.         final double[] dXdQ = posBodyDS.getX().getGradient();
  279.         final double[] dYdQ = posBodyDS.getY().getGradient();
  280.         final double[] dZdQ = posBodyDS.getZ().getGradient();

  282.         // Density with derivatives:
  283.         // - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
  284.         // - Others are set to 0.
  285.         final int p = position.getX().getFreeParameters();
  286.         final double[] rhoAll = new double[p];
  287.         for (int i = 0; i < 3; ++i) {
  288.             rhoAll[i] = dRhodX * dXdQ[i] + dRhodY * dYdQ[i] + dRhodZ * dZdQ[i];
  289.         }

  291.         return new Gradient(rho0, rhoAll);
  292.     }

  294. }
Created with JaCoCo 0.8.6.202009150832