DragForce.java

  1. /* Copyright 2002-2019 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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;

  18. import java.util.stream.Stream;

  19. import org.hipparchus.Field;
  20. import org.hipparchus.RealFieldElement;
  21. import org.hipparchus.analysis.differentiation.DSFactory;
  22. import org.hipparchus.analysis.differentiation.DerivativeStructure;
  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.propagation.SpacecraftState;
  31. import org.orekit.propagation.events.EventDetector;
  32. import org.orekit.propagation.events.FieldEventDetector;
  33. import org.orekit.time.AbsoluteDate;
  34. import org.orekit.time.FieldAbsoluteDate;
  35. import org.orekit.utils.FieldPVCoordinates;
  36. import org.orekit.utils.ParameterDriver;


  37. /** Atmospheric drag force model.
  38.  *
  39.  * The drag acceleration is computed as follows :
  40.  *
  41.  * γ = (1/2 * ρ * V² * S / Mass) * DragCoefVector
  42.  *
  43.  * With DragCoefVector = {C<sub>x</sub>, C<sub>y</sub>, C<sub>z</sub>} and S given by the user through the interface
  44.  * {@link DragSensitive}
  45.  *
  46.  * @author &Eacute;douard Delente
  47.  * @author Fabien Maussion
  48.  * @author V&eacute;ronique Pommier-Maurussane
  49.  * @author Pascal Parraud
  50.  */

  51. public class DragForce extends AbstractForceModel {

  52.     /** Atmospheric model. */
  53.     private final Atmosphere atmosphere;

  54.     /** Spacecraft. */
  55.     private final DragSensitive spacecraft;

  56.     /** Simple constructor.
  57.      * @param atmosphere atmospheric model
  58.      * @param spacecraft the object physical and geometrical information
  59.      */
  60.     public DragForce(final Atmosphere atmosphere, final DragSensitive spacecraft) {
  61.         this.atmosphere = atmosphere;
  62.         this.spacecraft = spacecraft;
  63.     }

  64.     /** {@inheritDoc} */
  65.     @Override
  66.     public boolean dependsOnPositionOnly() {
  67.         return false;
  68.     }

  69.     /** {@inheritDoc} */
  70.     @Override
  71.     public Vector3D acceleration(final SpacecraftState s, final double[] parameters) {

  72.         final AbsoluteDate date     = s.getDate();
  73.         final Frame        frame    = s.getFrame();
  74.         final Vector3D     position = s.getPVCoordinates().getPosition();

  75.         final double rho    = atmosphere.getDensity(date, position, frame);
  76.         final Vector3D vAtm = atmosphere.getVelocity(date, position, frame);
  77.         final Vector3D relativeVelocity = vAtm.subtract(s.getPVCoordinates().getVelocity());

  78.         return spacecraft.dragAcceleration(date, frame, position, s.getAttitude().getRotation(),
  79.                                            s.getMass(), rho, relativeVelocity, parameters);

  80.     }

  81.     /** {@inheritDoc} */
  82.     @Override
  83.     public <T extends RealFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
  84.                                                                          final T[] parameters) {

  85.         final FieldAbsoluteDate<T> date     = s.getDate();
  86.         final Frame                frame    = s.getFrame();
  87.         final FieldVector3D<T>     position = s.getPVCoordinates().getPosition();

  88.         // Density and its derivatives
  89.         final T rho;

  90.         // Check for faster computation dedicated to derivatives with respect to state
  91.         // Using finite differences instead of automatic differentiation as it seems to be much
  92.         // faster for the drag's derivatives' computation
  93.         if (isStateDerivative(s)) {
  94.             rho = this.getDensityWrtStateUsingFiniteDifferences(date.toAbsoluteDate(), frame, position);
  95.         } else {
  96.             rho    = atmosphere.getDensity(date, position, frame);
  97.         }

  98.         // Spacecraft relative velocity with respect to the atmosphere
  99.         final FieldVector3D<T> vAtm = atmosphere.getVelocity(date, position, frame);
  100.         final FieldVector3D<T> relativeVelocity = vAtm.subtract(s.getPVCoordinates().getVelocity());

  101.         // Drag acceleration along with its derivatives
  102.         return spacecraft.dragAcceleration(date, frame, position, s.getAttitude().getRotation(),
  103.                                            s.getMass(), rho, relativeVelocity, parameters);

  104.     }

  105.     /** {@inheritDoc} */
  106.     @Override
  107.     public Stream<EventDetector> getEventsDetectors() {
  108.         return Stream.empty();
  109.     }

  110.     /** {@inheritDoc} */
  111.     @Override
  112.     public <T extends RealFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventsDetectors(final Field<T> field) {
  113.         return Stream.empty();
  114.     }

  115.     /** {@inheritDoc} */
  116.     @Override
  117.     public ParameterDriver[] getParametersDrivers() {
  118.         return spacecraft.getDragParametersDrivers();
  119.     }

  120.     /** Check if a field state corresponds to derivatives with respect to state.
  121.      * @param state state to check
  122.      * @param <T> type of the field elements
  123.      * @return true if state corresponds to derivatives with respect to state
  124.      * @since 9.0
  125.      */
  126.     private <T extends RealFieldElement<T>> boolean isStateDerivative(final FieldSpacecraftState<T> state) {
  127.         try {
  128.             final DerivativeStructure dsMass = (DerivativeStructure) state.getMass();
  129.             final int o = dsMass.getOrder();
  130.             final int p = dsMass.getFreeParameters();

  131.             // To be in the desired case:
  132.             // Order must be 1 (first order derivatives only)
  133.             // Number of parameters must be 6 (PV), 7 (PV + drag coefficient) or 8 (PV + drag coefficient + lift ratio)
  134.             if (o != 1 || (p != 6 && p != 7 && p != 8)) {
  135.                 return false;
  136.             }

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

  150.     /** Check if a derivative represents a specified variable.
  151.      * @param ds derivative to check
  152.      * @param index index of the variable
  153.      * @return true if the derivative represents a specified variable
  154.      * @since 9.0
  155.      */
  156.     private boolean isVariable(final DerivativeStructure ds, final int index) {
  157.         final double[] derivatives = ds.getAllDerivatives();
  158.         boolean check = true;
  159.         for (int i = 1; i < derivatives.length; ++i) {
  160.             check &= derivatives[i] == ((index + 1 == i) ? 1.0 : 0.0);
  161.         }
  162.         return check;
  163.     }

  164.     /** Compute density and its derivatives.
  165.      * Using finite differences for the derivatives.
  166.      * And doing the actual computation only for the derivatives with respect to position (others are set to 0.).
  167.      * <p>
  168.      * From a theoretical point of view, this method computes the same values
  169.      * as {@link Atmosphere#getDensity(FieldAbsoluteDate, FieldVector3D, Frame)} in the
  170.      * specific case of {@link DerivativeStructure} with respect to state, so
  171.      * it is less general. However, it is *much* faster in this important case.
  172.      * <p>
  173.      * <p>
  174.      * The derivatives should be computed with respect to position. The input
  175.      * parameters already take into account the free parameters (6, 7 or 8 depending
  176.      * on derivation with respect to drag coefficient and lift ratio being considered or not)
  177.      * and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
  178.      * with respect to position. Free parameters at indices 3, 4 and 5 correspond
  179.      * to derivatives with respect to velocity (these derivatives will remain zero
  180.      * as the atmospheric density does not depend on velocity). Free parameter
  181.      * at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
  182.      * and/or lift ratio (one of these or both).
  183.      * This 2 last derivatives will remain zero as atmospheric density does not depend on them.
  184.      * </p>
  185.      * @param date current date
  186.      * @param frame inertial reference frame for state (both orbit and attitude)
  187.      * @param position position of spacecraft in inertial frame
  188.      * @param <T> type of the elements
  189.      * @return the density and its derivatives
  190.           * @since 9.0
  191.      */
  192.     private <T extends RealFieldElement<T>> T getDensityWrtStateUsingFiniteDifferences(final AbsoluteDate date,
  193.                                                                                        final Frame frame,
  194.                                                                                        final FieldVector3D<T> position) {

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

  200.         // Build a DerivativeStructure using only derivatives with respect to position
  201.         final DSFactory factory3 = new DSFactory(3, 1);
  202.         final FieldVector3D<DerivativeStructure> position3 =
  203.                         new FieldVector3D<>(factory3.variable(0, position.getX().getReal()),
  204.                                             factory3.variable(1,  position.getY().getReal()),
  205.                                             factory3.variable(2,  position.getZ().getReal()));

  206.         // Get atmosphere properties in atmosphere own frame
  207.         final Frame      atmFrame  = atmosphere.getFrame();
  208.         final Transform  toBody    = frame.getTransformTo(atmFrame, date);
  209.         final FieldVector3D<DerivativeStructure> posBodyDS = toBody.transformPosition(position3);
  210.         final Vector3D   posBody   = posBodyDS.toVector3D();

  211.         // Estimate density model by finite differences and composition
  212.         // Using a delta of 1m
  213.         final double delta  = 1.0;
  214.         final double x      = posBody.getX();
  215.         final double y      = posBody.getY();
  216.         final double z      = posBody.getZ();
  217.         final double rho0   = atmosphere.getDensity(date, posBody, atmFrame);
  218.         final double dRhodX = (atmosphere.getDensity(date, new Vector3D(x + delta, y,         z),         atmFrame) - rho0) / delta;
  219.         final double dRhodY = (atmosphere.getDensity(date, new Vector3D(x,         y + delta, z),         atmFrame) - rho0) / delta;
  220.         final double dRhodZ = (atmosphere.getDensity(date, new Vector3D(x,         y,         z + delta), atmFrame) - rho0) / delta;
  221.         final double[] dXdQ = posBodyDS.getX().getAllDerivatives();
  222.         final double[] dYdQ = posBodyDS.getY().getAllDerivatives();
  223.         final double[] dZdQ = posBodyDS.getZ().getAllDerivatives();

  224.         // Density with derivatives:
  225.         // - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
  226.         // - Others are set to 0.
  227.         final int p = factory.getCompiler().getFreeParameters();
  228.         final double[] rhoAll = new double[p + 1];
  229.         rhoAll[0] = rho0;
  230.         for (int i = 1; i < 4; ++i) {
  231.             rhoAll[i] = dRhodX * dXdQ[i] + dRhodY * dYdQ[i] + dRhodZ * dZdQ[i];
  232.         }
  233.         @SuppressWarnings("unchecked")
  234.         final T rho = (T) (factory.build(rhoAll));

  235.         return rho;
  236.     }
  237. }