SolarRadiationPressure.java

  1. /* Copyright 2002-2022 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;

  18. import java.util.ArrayList;
  19. import java.util.List;
  20. import java.util.Map;

  21. import org.hipparchus.CalculusFieldElement;
  22. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  23. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  24. import org.hipparchus.util.FastMath;
  25. import org.hipparchus.util.MathArrays;
  26. import org.hipparchus.util.Precision;
  27. import org.orekit.frames.Frame;
  28. import org.orekit.propagation.FieldSpacecraftState;
  29. import org.orekit.propagation.SpacecraftState;
  30. import org.orekit.time.AbsoluteDate;
  31. import org.orekit.time.FieldAbsoluteDate;
  32. import org.orekit.utils.Constants;
  33. import org.orekit.utils.ExtendedPVCoordinatesProvider;
  34. import org.orekit.utils.ParameterDriver;

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

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

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

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

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

  63.     /** Sun model. */
  64.     private final ExtendedPVCoordinatesProvider sun;

  65.     /** Spacecraft. */
  66.     private final RadiationSensitive spacecraft;

  67.     /** Simple constructor with default reference values.
  68.      * <p>When this constructor is used, the reference values are:</p>
  69.      * <ul>
  70.      *   <li>d<sub>ref</sub> = 149597870000.0 m</li>
  71.      *   <li>p<sub>ref</sub> = 4.56 10<sup>-6</sup> N/m²</li>
  72.      * </ul>
  73.      * @param sun Sun model
  74.      * @param equatorialRadius spherical shape model (for umbra/penumbra computation)
  75.      * @param spacecraft the object physical and geometrical information
  76.      * @since 9.2
  77.      */
  78.     public SolarRadiationPressure(final ExtendedPVCoordinatesProvider sun, final double equatorialRadius,
  79.                                   final RadiationSensitive spacecraft) {
  80.         this(D_REF, P_REF, sun, equatorialRadius, spacecraft);
  81.     }

  82.     /** Complete constructor.
  83.      * <p>Note that reference solar radiation pressure <code>pRef</code> in
  84.      * N/m² is linked to solar flux SF in W/m² using
  85.      * formula pRef = SF/c where c is the speed of light (299792458 m/s). So
  86.      * at 1UA a 1367 W/m² solar flux is a 4.56 10<sup>-6</sup>
  87.      * N/m² solar radiation pressure.</p>
  88.      * @param dRef reference distance for the solar radiation pressure (m)
  89.      * @param pRef reference solar radiation pressure at dRef (N/m²)
  90.      * @param sun Sun model
  91.      * @param equatorialRadius spherical shape model (for umbra/penumbra computation)
  92.      * @param spacecraft the object physical and geometrical information
  93.      * @since 9.2
  94.      */
  95.     public SolarRadiationPressure(final double dRef, final double pRef,
  96.                                   final ExtendedPVCoordinatesProvider sun,
  97.                                   final double equatorialRadius,
  98.                                   final RadiationSensitive spacecraft) {
  99.         super(sun, equatorialRadius);
  100.         this.kRef = pRef * dRef * dRef;
  101.         this.sun  = sun;
  102.         this.spacecraft = spacecraft;
  103.     }

  104.     /** {@inheritDoc} */
  105.     @Override
  106.     public Vector3D acceleration(final SpacecraftState s, final double[] parameters) {

  107.         final AbsoluteDate date         = s.getDate();
  108.         final Frame        frame        = s.getFrame();
  109.         final Vector3D     position     = s.getPVCoordinates().getPosition();
  110.         final Vector3D     sunSatVector = position.subtract(sun.getPVCoordinates(date, frame).getPosition());
  111.         final double       r2           = sunSatVector.getNormSq();

  112.         // compute flux
  113.         final double   ratio = getTotalLightingRatio(position, frame, date);
  114.         final double   rawP  = ratio  * kRef / r2;
  115.         final Vector3D flux  = new Vector3D(rawP / FastMath.sqrt(r2), sunSatVector);

  116.         return spacecraft.radiationPressureAcceleration(date, frame, position, s.getAttitude().getRotation(),
  117.                                                         s.getMass(), flux, parameters);

  118.     }

  119.     /** {@inheritDoc} */
  120.     @Override
  121.     public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
  122.                                                                          final T[] parameters) {

  123.         final FieldAbsoluteDate<T> date         = s.getDate();
  124.         final Frame                frame        = s.getFrame();
  125.         final FieldVector3D<T>     position     = s.getPVCoordinates().getPosition();
  126.         final FieldVector3D<T>     sunSatVector = position.subtract(sun.getPVCoordinates(date, frame).getPosition());
  127.         final T                    r2           = sunSatVector.getNormSq();

  128.         // compute flux
  129.         final T                ratio = getTotalLightingRatio(position, frame, date);
  130.         final T                rawP  = ratio.multiply(kRef).divide(r2);
  131.         final FieldVector3D<T> flux  = new FieldVector3D<>(rawP.divide(r2.sqrt()), sunSatVector);

  132.         return spacecraft.radiationPressureAcceleration(date, frame, position, s.getAttitude().getRotation(),
  133.                                                         s.getMass(), flux, parameters);

  134.     }

  135.     /** Get the lighting ratio ([0-1]).
  136.      * Considers only central body as occulting body.
  137.      * @param position the satellite's position in the selected frame.
  138.      * @param frame in which is defined the position
  139.      * @param date the date
  140.      * @return lighting ratio
  141.           * @since 7.1
  142.      */
  143.     public double getLightingRatio(final Vector3D position, final Frame frame, final AbsoluteDate date) {

  144.         final Vector3D sunPosition = sun.getPVCoordinates(date, frame).getPosition();
  145.         if (sunPosition.getNorm() < 2 * Constants.SUN_RADIUS) {
  146.             // we are in fact computing a trajectory around Sun (or solar system barycenter),
  147.             // not around a planet,we consider lighting ratio is always 1
  148.             return 1.0;
  149.         }

  150.         // Compute useful angles
  151.         final double[] angle = getEclipseAngles(sunPosition, position);

  152.         // Sat-Sun / Sat-CentralBody angle
  153.         final double sunSatCentralBodyAngle = angle[0];

  154.         // Central Body apparent radius
  155.         final double alphaCentral = angle[1];

  156.         // Sun apparent radius
  157.         final double alphaSun = angle[2];

  158.         double result = 1.0;

  159.         // Is the satellite in complete umbra ?
  160.         if (sunSatCentralBodyAngle - alphaCentral + alphaSun <= ANGULAR_MARGIN) {
  161.             result = 0.0;
  162.         } else if (sunSatCentralBodyAngle - alphaCentral - alphaSun < -ANGULAR_MARGIN) {
  163.             // Compute a lighting ratio in penumbra
  164.             final double sEA2    = sunSatCentralBodyAngle * sunSatCentralBodyAngle;
  165.             final double oo2sEA  = 1.0 / (2. * sunSatCentralBodyAngle);
  166.             final double aS2     = alphaSun * alphaSun;
  167.             final double aE2     = alphaCentral * alphaCentral;
  168.             final double aE2maS2 = aE2 - aS2;

  169.             final double alpha1  = (sEA2 - aE2maS2) * oo2sEA;
  170.             final double alpha2  = (sEA2 + aE2maS2) * oo2sEA;

  171.             // Protection against numerical inaccuracy at boundaries
  172.             final double almost0 = Precision.SAFE_MIN;
  173.             final double almost1 = FastMath.nextDown(1.0);
  174.             final double a1oaS   = FastMath.min(almost1, FastMath.max(-almost1, alpha1 / alphaSun));
  175.             final double aS2ma12 = FastMath.max(almost0, aS2 - alpha1 * alpha1);
  176.             final double a2oaE   = FastMath.min(almost1, FastMath.max(-almost1, alpha2 / alphaCentral));
  177.             final double aE2ma22 = FastMath.max(almost0, aE2 - alpha2 * alpha2);

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

  180.             result = 1. - (P1 + P2) / (FastMath.PI * aS2);
  181.         }

  182.         return result;

  183.     }

  184.     /** Get eclipse ratio between to bodies seen from a specific object.
  185.      * Ratio is in [0-1].
  186.      * @param position the satellite's position in the selected frame
  187.      * @param occultingPosition the position of the occulting object
  188.      * @param occultingRadius the mean radius of the occulting object
  189.      * @param occultedPosition the position of the occulted object
  190.      * @param occultedRadius the mean radius of the occulted object
  191.      * @return eclipse ratio
  192.      */
  193.     public double getGeneralEclipseRatio(final Vector3D position,
  194.                                   final Vector3D occultingPosition,
  195.                                   final double occultingRadius,
  196.                                   final Vector3D occultedPosition,
  197.                                   final double occultedRadius) {


  198.         // Compute useful angles
  199.         final double[] angle = getGeneralEclipseAngles(position, occultingPosition, occultingRadius, occultedPosition, occultedRadius);

  200.         // Sat-Occulted/ Sat-Occulting angle
  201.         final double occultedSatOcculting = angle[0];

  202.         // Occulting apparent radius
  203.         final double alphaOcculting = angle[1];

  204.         // Occulted apparent radius
  205.         final double alphaOcculted = angle[2];

  206.         double result = 1.0;

  207.         // Is the satellite in complete umbra ?
  208.         if (occultedSatOcculting - alphaOcculting + alphaOcculted <= ANGULAR_MARGIN) {
  209.             result = 0.0;
  210.         } else if (occultedSatOcculting - alphaOcculting - alphaOcculted < -ANGULAR_MARGIN) {
  211.             // Compute an eclipse ratio in penumbra
  212.             final double sEA2    = occultedSatOcculting * occultedSatOcculting;
  213.             final double oo2sEA  = 1.0 / (2. * occultedSatOcculting);
  214.             final double aS2     = alphaOcculted * alphaOcculted;
  215.             final double aE2     = alphaOcculting * alphaOcculting;
  216.             final double aE2maS2 = aE2 - aS2;

  217.             final double alpha1  = (sEA2 - aE2maS2) * oo2sEA;
  218.             final double alpha2  = (sEA2 + aE2maS2) * oo2sEA;

  219.             // Protection against numerical inaccuracy at boundaries
  220.             final double almost0 = Precision.SAFE_MIN;
  221.             final double almost1 = FastMath.nextDown(1.0);
  222.             final double a1oaS   = FastMath.min(almost1, FastMath.max(-almost1, alpha1 / alphaOcculted));
  223.             final double aS2ma12 = FastMath.max(almost0, aS2 - alpha1 * alpha1);
  224.             final double a2oaE   = FastMath.min(almost1, FastMath.max(-almost1, alpha2 / alphaOcculting));
  225.             final double aE2ma22 = FastMath.max(almost0, aE2 - alpha2 * alpha2);

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

  228.             result = 1. - (P1 + P2) / (FastMath.PI * aS2);
  229.         }

  230.         return result;
  231.     }

  232.     /** Get the total lighting ratio ([0-1]).
  233.      * This method considers every occulting bodies.
  234.      * @param position the satellite's position in the selected frame.
  235.      * @param frame in which is defined the position
  236.      * @param date the date
  237.      * @return lighting ratio
  238.      */
  239.     public double getTotalLightingRatio(final Vector3D position, final Frame frame, final AbsoluteDate date) {

  240.         double result = 0.0;
  241.         final Map<ExtendedPVCoordinatesProvider, Double> otherOccultingBodies = getOtherOccultingBodies();
  242.         final Vector3D sunPosition = sun.getPVCoordinates(date, frame).getPosition();
  243.         final int n = otherOccultingBodies.size() + 1;

  244.         if (n > 1) {

  245.             final Vector3D[] occultingBodyPositions = new Vector3D[n];
  246.             final double[] occultingBodyRadiuses = new double[n];

  247.             // Central body
  248.             occultingBodyPositions[0] = new Vector3D(0.0, 0.0, 0.0);
  249.             occultingBodyRadiuses[0] = getEquatorialRadius();

  250.             // Other occulting bodies
  251.             int k = 1;
  252.             for (Map.Entry<ExtendedPVCoordinatesProvider, Double> entry : otherOccultingBodies.entrySet()) {
  253.                 occultingBodyPositions[k] = entry.getKey().getPVCoordinates(date, frame).getPosition();
  254.                 occultingBodyRadiuses[k]  = entry.getValue();
  255.                 ++k;
  256.             }
  257.             for (int i = 0; i < n; ++i) {

  258.                 // Lighting ratio computations
  259.                 final double eclipseRatioI = getGeneralEclipseRatio(position,
  260.                                                                     occultingBodyPositions[i],
  261.                                                                     occultingBodyRadiuses[i],
  262.                                                                     sunPosition,
  263.                                                                     Constants.SUN_RADIUS);

  264.                 // First body totaly occults Sun, full eclipse is occuring.
  265.                 if (eclipseRatioI == 0.0) {
  266.                     return 0.0;
  267.                 }


  268.                 result += eclipseRatioI;


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

  271.                     final double eclipseRatioJ = getGeneralEclipseRatio(position,
  272.                                                                         occultingBodyPositions[j],
  273.                                                                         occultingBodyRadiuses[j],
  274.                                                                         sunPosition,
  275.                                                                         Constants.SUN_RADIUS);
  276.                     final double eclipseRatioIJ = getGeneralEclipseRatio(position,
  277.                                                                          occultingBodyPositions[i],
  278.                                                                          occultingBodyRadiuses[i],
  279.                                                                          occultingBodyPositions[j],
  280.                                                                          occultingBodyRadiuses[j]);

  281.                     final double alphaJ = getGeneralEclipseAngles(position,
  282.                                                                   occultingBodyPositions[i],
  283.                                                                   occultingBodyRadiuses[i],
  284.                                                                   occultingBodyPositions[j],
  285.                                                                   occultingBodyRadiuses[j])[2];

  286.                     final double alphaSun = getEclipseAngles(sunPosition, position)[2];
  287.                     final double alphaJSq = alphaJ * alphaJ;
  288.                     final double alphaSunSq = alphaSun * alphaSun;
  289.                     final double mutualEclipseCorrection = (1 - eclipseRatioIJ) * alphaJSq / alphaSunSq;

  290.                     // Secondary body totally occults Sun, no more computations are required, full eclipse is occuring.
  291.                     if (eclipseRatioJ ==  0.0 ) {
  292.                         return 0.0;
  293.                     }

  294.                     // Secondary body partially occults Sun
  295.                     else if (eclipseRatioJ != 1) {
  296.                         result -= mutualEclipseCorrection;
  297.                     }
  298.                 }
  299.             }
  300.             // Final term
  301.             result -= n - 1;
  302.         } else {
  303.             // only central body is considered
  304.             result = getLightingRatio(position, frame, date);
  305.         }
  306.         return result;
  307.     }


  308.     /** Get the lighting ratio ([0-1]).
  309.      * Considers only central body as occulting body.
  310.      * @param position the satellite's position in the selected frame.
  311.      * @param frame in which is defined the position
  312.      * @param date the date
  313.      * @param <T> extends CalculusFieldElement
  314.      * @return lighting ratio
  315.           * @since 7.1
  316.      */
  317.     public <T extends CalculusFieldElement<T>> T getLightingRatio(final FieldVector3D<T> position,
  318.                                                               final Frame frame,
  319.                                                               final FieldAbsoluteDate<T> date) {

  320.         final T one = date.getField().getOne();

  321.         final FieldVector3D<T> sunPosition = sun.getPVCoordinates(date, frame).getPosition();
  322.         if (sunPosition.getNorm().getReal() < 2 * Constants.SUN_RADIUS) {
  323.             // we are in fact computing a trajectory around Sun (or solar system barycenter),
  324.             // not around a planet,we consider lighting ratio is always 1
  325.             return one;
  326.         }

  327.         // Compute useful angles
  328.         final T[] angle = getEclipseAngles(sunPosition, position);

  329.         // Sat-Sun / Sat-CentralBody angle
  330.         final T sunsatCentralBodyAngle = angle[0];

  331.         // Central Body apparent radius
  332.         final T alphaCentral = angle[1];

  333.         // Sun apparent radius
  334.         final T alphaSun = angle[2];

  335.         T result = one;
  336.         // Is the satellite in complete umbra ?
  337.         if (sunsatCentralBodyAngle.getReal() - alphaCentral.getReal() + alphaSun.getReal() <= ANGULAR_MARGIN) {
  338.             result = date.getField().getZero();
  339.         } else if (sunsatCentralBodyAngle.getReal() - alphaCentral.getReal() - alphaSun.getReal() < -ANGULAR_MARGIN) {
  340.             // Compute a lighting ratio in penumbra
  341.             final T sEA2    = sunsatCentralBodyAngle.multiply(sunsatCentralBodyAngle);
  342.             final T oo2sEA  = sunsatCentralBodyAngle.multiply(2).reciprocal();
  343.             final T aS2     = alphaSun.multiply(alphaSun);
  344.             final T aE2     = alphaCentral.multiply(alphaCentral);
  345.             final T aE2maS2 = aE2.subtract(aS2);

  346.             final T alpha1  = sEA2.subtract(aE2maS2).multiply(oo2sEA);
  347.             final T alpha2  = sEA2.add(aE2maS2).multiply(oo2sEA);

  348.             // Protection against numerical inaccuracy at boundaries
  349.             final double almost0 = Precision.SAFE_MIN;
  350.             final double almost1 = FastMath.nextDown(1.0);
  351.             final T a1oaS   = min(almost1, max(-almost1, alpha1.divide(alphaSun)));
  352.             final T aS2ma12 = max(almost0, aS2.subtract(alpha1.multiply(alpha1)));
  353.             final T a2oaE   = min(almost1, max(-almost1, alpha2.divide(alphaCentral)));
  354.             final T aE2ma22 = max(almost0, aE2.subtract(alpha2.multiply(alpha2)));

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

  357.             result = one.subtract(P1.add(P2).divide(aS2.multiply(one.getPi())));
  358.         }

  359.         return result;
  360.     }

  361.     /** Get eclipse ratio between to bodies seen from a specific object.
  362.      * Ratio is in [0-1].
  363.      * @param position the satellite's position in the selected frame
  364.      * @param occultingPosition the position of the occulting object
  365.      * @param occultingRadius the mean radius of the occulting object
  366.      * @param occultedPosition the position of the occulted object
  367.      * @param occultedRadius the mean radius of the occulted object
  368.      * @param <T> extends RealFieldElement
  369.      * @return eclipse ratio
  370.      */
  371.     public <T extends CalculusFieldElement<T>> T getGeneralEclipseRatio(final FieldVector3D<T> position,
  372.                                                                         final FieldVector3D<T> occultingPosition,
  373.                                                                         final T occultingRadius,
  374.                                                                         final FieldVector3D<T> occultedPosition,
  375.                                                                         final T occultedRadius) {


  376.         final T one = occultingRadius.getField().getOne();

  377.         // Compute useful angles
  378.         final T[] angle = getGeneralEclipseAngles(position, occultingPosition, occultingRadius, occultedPosition, occultedRadius);

  379.         // Sat-Occulted/ Sat-Occulting angle
  380.         final T occultedSatOcculting = angle[0];

  381.         // Occulting apparent radius
  382.         final T alphaOcculting = angle[1];

  383.         // Occulted apparent radius
  384.         final T alphaOcculted = angle[2];

  385.         T result = one;

  386.         // Is the satellite in complete umbra ?
  387.         if (occultedSatOcculting.getReal() - alphaOcculting.getReal() + alphaOcculted.getReal() <= ANGULAR_MARGIN) {
  388.             result = occultingRadius.getField().getZero();
  389.         } else if (occultedSatOcculting.getReal() - alphaOcculting.getReal() - alphaOcculted.getReal() < -ANGULAR_MARGIN) {
  390.             // Compute a lighting ratio in penumbra
  391.             final T sEA2    = occultedSatOcculting.multiply(occultedSatOcculting);
  392.             final T oo2sEA  = occultedSatOcculting.multiply(2).reciprocal();
  393.             final T aS2     = alphaOcculted.multiply(alphaOcculted);
  394.             final T aE2     = alphaOcculting.multiply(alphaOcculting);
  395.             final T aE2maS2 = aE2.subtract(aS2);

  396.             final T alpha1  = sEA2.subtract(aE2maS2).multiply(oo2sEA);
  397.             final T alpha2  = sEA2.add(aE2maS2).multiply(oo2sEA);

  398.             // Protection against numerical inaccuracy at boundaries
  399.             final double almost0 = Precision.SAFE_MIN;
  400.             final double almost1 = FastMath.nextDown(1.0);
  401.             final T a1oaS   = min(almost1, max(-almost1, alpha1.divide(alphaOcculted)));
  402.             final T aS2ma12 = max(almost0, aS2.subtract(alpha1.multiply(alpha1)));
  403.             final T a2oaE   = min(almost1, max(-almost1, alpha2.divide(alphaOcculting)));
  404.             final T aE2ma22 = max(almost0, aE2.subtract(alpha2.multiply(alpha2)));

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

  407.             result = one.subtract(P1.add(P2).divide(aS2.multiply(one.getPi())));
  408.         }

  409.         return result;
  410.     }

  411.     /** Get the total lighting ratio ([0-1]).
  412.      * This method considers every occulting bodies.
  413.      * @param position the satellite's position in the selected frame.
  414.      * @param frame in which is defined the position
  415.      * @param date the date
  416.      * @param <T> extends RealFieldElement
  417.      * @return lighting rati
  418.      */
  419.     public  <T extends CalculusFieldElement<T>> T getTotalLightingRatio(final FieldVector3D<T> position, final Frame frame, final FieldAbsoluteDate<T> date) {

  420.         final T zero = position.getAlpha().getField().getZero();
  421.         T result = zero;
  422.         final Map<ExtendedPVCoordinatesProvider, Double> otherOccultingBodies = getOtherOccultingBodies();
  423.         final FieldVector3D<T> sunPosition = sun.getPVCoordinates(date, frame).getPosition();
  424.         final int n = otherOccultingBodies.size() + 1;

  425.         if (n > 1) {

  426.             final List<FieldVector3D<T>> occultingBodyPositions = new ArrayList<FieldVector3D<T>>(n);
  427.             final T[] occultingBodyRadiuses = MathArrays.buildArray(zero.getField(), n);

  428.             // Central body
  429.             occultingBodyPositions.add(0, new FieldVector3D<>(zero, zero, zero));
  430.             occultingBodyRadiuses[0] = zero.add(getEquatorialRadius());

  431.             // Other occulting bodies
  432.             int k = 1;
  433.             for (Map.Entry<ExtendedPVCoordinatesProvider, Double> entry: otherOccultingBodies.entrySet()) {
  434.                 occultingBodyPositions.add(k, entry.getKey().getPVCoordinates(date, frame).getPosition());
  435.                 occultingBodyRadiuses[k] = zero.newInstance(entry.getValue());
  436.                 ++k;
  437.             }
  438.             for (int i = 0; i < n; ++i) {

  439.                 // Lighting ratio computations
  440.                 final T eclipseRatioI = getGeneralEclipseRatio(position,
  441.                                                                occultingBodyPositions.get(i),
  442.                                                                occultingBodyRadiuses[i],
  443.                                                                sunPosition,
  444.                                                                zero.add(Constants.SUN_RADIUS));

  445.                 // First body totally occults Sun, full eclipse is occuring.
  446.                 if (eclipseRatioI.getReal() == 0.0) {
  447.                     return zero;
  448.                 }


  449.                 result = result.add(eclipseRatioI);


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

  452.                     final T eclipseRatioJ = getGeneralEclipseRatio(position,
  453.                                                                    occultingBodyPositions.get(i),
  454.                                                                    occultingBodyRadiuses[j],
  455.                                                                    sunPosition,
  456.                                                                    zero.add(Constants.SUN_RADIUS));
  457.                     final T eclipseRatioIJ = getGeneralEclipseRatio(position,
  458.                                                                     occultingBodyPositions.get(i),
  459.                                                                     occultingBodyRadiuses[i],
  460.                                                                     occultingBodyPositions.get(j),
  461.                                                                     occultingBodyRadiuses[j]);

  462.                     final T alphaJ = getGeneralEclipseAngles(position,
  463.                                                              occultingBodyPositions.get(i),
  464.                                                              occultingBodyRadiuses[i],
  465.                                                              occultingBodyPositions.get(j),
  466.                                                              occultingBodyRadiuses[j])[2];

  467.                     final T alphaSun = getEclipseAngles(sunPosition, position)[2];
  468.                     final T alphaJSq = alphaJ.multiply(alphaJ);
  469.                     final T alphaSunSq = alphaSun.multiply(alphaSun);
  470.                     final T mutualEclipseCorrection = eclipseRatioIJ.negate().add(1).multiply(alphaJSq).divide(alphaSunSq);

  471.                     // Secondary body totaly occults Sun, no more computations are required, full eclipse is occuring.
  472.                     if (eclipseRatioJ.getReal() ==  0.0 ) {
  473.                         return zero;
  474.                     }

  475.                     // Secondary body partially occults Sun
  476.                     else if (eclipseRatioJ.getReal() != 1) {
  477.                         result = result.subtract(mutualEclipseCorrection);
  478.                     }
  479.                 }
  480.             }
  481.             // Final term
  482.             result = result.subtract(n - 1);
  483.         } else {
  484.             // only central body is considered
  485.             result = getLightingRatio(position, frame, date);
  486.         }
  487.         return result;
  488.     }


  489.     /** {@inheritDoc} */
  490.     @Override
  491.     public List<ParameterDriver> getParametersDrivers() {
  492.         return spacecraft.getRadiationParametersDrivers();
  493.     }

  494.     /** Compute min of two values, one double and one field element.
  495.      * @param d double value
  496.      * @param f field element
  497.      * @param <T> type fo the field elements
  498.      * @return min value
  499.      */
  500.     private <T extends CalculusFieldElement<T>> T min(final double d, final T f) {
  501.         return (f.getReal() > d) ? f.getField().getZero().newInstance(d) : f;
  502.     }

  503.     /** Compute max of two values, one double and one field element.
  504.      * @param d double value
  505.      * @param f field element
  506.      * @param <T> type fo the field elements
  507.      * @return max value
  508.      */
  509.     private <T extends CalculusFieldElement<T>> T max(final double d, final T f) {
  510.         return (f.getReal() <= d) ? f.getField().getZero().newInstance(d) : f;
  511.     }

  512. }