DSSTAtmosphericDrag.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.propagation.semianalytical.dsst.forces;

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

  22. import org.hipparchus.CalculusFieldElement;
  23. import org.hipparchus.Field;
  24. import org.hipparchus.util.FastMath;
  25. import org.hipparchus.util.MathArrays;
  26. import org.hipparchus.util.MathUtils;
  27. import org.orekit.forces.drag.DragForce;
  28. import org.orekit.forces.drag.DragSensitive;
  29. import org.orekit.forces.drag.IsotropicDrag;
  30. import org.orekit.models.earth.atmosphere.Atmosphere;
  31. import org.orekit.propagation.FieldSpacecraftState;
  32. import org.orekit.propagation.SpacecraftState;
  33. import org.orekit.propagation.events.EventDetector;
  34. import org.orekit.propagation.events.FieldEventDetector;
  35. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  36. import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
  37. import org.orekit.utils.Constants;
  38. import org.orekit.utils.ParameterDriver;

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

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

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

  57.     /** Atmospheric drag force model. */
  58.     private final DragForce drag;

  60.     /** Critical distance from the center of the central body for
  61.      * entering/leaving the atmosphere, i.e. upper limit of atmosphere.
  62.      */
  63.     private double rbar;

  65.     /** Simple constructor with custom force.
  66.      * @param force atmospheric drag force model
  67.      * @param mu central attraction coefficient
  68.      */
  69.     public DSSTAtmosphericDrag(final DragForce force, final double mu) {
  70.         //Call to the constructor from superclass using the numerical drag model as ForceModel
  71.         super("DSST-drag-", GAUSS_THRESHOLD, force, mu);
  72.         this.drag = force;
  73.         this.rbar = DEFAULT_MAX_ATMOSPHERE_ALTITUDE + Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
  74.     }

  76.     /** Simple constructor assuming spherical spacecraft.
  77.      * @param atmosphere atmospheric model
  78.      * @param cd drag coefficient
  79.      * @param area cross sectionnal area of satellite
  80.      * @param mu central attraction coefficient
  81.      */
  82.     public DSSTAtmosphericDrag(final Atmosphere atmosphere, final double cd,
  83.                                final double area, final double mu) {
  84.         this(atmosphere, new IsotropicDrag(area, cd), mu);
  85.     }

  87.     /** Simple constructor with custom spacecraft.
  88.      * @param atmosphere atmospheric model
  89.      * @param spacecraft spacecraft model
  90.      * @param mu central attraction coefficient
  91.      */
  92.     public DSSTAtmosphericDrag(final Atmosphere atmosphere, final DragSensitive spacecraft, final double mu) {

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

  98.     /** Get the atmospheric model.
  99.      * @return atmosphere model
  100.      */
  101.     public Atmosphere getAtmosphere() {
  102.         return drag.getAtmosphere();
  103.     }

  105.     /** Get the critical distance.
  106.      *  <p>
  107.      *  The critical distance from the center of the central body aims at
  108.      *  defining the atmosphere entry/exit.
  109.      *  </p>
  110.      *  @return the critical distance from the center of the central body (m)
  111.      */
  112.     public double getRbar() {
  113.         return rbar;
  114.     }

  116.     /** Set the critical distance from the center of the central body at which
  117.      * the atmosphere is considered to end, i.e. beyond this distance
  118.      * atmospheric drag is not considered.
  119.      *
  120.      *  @param rbar the critical distance from the center of the central body (m)
  121.      */
  122.     public void setRbar(final double rbar) {
  123.         this.rbar = rbar;
  124.     }

  126.     /** {@inheritDoc} */
  127.     public Stream<EventDetector> getEventDetectors() {
  128.         return drag.getEventDetectors();
  129.     }

  131.     /** {@inheritDoc} */
  132.     @Override
  133.     public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventDetectors(final Field<T> field) {
  134.         return drag.getFieldEventDetectors(field);
  135.     }

  137.     /** {@inheritDoc} */
  138.     @Override
  139.     protected double[] getLLimits(final SpacecraftState state, final AuxiliaryElements auxiliaryElements) {

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

  143.         // Trajectory entirely out of the atmosphere
  144.         if (perigee > rbar) {
  145.             return new double[2];
  146.         }
  147.         final double apogee  = auxiliaryElements.getSma() * (1. + auxiliaryElements.getEcc());
  148.         // Trajectory entirely within of the atmosphere
  149.         if (apogee < rbar) {
  150.             return new double[] { -FastMath.PI + MathUtils.normalizeAngle(state.getOrbit().getLv(), 0),
  151.                                   FastMath.PI + MathUtils.normalizeAngle(state.getOrbit().getLv(), 0) };
  152.         }
  153.         // Else, trajectory partialy within of the atmosphere
  154.         final double fb = FastMath.acos(((auxiliaryElements.getSma() * (1. - auxiliaryElements.getEcc() * auxiliaryElements.getEcc()) / rbar) - 1.) / auxiliaryElements.getEcc());
  155.         final double wW = FastMath.atan2(auxiliaryElements.getH(), auxiliaryElements.getK());
  156.         return new double[] {wW - fb, wW + fb};
  157.     }

  159.     /** {@inheritDoc} */
  160.     @Override
  161.     protected <T extends CalculusFieldElement<T>> T[] getLLimits(final FieldSpacecraftState<T> state,
  162.                                                                  final FieldAuxiliaryElements<T> auxiliaryElements) {

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

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

  168.         final T perigee = auxiliaryElements.getSma().multiply(auxiliaryElements.getEcc().negate().add(1.));
  169.         // Trajectory entirely out of the atmosphere
  170.         if (perigee.getReal() > rbar) {
  171.             return tab;
  172.         }
  173.         final T apogee  = auxiliaryElements.getSma().multiply(auxiliaryElements.getEcc().add(1.));
  174.         // Trajectory entirely within of the atmosphere
  175.         if (apogee.getReal() < rbar) {
  176.             final T zero = field.getZero();
  177.             final T pi   = zero.getPi();
  178.             tab[0] = MathUtils.normalizeAngle(state.getOrbit().getLv(), zero).subtract(pi);
  179.             tab[1] = MathUtils.normalizeAngle(state.getOrbit().getLv(), zero).add(pi);
  180.             return tab;
  181.         }
  182.         // Else, trajectory partialy within of the atmosphere
  183.         final T fb = FastMath.acos(((auxiliaryElements.getSma().multiply(auxiliaryElements.getEcc().multiply(auxiliaryElements.getEcc()).negate().add(1.)).divide(rbar)).subtract(1.)).divide(auxiliaryElements.getEcc()));
  184.         final T wW = FastMath.atan2(auxiliaryElements.getH(), auxiliaryElements.getK());

  186.         tab[0] = wW.subtract(fb);
  187.         tab[1] = wW.add(fb);
  188.         return tab;
  189.     }

  191.     /** {@inheritDoc} */
  192.     @Override
  193.     protected List<ParameterDriver> getParametersDriversWithoutMu() {
  194.         return drag.getParametersDrivers();
  195.     }

  197.     /** Get spacecraft shape.
  198.      *
  199.      * @return spacecraft shape
  200.      */
  201.     public DragSensitive getSpacecraft() {
  202.         return drag.getSpacecraft();
  203.     }

  205.     /** Get drag force.
  206.      *
  207.      * @return drag force
  208.      */
  209.     public DragForce getDrag() {
  210.         return drag;
  211.     }
  212. }
Created with JaCoCo 0.8.12.202403310830