SolarRadiationPressure.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.radiation;

  19. import java.lang.reflect.Array;
  20. import java.util.List;

  22. import org.hipparchus.CalculusFieldElement;
  23. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  24. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  25. import org.hipparchus.util.FastMath;
  26. import org.hipparchus.util.Precision;
  27. import org.orekit.bodies.OneAxisEllipsoid;
  28. import org.orekit.frames.Frame;
  29. import org.orekit.propagation.FieldSpacecraftState;
  30. import org.orekit.propagation.SpacecraftState;
  31. import org.orekit.propagation.events.EventDetectionSettings;
  32. import org.orekit.time.AbsoluteDate;
  33. import org.orekit.time.FieldAbsoluteDate;
  34. import org.orekit.utils.Constants;
  35. import org.orekit.utils.ExtendedPositionProvider;
  36. import org.orekit.utils.FrameAdapter;
  37. import org.orekit.utils.OccultationEngine;
  38. import org.orekit.utils.ParameterDriver;

  40. /** Solar radiation pressure force model.
  41.  * <p>
  42.  * Since Orekit 11.0, it is possible to take into account
  43.  * the eclipses generated by Moon in the solar radiation
  44.  * pressure force model using the
  45.  * {@link #addOccultingBody(ExtendedPositionProvider, double)}
  46.  * method.
  47.  * <p>
  48.  * Example:<br>
  49.  * <code> SolarRadiationPressure srp = </code>
  50.  * <code>                      new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,</code>
  51.  * <code>                                     new IsotropicRadiationClassicalConvention(50.0, 0.5, 0.5));</code><br>
  52.  * <code> srp.addOccultingBody(CelestialBodyFactory.getMoon(), Constants.MOON_EQUATORIAL_RADIUS);</code><br>
  53.  *
  54.  * @author Fabien Maussion
  55.  * @author &Eacute;douard Delente
  56.  * @author V&eacute;ronique Pommier-Maurussane
  57.  * @author Pascal Parraud
  58.  */
  59. public class SolarRadiationPressure extends AbstractRadiationForceModel {

  61.     /** Reference distance for the solar radiation pressure (m). */
  62.     private static final double D_REF = 149597870000.0;

  64.     /** Reference solar radiation pressure at D_REF (N/m²). */
  65.     private static final double P_REF = 4.56e-6;

  67.     /** Margin to force recompute lighting ratio derivatives when we are really inside penumbra. */
  68.     private static final double ANGULAR_MARGIN = 1.0e-10;

  70.     /** Threshold to decide whether the S/C frame is Sun-centered. */
  71.     private static final double SUN_CENTERED_FRAME_THRESHOLD = 2. * Constants.SUN_RADIUS;

  73.     /** Reference flux normalized for a 1m distance (N). */
  74.     private final double kRef;

  76.     /** Sun model. */
  77.     private final ExtendedPositionProvider sun;

  79.     /** Spacecraft. */
  80.     private final RadiationSensitive spacecraft;

  82.     /** Simple constructor with default reference values.
  83.      * <p>When this constructor is used, the reference values are:</p>
  84.      * <ul>
  85.      *   <li>d<sub>ref</sub> = 149597870000.0 m</li>
  86.      *   <li>p<sub>ref</sub> = 4.56 10<sup>-6</sup> N/m²</li>
  87.      * </ul>
  88.      * @param sun Sun model
  89.      * @param centralBody central body shape model (for umbra/penumbra computation)
  90.      * @param spacecraft the object physical and geometrical information
  91.      * @since 12.0
  92.      */
  93.     public SolarRadiationPressure(final ExtendedPositionProvider sun,
  94.                                   final OneAxisEllipsoid centralBody,
  95.                                   final RadiationSensitive spacecraft) {
  96.         this(D_REF, P_REF, sun, centralBody, spacecraft);
  97.     }

  99.     /** Constructor with default eclipse detection settings.
  100.      * <p>Note that reference solar radiation pressure <code>pRef</code> in
  101.      * N/m² is linked to solar flux SF in W/m² using
  102.      * formula pRef = SF/c where c is the speed of light (299792458 m/s). So
  103.      * at 1UA a 1367 W/m² solar flux is a 4.56 10<sup>-6</sup>
  104.      * N/m² solar radiation pressure.</p>
  105.      * @param dRef reference distance for the solar radiation pressure (m)
  106.      * @param pRef reference solar radiation pressure at dRef (N/m²)
  107.      * @param sun Sun model
  108.      * @param centralBody central body shape model (for umbra/penumbra computation)
  109.      * @param spacecraft the object physical and geometrical information
  110.      * @since 12.0
  111.      */
  112.     public SolarRadiationPressure(final double dRef, final double pRef,
  113.                                   final ExtendedPositionProvider sun,
  114.                                   final OneAxisEllipsoid centralBody,
  115.                                   final RadiationSensitive spacecraft) {
  116.         this(dRef, pRef, sun, centralBody, spacecraft, getDefaultEclipseDetectionSettings());
  117.     }

  119.     /** Complete constructor.
  120.      * <p>Note that reference solar radiation pressure <code>pRef</code> in
  121.      * N/m² is linked to solar flux SF in W/m² using
  122.      * formula pRef = SF/c where c is the speed of light (299792458 m/s). So
  123.      * at 1UA a 1367 W/m² solar flux is a 4.56 10<sup>-6</sup>
  124.      * N/m² solar radiation pressure.</p>
  125.      * @param dRef reference distance for the solar radiation pressure (m)
  126.      * @param pRef reference solar radiation pressure at dRef (N/m²)
  127.      * @param sun Sun model
  128.      * @param centralBody central body shape model (for umbra/penumbra computation)
  129.      * @param spacecraft the object physical and geometrical information
  130.      * @param eclipseDetectionSettings event detection settings for penumbra and umbra
  131.      * @since 13.0
  132.      */
  133.     public SolarRadiationPressure(final double dRef, final double pRef,
  134.                                   final ExtendedPositionProvider sun,
  135.                                   final OneAxisEllipsoid centralBody,
  136.                                   final RadiationSensitive spacecraft,
  137.                                   final EventDetectionSettings eclipseDetectionSettings) {
  138.         super(sun, centralBody, eclipseDetectionSettings);
  139.         this.kRef = pRef * dRef * dRef;
  140.         this.sun  = sun;
  141.         this.spacecraft = spacecraft;
  142.     }

  144.     /**
  145.      * Getter for radiation-sensitive spacecraft.
  146.      * @return radiation-sensitive model
  147.      * @since 12.1
  148.      */
  149.     public RadiationSensitive getRadiationSensitiveSpacecraft() {
  150.         return spacecraft;
  151.     }

  153.     /** {@inheritDoc} */
  154.     @Override
  155.     public boolean dependsOnPositionOnly() {
  156.         return spacecraft instanceof IsotropicRadiationClassicalConvention || spacecraft instanceof IsotropicRadiationCNES95Convention || spacecraft instanceof IsotropicRadiationSingleCoefficient;
  157.     }

  159.     /** {@inheritDoc} */
  160.     @Override
  161.     public Vector3D acceleration(final SpacecraftState s, final double[] parameters) {

  163.         final AbsoluteDate date         = s.getDate();
  164.         final Frame        frame        = s.getFrame();
  165.         final Vector3D     position     = s.getPosition();
  166.         final Vector3D     sunPosition  = sun.getPosition(date, frame);
  167.         final Vector3D     sunSatVector = position.subtract(sunPosition);
  168.         final double       r2           = sunSatVector.getNormSq();

  170.         // compute flux
  171.         final double   ratio = getLightingRatio(s, sunPosition);
  172.         final double   rawP  = ratio  * kRef / r2;
  173.         final Vector3D flux  = new Vector3D(rawP / FastMath.sqrt(r2), sunSatVector);

  175.         return spacecraft.radiationPressureAcceleration(s, flux, parameters);

  177.     }

  179.     /** {@inheritDoc} */
  180.     @Override
  181.     public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
  182.                                                                              final T[] parameters) {

  184.         final FieldAbsoluteDate<T> date         = s.getDate();
  185.         final Frame                frame        = s.getFrame();
  186.         final FieldVector3D<T>     position     = s.getPosition();
  187.         final FieldVector3D<T>     sunPosition  = sun.getPosition(date, frame);
  188.         final FieldVector3D<T>     sunSatVector = position.subtract(sunPosition);
  189.         final T                    r2           = sunSatVector.getNormSq();

  191.         // compute flux
  192.         final T                ratio = getLightingRatio(s, sunPosition);
  193.         final T                rawP  = ratio.multiply(kRef).divide(r2);
  194.         final FieldVector3D<T> flux  = new FieldVector3D<>(rawP.divide(r2.sqrt()), sunSatVector);

  196.         return spacecraft.radiationPressureAcceleration(s, flux, parameters);

  198.     }

  200.     /** Check whether the frame is considerer Sun-centered.
  201.      *
  202.      * @param sunPositionInFrame Sun position in frame to test
  203.      * @return true if frame is considered Sun-centered
  204.      * @since 12.0
  205.      */
  206.     private boolean isSunCenteredFrame(final Vector3D sunPositionInFrame) {
  207.         // Frame is considered Sun-centered if Sun (or Solar System barycenter) position
  208.         // in that frame is smaller than SUN_CENTERED_FRAME_THRESHOLD
  209.         return sunPositionInFrame.getNorm() < SUN_CENTERED_FRAME_THRESHOLD;
  210.     }


  213.     /** Get the lighting ratio ([0-1]).
  214.      * @param state spacecraft state
  215.      * @param sunPosition Sun position in S/C frame at S/C date
  216.      * @return lighting ratio
  217.      * @since 12.0 added to avoid numerous call to sun.getPosition(...)
  218.      */
  219.     private double getLightingRatio(final SpacecraftState state, final Vector3D sunPosition) {

  221.         // Check if S/C frame is Sun-centered
  222.         if (isSunCenteredFrame(sunPosition)) {
  223.             // We are in fact computing a trajectory around Sun (or solar system barycenter),
  224.             // not around a planet, we consider lighting ratio will always be 1
  225.             return 1.0;
  226.         }

  228.         final List<OccultationEngine> occultingBodies = getOccultingBodies();
  229.         final int n = occultingBodies.size();

  231.         final OccultationEngine.OccultationAngles[] angles = new OccultationEngine.OccultationAngles[n];
  232.         for (int i = 0; i < n; ++i) {
  233.             angles[i] = occultingBodies.get(i).angles(state);
  234.         }
  235.         final double alphaSunSq = angles[0].getOccultedApparentRadius() * angles[0].getOccultedApparentRadius();

  237.         double result = 0.0;
  238.         for (int i = 0; i < n; ++i) {

  240.             // compute lighting ratio considering one occulting body only
  241.             final OccultationEngine oi  = occultingBodies.get(i);
  242.             final double lightingRatioI = maskingRatio(angles[i]);
  243.             if (lightingRatioI == 0.0) {
  244.                 // body totally occults Sun, total eclipse is occurring.
  245.                 return 0.0;
  246.             }
  247.             result += lightingRatioI;

  249.             // Mutual occulting body eclipse ratio computations between first and secondary bodies
  250.             for (int j = i + 1; j < n; ++j) {

  252.                 final OccultationEngine oj = occultingBodies.get(j);
  253.                 final double lightingRatioJ = maskingRatio(angles[j]);
  254.                 if (lightingRatioJ == 0.0) {
  255.                     // Secondary body totally occults Sun, no more computations are required, total eclipse is occurring.
  256.                     return 0.0;
  257.                 } else if (lightingRatioJ != 1) {
  258.                     // Secondary body partially occults Sun

  260.                     final OccultationEngine oij = new OccultationEngine(new FrameAdapter(oi.getOcculting().getBodyFrame()),
  261.                                                                         oi.getOcculting().getEquatorialRadius(),
  262.                                                                         oj.getOcculting());
  263.                     final OccultationEngine.OccultationAngles aij = oij.angles(state);
  264.                     final double maskingRatioIJ = maskingRatio(aij);
  265.                     final double alphaJSq       = aij.getOccultedApparentRadius() * aij.getOccultedApparentRadius();

  267.                     final double mutualEclipseCorrection = (1 - maskingRatioIJ) * alphaJSq / alphaSunSq;
  268.                     result -= mutualEclipseCorrection;

  270.                 }

  272.             }
  273.         }

  275.         // Final term
  276.         result -= n - 1;

  278.         return result;
  279.     }

  281.     /** Get the lighting ratio ([0-1]).
  282.      * @param state spacecraft state
  283.      * @return lighting ratio
  284.      * @since 7.1
  285.      */
  286.     public double getLightingRatio(final SpacecraftState state) {
  287.         return getLightingRatio(state, sun.getPosition(state.getDate(), state.getFrame()));
  288.     }

  290.     /** Get the masking ratio ([0-1]) considering one pair of bodies.
  291.      * @param angles occultation angles
  292.      * @return masking ratio: 0.0 body fully masked, 1.0 body fully visible
  293.      * @since 12.0
  294.      */
  295.     private double maskingRatio(final OccultationEngine.OccultationAngles angles) {

  297.         // Sat-Occulted/ Sat-Occulting angle
  298.         final double sunSatCentralBodyAngle = angles.getSeparation();

  300.         // Occulting apparent radius
  301.         final double alphaCentral = angles.getLimbRadius();

  303.         // Occulted apparent radius
  304.         final double alphaSun = angles.getOccultedApparentRadius();

  306.         // Is the satellite in complete umbra ?
  307.         if (sunSatCentralBodyAngle - alphaCentral + alphaSun <= ANGULAR_MARGIN) {
  308.             return 0.0;
  309.         } else if (sunSatCentralBodyAngle - alphaCentral - alphaSun < -ANGULAR_MARGIN) {
  310.             // Compute a masking ratio in penumbra
  311.             final double sEA2    = sunSatCentralBodyAngle * sunSatCentralBodyAngle;
  312.             final double oo2sEA  = 1.0 / (2. * sunSatCentralBodyAngle);
  313.             final double aS2     = alphaSun * alphaSun;
  314.             final double aE2     = alphaCentral * alphaCentral;
  315.             final double aE2maS2 = aE2 - aS2;

  317.             final double alpha1  = (sEA2 - aE2maS2) * oo2sEA;
  318.             final double alpha2  = (sEA2 + aE2maS2) * oo2sEA;

  320.             // Protection against numerical inaccuracy at boundaries
  321.             final double almost0 = Precision.SAFE_MIN;
  322.             final double almost1 = FastMath.nextDown(1.0);
  323.             final double a1oaS   = FastMath.min(almost1, FastMath.max(-almost1, alpha1 / alphaSun));
  324.             final double aS2ma12 = FastMath.max(almost0, aS2 - alpha1 * alpha1);
  325.             final double a2oaE   = FastMath.min(almost1, FastMath.max(-almost1, alpha2 / alphaCentral));
  326.             final double aE2ma22 = FastMath.max(almost0, aE2 - alpha2 * alpha2);

  328.             final double P1 = aS2 * FastMath.acos(a1oaS) - alpha1 * FastMath.sqrt(aS2ma12);
  329.             final double P2 = aE2 * FastMath.acos(a2oaE) - alpha2 * FastMath.sqrt(aE2ma22);

  331.             return 1. - (P1 + P2) / (FastMath.PI * aS2);
  332.         } else {
  333.             return 1.0;
  334.         }

  336.     }

  338.     /** Get the lighting ratio ([0-1]).
  339.      * @param state spacecraft state
  340.      * @param sunPosition Sun position in S/C frame at S/C date
  341.      * @param <T> extends CalculusFieldElement
  342.      * @return lighting ratio
  343.      * @since 12.0 added to avoid numerous call to sun.getPosition(...)
  344.      */
  345.     private <T extends CalculusFieldElement<T>> T getLightingRatio(final FieldSpacecraftState<T> state, final FieldVector3D<T> sunPosition) {

  347.         final T one  = state.getDate().getField().getOne();
  348.         if (isSunCenteredFrame(sunPosition.toVector3D())) {
  349.             // We are in fact computing a trajectory around Sun (or solar system barycenter),
  350.             // not around a planet, we consider lighting ratio will always be 1
  351.             return one;
  352.         }
  353.         final T zero = state.getDate().getField().getZero();
  354.         final List<OccultationEngine> occultingBodies = getOccultingBodies();
  355.         final int n = occultingBodies.size();

  357.         @SuppressWarnings("unchecked")
  358.         final OccultationEngine.FieldOccultationAngles<T>[] angles =
  359.             (OccultationEngine.FieldOccultationAngles<T>[]) Array.newInstance(OccultationEngine.FieldOccultationAngles.class, n);
  360.         for (int i = 0; i < n; ++i) {
  361.             angles[i] = occultingBodies.get(i).angles(state);
  362.         }
  363.         final T alphaSunSq = angles[0].getOccultedApparentRadius().multiply(angles[0].getOccultedApparentRadius());

  365.         T result = state.getDate().getField().getZero();
  366.         for (int i = 0; i < n; ++i) {

  368.             // compute lighting ratio considering one occulting body only
  369.             final OccultationEngine oi  = occultingBodies.get(i);
  370.             final T lightingRatioI = maskingRatio(angles[i]);
  371.             if (lightingRatioI.isZero()) {
  372.                 // body totally occults Sun, total eclipse is occurring.
  373.                 return zero;
  374.             }
  375.             result = result.add(lightingRatioI);

  377.             // Mutual occulting body eclipse ratio computations between first and secondary bodies
  378.             for (int j = i + 1; j < n; ++j) {

  380.                 final OccultationEngine oj = occultingBodies.get(j);
  381.                 final T lightingRatioJ = maskingRatio(angles[j]);
  382.                 if (lightingRatioJ.isZero()) {
  383.                     // Secondary body totally occults Sun, no more computations are required, total eclipse is occurring.
  384.                     return zero;
  385.                 } else if (lightingRatioJ.getReal() != 1) {
  386.                     // Secondary body partially occults Sun

  388.                     final OccultationEngine oij = new OccultationEngine(new FrameAdapter(oi.getOcculting().getBodyFrame()),
  389.                                                                         oi.getOcculting().getEquatorialRadius(),
  390.                                                                         oj.getOcculting());
  391.                     final OccultationEngine.FieldOccultationAngles<T> aij = oij.angles(state);
  392.                     final T maskingRatioIJ = maskingRatio(aij);
  393.                     final T alphaJSq       = aij.getOccultedApparentRadius().multiply(aij.getOccultedApparentRadius());

  395.                     final T mutualEclipseCorrection = one.subtract(maskingRatioIJ).multiply(alphaJSq).divide(alphaSunSq);
  396.                     result = result.subtract(mutualEclipseCorrection);

  398.                 }

  400.             }
  401.         }

  403.         // Final term
  404.         result = result.subtract(n - 1);

  406.         return result;
  407.     }

  409.     /** Get the lighting ratio ([0-1]).
  410.      * @param state spacecraft state
  411.      * @param <T> extends CalculusFieldElement
  412.      * @return lighting ratio
  413.      * @since 7.1
  414.      */
  415.     public <T extends CalculusFieldElement<T>> T getLightingRatio(final FieldSpacecraftState<T> state) {
  416.         return getLightingRatio(state, sun.getPosition(state.getDate(), state.getFrame()));
  417.     }

  419.     /** Get the masking ratio ([0-1]) considering one pair of bodies.
  420.      * @param angles occultation angles
  421.      * @param <T> type of the field elements
  422.      * @return masking ratio: 0.0 body fully masked, 1.0 body fully visible
  423.      * @since 12.0
  424.      */
  425.     private <T extends CalculusFieldElement<T>> T maskingRatio(final OccultationEngine.FieldOccultationAngles<T> angles) {


  428.         // Sat-Occulted/ Sat-Occulting angle
  429.         final T occultedSatOcculting = angles.getSeparation();

  431.         // Occulting apparent radius
  432.         final T alphaOcculting = angles.getLimbRadius();

  434.         // Occulted apparent radius
  435.         final T alphaOcculted = angles.getOccultedApparentRadius();

  437.         // Is the satellite in complete umbra ?
  438.         if (occultedSatOcculting.getReal() - alphaOcculting.getReal() + alphaOcculted.getReal() <= ANGULAR_MARGIN) {
  439.             return occultedSatOcculting.getField().getZero();
  440.         } else if (occultedSatOcculting.getReal() - alphaOcculting.getReal() - alphaOcculted.getReal() < -ANGULAR_MARGIN) {
  441.             // Compute a masking ratio in penumbra
  442.             final T sEA2    = occultedSatOcculting.multiply(occultedSatOcculting);
  443.             final T oo2sEA  = occultedSatOcculting.multiply(2).reciprocal();
  444.             final T aS2     = alphaOcculted.multiply(alphaOcculted);
  445.             final T aE2     = alphaOcculting.multiply(alphaOcculting);
  446.             final T aE2maS2 = aE2.subtract(aS2);

  448.             final T alpha1  = sEA2.subtract(aE2maS2).multiply(oo2sEA);
  449.             final T alpha2  = sEA2.add(aE2maS2).multiply(oo2sEA);

  451.             // Protection against numerical inaccuracy at boundaries
  452.             final double almost0 = Precision.SAFE_MIN;
  453.             final double almost1 = FastMath.nextDown(1.0);
  454.             final T a1oaS   = min(almost1, max(-almost1, alpha1.divide(alphaOcculted)));
  455.             final T aS2ma12 = max(almost0, aS2.subtract(alpha1.multiply(alpha1)));
  456.             final T a2oaE   = min(almost1, max(-almost1, alpha2.divide(alphaOcculting)));
  457.             final T aE2ma22 = max(almost0, aE2.subtract(alpha2.multiply(alpha2)));

  459.             final T P1 = aS2.multiply(a1oaS.acos()).subtract(alpha1.multiply(aS2ma12.sqrt()));
  460.             final T P2 = aE2.multiply(a2oaE.acos()).subtract(alpha2.multiply(aE2ma22.sqrt()));

  462.             return occultedSatOcculting.getField().getOne().subtract(P1.add(P2).divide(aS2.multiply(occultedSatOcculting.getPi())));
  463.         } else {
  464.             return occultedSatOcculting.getField().getOne();
  465.         }

  467.     }

  469.     /** {@inheritDoc} */
  470.     @Override
  471.     public List<ParameterDriver> getParametersDrivers() {
  472.         return spacecraft.getRadiationParametersDrivers();
  473.     }

  475.     /** Compute min of two values, one double and one field element.
  476.      * @param d double value
  477.      * @param f field element
  478.      * @param <T> type of the field elements
  479.      * @return min value
  480.      */
  481.     private <T extends CalculusFieldElement<T>> T min(final double d, final T f) {
  482.         return (f.getReal() > d) ? f.getField().getZero().newInstance(d) : f;
  483.     }

  485.     /** Compute max of two values, one double and one field element.
  486.      * @param d double value
  487.      * @param f field element
  488.      * @param <T> type of the field elements
  489.      * @return max value
  490.      */
  491.     private <T extends CalculusFieldElement<T>> T max(final double d, final T f) {
  492.         return (f.getReal() <= d) ? f.getField().getZero().newInstance(d) : f;
  493.     }

  495. }
Created with JaCoCo 0.8.12.202403310830