DSSTAtmosphericDrag.java

  1. /* Copyright 2002-2024 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.propagation.semianalytical.dsst.forces;

  18. import java.util.List;
  19. import java.util.stream.Stream;

  20. import org.hipparchus.CalculusFieldElement;
  21. import org.hipparchus.Field;
  22. import org.hipparchus.util.FastMath;
  23. import org.hipparchus.util.MathArrays;
  24. import org.hipparchus.util.MathUtils;
  25. import org.orekit.forces.drag.DragForce;
  26. import org.orekit.forces.drag.DragSensitive;
  27. import org.orekit.forces.drag.IsotropicDrag;
  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.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  34. import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
  35. import org.orekit.utils.Constants;
  36. import org.orekit.utils.ParameterDriver;

  37. /** Atmospheric drag contribution to the
  38.  *  {@link org.orekit.propagation.semianalytical.dsst.DSSTPropagator DSSTPropagator}.
  39.  *  <p>
  40.  *  The drag acceleration is computed through the acceleration model of
  41.  *  {@link org.orekit.forces.drag.DragForce DragForce}.
  42.  *  </p>
  43.  *
  44.  * @author Pascal Parraud
  45.  */
  46. public class DSSTAtmosphericDrag extends AbstractGaussianContribution {

  47.     /** Threshold for the choice of the Gauss quadrature order. */
  48.     private static final double GAUSS_THRESHOLD = 6.0e-10;

  49.     /** Default upper limit for atmospheric drag (m) . */
  50.     private static final double DEFAULT_MAX_ATMOSPHERE_ALTITUDE = 1000000.;

  51.     /** Atmospheric drag force model. */
  52.     private final DragForce drag;

  53.     /** Critical distance from the center of the central body for
  54.      * entering/leaving the atmosphere, i.e. upper limit of atmosphere.
  55.      */
  56.     private double rbar;

  57.     /** Simple constructor with custom force.
  58.      * @param force atmospheric drag force model
  59.      * @param mu central attraction coefficient
  60.      */
  61.     public DSSTAtmosphericDrag(final DragForce force, final double mu) {
  62.         //Call to the constructor from superclass using the numerical drag model as ForceModel
  63.         super("DSST-drag-", GAUSS_THRESHOLD, force, mu);
  64.         this.drag = force;
  65.         this.rbar = DEFAULT_MAX_ATMOSPHERE_ALTITUDE + Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
  66.     }

  67.     /** Simple constructor assuming spherical spacecraft.
  68.      * @param atmosphere atmospheric model
  69.      * @param cd drag coefficient
  70.      * @param area cross sectionnal area of satellite
  71.      * @param mu central attraction coefficient
  72.      */
  73.     public DSSTAtmosphericDrag(final Atmosphere atmosphere, final double cd,
  74.                                final double area, final double mu) {
  75.         this(atmosphere, new IsotropicDrag(area, cd), mu);
  76.     }

  77.     /** Simple constructor with custom spacecraft.
  78.      * @param atmosphere atmospheric model
  79.      * @param spacecraft spacecraft model
  80.      * @param mu central attraction coefficient
  81.      */
  82.     public DSSTAtmosphericDrag(final Atmosphere atmosphere, final DragSensitive spacecraft, final double mu) {

  83.         //Call to the constructor from superclass using the numerical drag model as ForceModel
  84.         this(new DragForce(atmosphere, spacecraft), mu);
  85.     }

  86.     /** Get the atmospheric model.
  87.      * @return atmosphere model
  88.      */
  89.     public Atmosphere getAtmosphere() {
  90.         return drag.getAtmosphere();
  91.     }

  92.     /** Get the critical distance.
  93.      *  <p>
  94.      *  The critical distance from the center of the central body aims at
  95.      *  defining the atmosphere entry/exit.
  96.      *  </p>
  97.      *  @return the critical distance from the center of the central body (m)
  98.      */
  99.     public double getRbar() {
  100.         return rbar;
  101.     }

  102.     /** Set the critical distance from the center of the central body at which
  103.      * the atmosphere is considered to end, i.e. beyond this distance
  104.      * atmospheric drag is not considered.
  105.      *
  106.      *  @param rbar the critical distance from the center of the central body (m)
  107.      */
  108.     public void setRbar(final double rbar) {
  109.         this.rbar = rbar;
  110.     }

  111.     /** {@inheritDoc} */
  112.     public Stream<EventDetector> getEventDetectors() {
  113.         return drag.getEventDetectors();
  114.     }

  115.     /** {@inheritDoc} */
  116.     @Override
  117.     public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventDetectors(final Field<T> field) {
  118.         return drag.getFieldEventDetectors(field);
  119.     }

  120.     /** {@inheritDoc} */
  121.     protected double[] getLLimits(final SpacecraftState state, final AuxiliaryElements auxiliaryElements) {

  122.         final double perigee = auxiliaryElements.getSma() * (1. - auxiliaryElements.getEcc());

  123.         // Trajectory entirely out of the atmosphere
  124.         if (perigee > rbar) {
  125.             return new double[2];
  126.         }
  127.         final double apogee  = auxiliaryElements.getSma() * (1. + auxiliaryElements.getEcc());
  128.         // Trajectory entirely within of the atmosphere
  129.         if (apogee < rbar) {
  130.             return new double[] { -FastMath.PI + MathUtils.normalizeAngle(state.getLv(), 0),
  131.                                   FastMath.PI + MathUtils.normalizeAngle(state.getLv(), 0) };
  132.         }
  133.         // Else, trajectory partialy within of the atmosphere
  134.         final double fb = FastMath.acos(((auxiliaryElements.getSma() * (1. - auxiliaryElements.getEcc() * auxiliaryElements.getEcc()) / rbar) - 1.) / auxiliaryElements.getEcc());
  135.         final double wW = FastMath.atan2(auxiliaryElements.getH(), auxiliaryElements.getK());
  136.         return new double[] {wW - fb, wW + fb};
  137.     }

  138.     /** {@inheritDoc} */
  139.     protected <T extends CalculusFieldElement<T>> T[] getLLimits(final FieldSpacecraftState<T> state,
  140.                                                              final FieldAuxiliaryElements<T> auxiliaryElements) {

  141.         final Field<T> field = state.getDate().getField();

  142.         final T[] tab = MathArrays.buildArray(field, 2);

  143.         final T perigee = auxiliaryElements.getSma().multiply(auxiliaryElements.getEcc().negate().add(1.));
  144.         // Trajectory entirely out of the atmosphere
  145.         if (perigee.getReal() > rbar) {
  146.             return tab;
  147.         }
  148.         final T apogee  = auxiliaryElements.getSma().multiply(auxiliaryElements.getEcc().add(1.));
  149.         // Trajectory entirely within of the atmosphere
  150.         if (apogee.getReal() < rbar) {
  151.             final T zero = field.getZero();
  152.             final T pi   = zero.getPi();
  153.             tab[0] = MathUtils.normalizeAngle(state.getLv(), zero).subtract(pi);
  154.             tab[1] = MathUtils.normalizeAngle(state.getLv(), zero).add(pi);
  155.             return tab;
  156.         }
  157.         // Else, trajectory partialy within of the atmosphere
  158.         final T fb = FastMath.acos(((auxiliaryElements.getSma().multiply(auxiliaryElements.getEcc().multiply(auxiliaryElements.getEcc()).negate().add(1.)).divide(rbar)).subtract(1.)).divide(auxiliaryElements.getEcc()));
  159.         final T wW = FastMath.atan2(auxiliaryElements.getH(), auxiliaryElements.getK());

  160.         tab[0] = wW.subtract(fb);
  161.         tab[1] = wW.add(fb);
  162.         return tab;
  163.     }

  164.     /** {@inheritDoc} */
  165.     @Override
  166.     protected List<ParameterDriver> getParametersDriversWithoutMu() {
  167.         return drag.getParametersDrivers();
  168.     }

  169.     /** Get spacecraft shape.
  170.      *
  171.      * @return spacecraft shape
  172.      */
  173.     public DragSensitive getSpacecraft() {
  174.         return drag.getSpacecraft();
  175.     }

  176.     /** Get drag force.
  177.      *
  178.      * @return drag force
  179.      */
  180.     public DragForce getDrag() {
  181.         return drag;
  182.     }
  183. }