SolarRadiationPressure.java

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

  18. import java.util.stream.Stream;

  19. import org.hipparchus.Field;
  20. import org.hipparchus.RealFieldElement;
  21. import org.hipparchus.exception.LocalizedCoreFormats;
  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.errors.OrekitException;
  28. import org.orekit.errors.OrekitIllegalArgumentException;
  29. import org.orekit.errors.OrekitMessages;
  30. import org.orekit.forces.AbstractForceModel;
  31. import org.orekit.frames.Frame;
  32. import org.orekit.propagation.FieldSpacecraftState;
  33. import org.orekit.propagation.SpacecraftState;
  34. import org.orekit.propagation.events.AbstractDetector;
  35. import org.orekit.propagation.events.EventDetector;
  36. import org.orekit.propagation.events.FieldAbstractDetector;
  37. import org.orekit.propagation.events.FieldEventDetector;
  38. import org.orekit.propagation.events.handlers.EventHandler;
  39. import org.orekit.propagation.events.handlers.FieldEventHandler;
  40. import org.orekit.time.AbsoluteDate;
  41. import org.orekit.time.FieldAbsoluteDate;
  42. import org.orekit.utils.Constants;
  43. import org.orekit.utils.ExtendedPVCoordinatesProvider;
  44. import org.orekit.utils.PVCoordinatesProvider;
  45. import org.orekit.utils.ParameterDriver;
  46. import org.orekit.utils.TimeStampedFieldPVCoordinates;
  47. import org.orekit.utils.TimeStampedPVCoordinates;

  48. /** Solar radiation pressure force model.
  49.  *
  50.  * @author Fabien Maussion
  51.  * @author Édouard Delente
  52.  * @author Véronique Pommier-Maurussane
  53.  * @author Pascal Parraud
  54.  */
  55. public class SolarRadiationPressure extends AbstractForceModel {

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

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

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

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

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

  66.     /** Central body model. */
  67.     private final double equatorialRadius;

  68.     /** Spacecraft. */
  69.     private final RadiationSensitive spacecraft;

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

  87.     /** Simple constructor with default reference values.
  88.      * <p>When this constructor is used, the reference values are:</p>
  89.      * <ul>
  90.      *   <li>d<sub>ref</sub> = 149597870000.0 m</li>
  91.      *   <li>p<sub>ref</sub> = 4.56 10<sup>-6</sup> N/m²</li>
  92.      * </ul>
  93.      * @param sun Sun model
  94.      * @param equatorialRadius spherical shape model (for umbra/penumbra computation)
  95.      * @param spacecraft the object physical and geometrical information
  96.      * @since 9.2
  97.      */
  98.     public SolarRadiationPressure(final ExtendedPVCoordinatesProvider sun, final double equatorialRadius,
  99.                                   final RadiationSensitive spacecraft) {
  100.         this(D_REF, P_REF, sun, equatorialRadius, spacecraft);
  101.     }

  102.     /** Complete constructor.
  103.      * <p>Note that reference solar radiation pressure <code>pRef</code> in
  104.      * N/m² is linked to solar flux SF in W/m² using
  105.      * formula pRef = SF/c where c is the speed of light (299792458 m/s). So
  106.      * at 1UA a 1367 W/m² solar flux is a 4.56 10<sup>-6</sup>
  107.      * N/m² solar radiation pressure.</p>
  108.      * @param dRef reference distance for the solar radiation pressure (m)
  109.      * @param pRef reference solar radiation pressure at dRef (N/m²)
  110.      * @param sun Sun model
  111.      * @param equatorialRadius spherical shape model (for umbra/penumbra computation)
  112.      * @param spacecraft the object physical and geometrical information
  113.      * @deprecated as of 9.2 replaced by {@link #SolarRadiationPressure(double, double,
  114.      * ExtendedPVCoordinatesProvider, double, RadiationSensitive)}
  115.      */
  116.     @Deprecated
  117.     public SolarRadiationPressure(final double dRef, final double pRef,
  118.                                   final PVCoordinatesProvider sun,
  119.                                   final double equatorialRadius,
  120.                                   final RadiationSensitive spacecraft) {
  121.         this.kRef = pRef * dRef * dRef;
  122.         if (sun instanceof ExtendedPVCoordinatesProvider) {
  123.             this.sun = (ExtendedPVCoordinatesProvider) sun;
  124.         } else {
  125.             this.sun = new ExtendedPVCoordinatesProvider() {

  126.                 /** {@inheritDoc} */
  127.                 @Override
  128.                 public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate date, final Frame frame)
  129.                     throws OrekitException {
  130.                     // delegate to raw Sun provider
  131.                     return sun.getPVCoordinates(date, frame);
  132.                 }

  133.                 /** {@inheritDoc} */
  134.                 @Override
  135.                 public <T extends RealFieldElement<T>> TimeStampedFieldPVCoordinates<T>
  136.                     getPVCoordinates(final FieldAbsoluteDate<T> date, final Frame frame)
  137.                         throws OrekitException {
  138.                     // SRP was created with a provider that does not support fields,
  139.                     // but the fields methods are called
  140.                     throw new OrekitIllegalArgumentException(LocalizedCoreFormats.UNSUPPORTED_OPERATION);
  141.                 }

  142.             };
  143.         };
  144.         this.equatorialRadius = equatorialRadius;
  145.         this.spacecraft = spacecraft;
  146.     }

  147.     /** Complete constructor.
  148.      * <p>Note that reference solar radiation pressure <code>pRef</code> in
  149.      * N/m² is linked to solar flux SF in W/m² using
  150.      * formula pRef = SF/c where c is the speed of light (299792458 m/s). So
  151.      * at 1UA a 1367 W/m² solar flux is a 4.56 10<sup>-6</sup>
  152.      * N/m² solar radiation pressure.</p>
  153.      * @param dRef reference distance for the solar radiation pressure (m)
  154.      * @param pRef reference solar radiation pressure at dRef (N/m²)
  155.      * @param sun Sun model
  156.      * @param equatorialRadius spherical shape model (for umbra/penumbra computation)
  157.      * @param spacecraft the object physical and geometrical information
  158.      * @since 9.2
  159.      */
  160.     public SolarRadiationPressure(final double dRef, final double pRef,
  161.                                   final ExtendedPVCoordinatesProvider sun,
  162.                                   final double equatorialRadius,
  163.                                   final RadiationSensitive spacecraft) {
  164.         this.kRef = pRef * dRef * dRef;
  165.         this.sun  = sun;
  166.         this.equatorialRadius = equatorialRadius;
  167.         this.spacecraft = spacecraft;
  168.     }

  169.     /** {@inheritDoc} */
  170.     @Override
  171.     public boolean dependsOnPositionOnly() {
  172.         return false;
  173.     }

  174.     /** {@inheritDoc} */
  175.     @Override
  176.     public Vector3D acceleration(final SpacecraftState s, final double[] parameters)
  177.         throws OrekitException {

  178.         final AbsoluteDate date         = s.getDate();
  179.         final Frame        frame        = s.getFrame();
  180.         final Vector3D     position     = s.getPVCoordinates().getPosition();
  181.         final Vector3D     sunSatVector = position.subtract(sun.getPVCoordinates(date, frame).getPosition());
  182.         final double       r2           = sunSatVector.getNormSq();

  183.         // compute flux
  184.         final double   ratio = getLightingRatio(position, frame, date);
  185.         final double   rawP  = ratio  * kRef / r2;
  186.         final Vector3D flux  = new Vector3D(rawP / FastMath.sqrt(r2), sunSatVector);

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

  189.     }

  190.     /** {@inheritDoc} */
  191.     @Override
  192.     public <T extends RealFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
  193.                                                                          final T[] parameters)
  194.         throws OrekitException {

  195.         final FieldAbsoluteDate<T> date         = s.getDate();
  196.         final Frame                frame        = s.getFrame();
  197.         final FieldVector3D<T>     position     = s.getPVCoordinates().getPosition();
  198.         final FieldVector3D<T>     sunSatVector = position.subtract(sun.getPVCoordinates(date.toAbsoluteDate(), frame).getPosition());
  199.         final T                    r2           = sunSatVector.getNormSq();

  200.         // compute flux
  201.         final T                ratio = getLightingRatio(position, frame, date);
  202.         final T                rawP  = ratio.divide(r2).multiply(kRef);
  203.         final FieldVector3D<T> flux  = new FieldVector3D<>(rawP.divide(r2.sqrt()), sunSatVector);

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

  206.     }

  207.     /** Get the lighting ratio ([0-1]).
  208.      * @param position the satellite's position in the selected frame.
  209.      * @param frame in which is defined the position
  210.      * @param date the date
  211.      * @return lighting ratio
  212.      * @exception OrekitException if the trajectory is inside the central body
  213.      * @since 7.1
  214.      */
  215.     public double getLightingRatio(final Vector3D position, final Frame frame, final AbsoluteDate date)
  216.         throws OrekitException {

  217.         final Vector3D sunPosition = sun.getPVCoordinates(date, frame).getPosition();
  218.         if (sunPosition.getNorm() < 2 * Constants.SUN_RADIUS) {
  219.             // we are in fact computing a trajectory around Sun (or solar system barycenter),
  220.             // not around a planet,we consider lighting ratio is always 1
  221.             return 1.0;
  222.         }

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

  225.         // Sat-Sun / Sat-CentralBody angle
  226.         final double sunSatCentralBodyAngle = angle[0];

  227.         // Central Body apparent radius
  228.         final double alphaCentral = angle[1];

  229.         // Sun apparent radius
  230.         final double alphaSun = angle[2];

  231.         double result = 1.0;

  232.         // Is the satellite in complete umbra ?
  233.         if (sunSatCentralBodyAngle - alphaCentral + alphaSun <= ANGULAR_MARGIN) {
  234.             result = 0.0;
  235.         } else if (sunSatCentralBodyAngle - alphaCentral - alphaSun < -ANGULAR_MARGIN) {
  236.             // Compute a lighting ratio in penumbra
  237.             final double sEA2    = sunSatCentralBodyAngle * sunSatCentralBodyAngle;
  238.             final double oo2sEA  = 1.0 / (2. * sunSatCentralBodyAngle);
  239.             final double aS2     = alphaSun * alphaSun;
  240.             final double aE2     = alphaCentral * alphaCentral;
  241.             final double aE2maS2 = aE2 - aS2;

  242.             final double alpha1  = (sEA2 - aE2maS2) * oo2sEA;
  243.             final double alpha2  = (sEA2 + aE2maS2) * oo2sEA;

  244.             // Protection against numerical inaccuracy at boundaries
  245.             final double almost0 = Precision.SAFE_MIN;
  246.             final double almost1 = FastMath.nextDown(1.0);
  247.             final double a1oaS   = FastMath.min(almost1, FastMath.max(-almost1, alpha1 / alphaSun));
  248.             final double aS2ma12 = FastMath.max(almost0, aS2 - alpha1 * alpha1);
  249.             final double a2oaE   = FastMath.min(almost1, FastMath.max(-almost1, alpha2 / alphaCentral));
  250.             final double aE2ma22 = FastMath.max(almost0, aE2 - alpha2 * alpha2);

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

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

  255.         return result;

  256.     }

  257.     /** Get the lighting ratio ([0-1]).
  258.      * @param position the satellite's position in the selected frame.
  259.      * @param frame in which is defined the position
  260.      * @param date the date
  261.      * @param <T> extends RealFieldElement
  262.      * @return lighting ratio
  263.      * @exception OrekitException if the trajectory is inside the central body
  264.      * @since 7.1
  265.      */
  266.     public <T extends RealFieldElement<T>> T getLightingRatio(final FieldVector3D<T> position,
  267.                                                               final Frame frame,
  268.                                                               final FieldAbsoluteDate<T> date)
  269.         throws OrekitException {

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

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

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

  279.         // Sat-Sun / Sat-CentralBody angle
  280.         final T sunsatCentralBodyAngle = angle[0];

  281.         // Central Body apparent radius
  282.         final T alphaCentral = angle[1];

  283.         // Sun apparent radius
  284.         final T alphaSun = angle[2];

  285.         T result = one;
  286.         // Is the satellite in complete umbra ?
  287.         if (sunsatCentralBodyAngle.getReal() - alphaCentral.getReal() + alphaSun.getReal() <= ANGULAR_MARGIN) {
  288.             result = date.getField().getZero();
  289.         } else if (sunsatCentralBodyAngle.getReal() - alphaCentral.getReal() - alphaSun.getReal() < -ANGULAR_MARGIN) {
  290.             // Compute a lighting ratio in penumbra
  291.             final T sEA2    = sunsatCentralBodyAngle.multiply(sunsatCentralBodyAngle);
  292.             final T oo2sEA  = sunsatCentralBodyAngle.multiply(2).reciprocal();
  293.             final T aS2     = alphaSun.multiply(alphaSun);
  294.             final T aE2     = alphaCentral.multiply(alphaCentral);
  295.             final T aE2maS2 = aE2.subtract(aS2);

  296.             final T alpha1  = sEA2.subtract(aE2maS2).multiply(oo2sEA);
  297.             final T alpha2  = sEA2.add(aE2maS2).multiply(oo2sEA);

  298.             // Protection against numerical inaccuracy at boundaries
  299.             final double almost0 = Precision.SAFE_MIN;
  300.             final double almost1 = FastMath.nextDown(1.0);
  301.             final T a1oaS   = min(almost1, max(-almost1, alpha1.divide(alphaSun)));
  302.             final T aS2ma12 = max(almost0, aS2.subtract(alpha1.multiply(alpha1)));
  303.             final T a2oaE   = min(almost1, max(-almost1, alpha2.divide(alphaCentral)));
  304.             final T aE2ma22 = max(almost0, aE2.subtract(alpha2.multiply(alpha2)));

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

  307.             result = one.subtract(P1.add(P2).divide(aS2.multiply(FastMath.PI)));
  308.         }

  309.         return result;
  310.     }

  311.     /** {@inheritDoc} */
  312.     @Override
  313.     public Stream<EventDetector> getEventsDetectors() {
  314.         return Stream.of(new UmbraDetector(), new PenumbraDetector());
  315.     }

  316.     /** {@inheritDoc} */
  317.     @Override
  318.     public <T extends RealFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventsDetectors(final Field<T> field) {
  319.         return Stream.of(new FieldUmbraDetector<>(field), new FieldPenumbraDetector<>(field));
  320.     }

  321.     /** {@inheritDoc} */
  322.     @Override
  323.     public ParameterDriver[] getParametersDrivers() {
  324.         return spacecraft.getRadiationParametersDrivers();
  325.     }

  326.     /** Get the useful angles for eclipse computation.
  327.      * @param sunPosition Sun position in the selected frame
  328.      * @param position the satellite's position in the selected frame
  329.      * @return the 3 angles {(satCentral, satSun), Central body apparent radius, Sun apparent radius}
  330.      * @exception OrekitException if the trajectory is inside the central body
  331.      */
  332.     private double[] getEclipseAngles(final Vector3D sunPosition, final Vector3D position)
  333.         throws OrekitException {
  334.         final double[] angle = new double[3];

  335.         final Vector3D satSunVector = sunPosition.subtract(position);

  336.         // Sat-Sun / Sat-CentralBody angle
  337.         angle[0] = Vector3D.angle(satSunVector, position.negate());

  338.         // Central body apparent radius
  339.         final double r = position.getNorm();
  340.         if (r <= equatorialRadius) {
  341.             throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE, r);
  342.         }
  343.         angle[1] = FastMath.asin(equatorialRadius / r);

  344.         // Sun apparent radius
  345.         angle[2] = FastMath.asin(Constants.SUN_RADIUS / satSunVector.getNorm());

  346.         return angle;
  347.     }

  348.     /** Get the useful angles for eclipse computation.
  349.      * @param sunPosition Sun position in the selected frame
  350.      * @param position the satellite's position in the selected frame.
  351.      * @param <T> extends RealFieldElement
  352.      * @return the 3 angles {(satCentral, satSun), Central body apparent radius, Sun apparent radius}
  353.      * @exception OrekitException if the trajectory is inside the central body
  354.      */
  355.     private <T extends RealFieldElement<T>> T[] getEclipseAngles(final FieldVector3D<T> sunPosition, final FieldVector3D<T> position)
  356.         throws OrekitException {
  357.         final T[] angle = MathArrays.buildArray(position.getX().getField(), 3);

  358.         final FieldVector3D<T> mP           = position.negate();
  359.         final FieldVector3D<T> satSunVector = mP.add(sunPosition);

  360.         // Sat-Sun / Sat-CentralBody angle
  361.         angle[0] = FieldVector3D.angle(satSunVector, mP);

  362.         // Central body apparent radius
  363.         final T r = position.getNorm();
  364.         if (r.getReal() <= equatorialRadius) {
  365.             throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE, r);
  366.         }
  367.         angle[1] = r.reciprocal().multiply(equatorialRadius).asin();

  368.         // Sun apparent radius
  369.         angle[2] = satSunVector.getNorm().reciprocal().multiply(Constants.SUN_RADIUS).asin();

  370.         return angle;
  371.     }

  372.     /** Compute min of two values, one double and one field element.
  373.      * @param d double value
  374.      * @param f field element
  375.      * @param <T> type fo the field elements
  376.      * @return min value
  377.      */
  378.     private <T extends RealFieldElement<T>> T min(final double d, final T f) {
  379.         return (f.getReal() > d) ? f.getField().getZero().add(d) : f;
  380.     }

  381.     /** Compute max of two values, one double and one field element.
  382.      * @param d double value
  383.      * @param f field element
  384.      * @param <T> type fo the field elements
  385.      * @return max value
  386.      */
  387.     private <T extends RealFieldElement<T>> T max(final double d, final T f) {
  388.         return (f.getReal() <= d) ? f.getField().getZero().add(d) : f;
  389.     }

  390.     /** This class defines the umbra entry/exit detector. */
  391.     private class UmbraDetector extends AbstractDetector<UmbraDetector> {

  392.         /** Serializable UID. */
  393.         private static final long serialVersionUID = 20141228L;

  394.         /** Build a new instance. */
  395.         UmbraDetector() {
  396.             super(60.0, 1.0e-3, DEFAULT_MAX_ITER, new EventHandler<UmbraDetector>() {

  397.                 /** {@inheritDoc} */
  398.                 public Action eventOccurred(final SpacecraftState s, final UmbraDetector detector,
  399.                                             final boolean increasing) {
  400.                     return Action.RESET_DERIVATIVES;
  401.                 }

  402.             });
  403.         }

  404.         /** Private constructor with full parameters.
  405.          * <p>
  406.          * This constructor is private as users are expected to use the builder
  407.          * API with the various {@code withXxx()} methods to set up the instance
  408.          * in a readable manner without using a huge amount of parameters.
  409.          * </p>
  410.          * @param maxCheck maximum checking interval (s)
  411.          * @param threshold convergence threshold (s)
  412.          * @param maxIter maximum number of iterations in the event time search
  413.          * @param handler event handler to call at event occurrences
  414.          * @since 6.1
  415.          */
  416.         private UmbraDetector(final double maxCheck, final double threshold,
  417.                               final int maxIter, final EventHandler<? super UmbraDetector> handler) {
  418.             super(maxCheck, threshold, maxIter, handler);
  419.         }

  420.         /** {@inheritDoc} */
  421.         @Override
  422.         protected UmbraDetector create(final double newMaxCheck, final double newThreshold,
  423.                                        final int newMaxIter, final EventHandler<? super UmbraDetector> newHandler) {
  424.             return new UmbraDetector(newMaxCheck, newThreshold, newMaxIter, newHandler);
  425.         }

  426.         /** The G-function is the difference between the Sun-Sat-Central-Body angle and
  427.          * the central body apparent radius.
  428.          * @param s the current state information : date, kinematics, attitude
  429.          * @return value of the g function
  430.          * @exception OrekitException if sun or spacecraft position cannot be computed
  431.          */
  432.         public double g(final SpacecraftState s) throws OrekitException {
  433.             final double[] angle = getEclipseAngles(sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
  434.                                                     s.getPVCoordinates().getPosition());
  435.             return angle[0] - angle[1] + angle[2] - ANGULAR_MARGIN;
  436.         }

  437.     }

  438.     /** This class defines the penumbra entry/exit detector. */
  439.     private class PenumbraDetector extends AbstractDetector<PenumbraDetector> {

  440.         /** Serializable UID. */
  441.         private static final long serialVersionUID = 20141228L;

  442.         /** Build a new instance. */
  443.         PenumbraDetector() {
  444.             super(60.0, 1.0e-3, DEFAULT_MAX_ITER, new EventHandler<PenumbraDetector>() {

  445.                 /** {@inheritDoc} */
  446.                 public Action eventOccurred(final SpacecraftState s, final PenumbraDetector detector,
  447.                                             final boolean increasing) {
  448.                     return Action.RESET_DERIVATIVES;
  449.                 }

  450.             });
  451.         }

  452.         /** Private constructor with full parameters.
  453.          * <p>
  454.          * This constructor is private as users are expected to use the builder
  455.          * API with the various {@code withXxx()} methods to set up the instance
  456.          * in a readable manner without using a huge amount of parameters.
  457.          * </p>
  458.          * @param maxCheck maximum checking interval (s)
  459.          * @param threshold convergence threshold (s)
  460.          * @param maxIter maximum number of iterations in the event time search
  461.          * @param handler event handler to call at event occurrences
  462.          * @since 6.1
  463.          */
  464.         private PenumbraDetector(final double maxCheck, final double threshold,
  465.                                  final int maxIter, final EventHandler<? super PenumbraDetector> handler) {
  466.             super(maxCheck, threshold, maxIter, handler);
  467.         }

  468.         /** {@inheritDoc} */
  469.         @Override
  470.         protected PenumbraDetector create(final double newMaxCheck, final double newThreshold,
  471.                                           final int newMaxIter, final EventHandler<? super PenumbraDetector> newHandler) {
  472.             return new PenumbraDetector(newMaxCheck, newThreshold, newMaxIter, newHandler);
  473.         }

  474.         /** The G-function is the difference between the Sun-Sat-Central-Body angle and
  475.          * the sum of the central body and Sun's apparent radius.
  476.          * @param s the current state information : date, kinematics, attitude
  477.          * @return value of the g function
  478.          * @exception OrekitException if sun or spacecraft position cannot be computed
  479.          */
  480.         public double g(final SpacecraftState s) throws OrekitException {
  481.             final double[] angle = getEclipseAngles(sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
  482.                                                     s.getPVCoordinates().getPosition());
  483.             return angle[0] - angle[1] - angle[2] + ANGULAR_MARGIN;
  484.         }

  485.     }

  486.     /** This class defines the umbra entry/exit detector.
  487.      * @since 9.2
  488.      */
  489.     private class FieldUmbraDetector<T extends RealFieldElement<T>>
  490.         extends FieldAbstractDetector<FieldUmbraDetector<T>, T> {

  491.         /** Build a new instance.
  492.          * @param field field to which elements belong
  493.          */
  494.         FieldUmbraDetector(final Field<T> field) {
  495.             super(field.getZero().add(60.0), field.getZero().add(1.0e-3),
  496.                   DEFAULT_MAX_ITER, new FieldEventHandler<FieldUmbraDetector<T>, T>() {

  497.                       /** {@inheritDoc} */
  498.                       public Action eventOccurred(final FieldSpacecraftState<T> s,
  499.                                                   final FieldUmbraDetector<T> detector,
  500.                                                   final boolean increasing) {
  501.                           return Action.RESET_DERIVATIVES;
  502.                       }

  503.                   });
  504.         }

  505.         /** Private constructor with full parameters.
  506.          * <p>
  507.          * This constructor is private as users are expected to use the builder
  508.          * API with the various {@code withXxx()} methods to set up the instance
  509.          * in a readable manner without using a huge amount of parameters.
  510.          * </p>
  511.          * @param maxCheck maximum checking interval (s)
  512.          * @param threshold convergence threshold (s)
  513.          * @param maxIter maximum number of iterations in the event time search
  514.          * @param handler event handler to call at event occurrences
  515.          */
  516.         private FieldUmbraDetector(final T maxCheck, final T threshold,
  517.                                    final int maxIter,
  518.                                    final FieldEventHandler<? super FieldUmbraDetector<T>, T> handler) {
  519.             super(maxCheck, threshold, maxIter, handler);
  520.         }

  521.         /** {@inheritDoc} */
  522.         @Override
  523.         protected FieldUmbraDetector<T> create(final T newMaxCheck, final T newThreshold,
  524.                                                final int newMaxIter,
  525.                                                final FieldEventHandler<? super FieldUmbraDetector<T>, T> newHandler) {
  526.             return new FieldUmbraDetector<>(newMaxCheck, newThreshold, newMaxIter, newHandler);
  527.         }

  528.         /** The G-function is the difference between the Sun-Sat-Central-Body angle and
  529.          * the central body apparent radius.
  530.          * @param s the current state information : date, kinematics, attitude
  531.          * @return value of the g function
  532.          * @exception OrekitException if sun or spacecraft position cannot be computed
  533.          */
  534.         public T g(final FieldSpacecraftState<T> s) throws OrekitException {
  535.             final T[] angle = getEclipseAngles(sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
  536.                                                s.getPVCoordinates().getPosition());
  537.             return angle[0].subtract(angle[1]).add(angle[2]).subtract(ANGULAR_MARGIN);
  538.         }

  539.     }

  540.     /** This class defines the penumbra entry/exit detector.
  541.      * @since 9.2
  542.      */
  543.     private class FieldPenumbraDetector<T extends RealFieldElement<T>>
  544.           extends FieldAbstractDetector<FieldPenumbraDetector<T>, T> {

  545.         /** Build a new instance.
  546.          * @param field field to which elements belong
  547.          */
  548.         FieldPenumbraDetector(final Field<T> field) {
  549.             super(field.getZero().add(60.0), field.getZero().add(1.0e-3),
  550.                   DEFAULT_MAX_ITER, new FieldEventHandler<FieldPenumbraDetector<T>, T>() {

  551.                       /** {@inheritDoc} */
  552.                       public Action eventOccurred(final FieldSpacecraftState<T> s,
  553.                                                   final FieldPenumbraDetector<T> detector,
  554.                                                   final boolean increasing) {
  555.                           return Action.RESET_DERIVATIVES;
  556.                       }

  557.                   });
  558.         }

  559.         /** Private constructor with full parameters.
  560.          * <p>
  561.          * This constructor is private as users are expected to use the builder
  562.          * API with the various {@code withXxx()} methods to set up the instance
  563.          * in a readable manner without using a huge amount of parameters.
  564.          * </p>
  565.          * @param maxCheck maximum checking interval (s)
  566.          * @param threshold convergence threshold (s)
  567.          * @param maxIter maximum number of iterations in the event time search
  568.          * @param handler event handler to call at event occurrences
  569.          */
  570.         private FieldPenumbraDetector(final T maxCheck, final T threshold,
  571.                                       final int maxIter,
  572.                                       final FieldEventHandler<? super FieldPenumbraDetector<T>, T> handler) {
  573.             super(maxCheck, threshold, maxIter, handler);
  574.         }

  575.         /** {@inheritDoc} */
  576.         @Override
  577.         protected FieldPenumbraDetector<T> create(final T newMaxCheck, final T newThreshold,
  578.                                                   final int newMaxIter,
  579.                                                   final FieldEventHandler<? super FieldPenumbraDetector<T>, T> newHandler) {
  580.             return new FieldPenumbraDetector<>(newMaxCheck, newThreshold, newMaxIter, newHandler);
  581.         }

  582.         /** The G-function is the difference between the Sun-Sat-Central-Body angle and
  583.          * the sum of the central body and Sun's apparent radius.
  584.          * @param s the current state information : date, kinematics, attitude
  585.          * @return value of the g function
  586.          * @exception OrekitException if sun or spacecraft position cannot be computed
  587.          */
  588.         public T g(final FieldSpacecraftState<T> s) throws OrekitException {
  589.             final T[] angle = getEclipseAngles(sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
  590.                                                s.getPVCoordinates().getPosition());
  591.             return angle[0].subtract(angle[1]).subtract(angle[2]).add(ANGULAR_MARGIN);
  592.         }

  593.     }

  594. }