diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 8d99e023887b6cb138aa78cb06b880624a32e9c8..de8b9bdc77a2bc000f4b5162b4a9824a3ceac24a 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -21,6 +21,9 @@
</properties>
<body>
<release version="11.0" date="TBD" description="TBD">
+ <action dev="thomas" type="fix" issue="702">
+ Added possibility to take in account several bodies while computing SRP perturbation.
+ </action>
<action dev="bryan" type="update" issue="793">
Updated SP3File visibility to public.
</action>
diff --git a/src/main/java/org/orekit/forces/radiation/AbstractRadiationForceModel.java b/src/main/java/org/orekit/forces/radiation/AbstractRadiationForceModel.java
index 52775f4d47def0aaa0672e23c02fbadf314a88c9..aab9a7b1efad77423cca300224695ebeaf2898f5 100644
--- a/src/main/java/org/orekit/forces/radiation/AbstractRadiationForceModel.java
+++ b/src/main/java/org/orekit/forces/radiation/AbstractRadiationForceModel.java
@@ -16,6 +16,9 @@
*/
package org.orekit.forces.radiation;
+import java.lang.reflect.Array;
+import java.util.HashMap;
+import java.util.Map;
import java.util.stream.Stream;
import org.hipparchus.Field;
@@ -56,14 +59,19 @@ public abstract class AbstractRadiationForceModel extends AbstractForceModel {
/** Sun model. */
private final ExtendedPVCoordinatesProvider sun;
+ /** Other occulting bodies to consider. The Moon for instance. */
+ private final Map<ExtendedPVCoordinatesProvider, Double> otherOccultingBodies;
+
/**
- * Constructor.
+ * Default constructor.
+ * Only central body is considered.
* @param sun Sun model
- * @param equatorialRadius spherical shape model (for umbra/penumbra computation)
+ * @param equatorialRadius central body spherical shape model (for umbra/penumbra computation)
*/
protected AbstractRadiationForceModel(final ExtendedPVCoordinatesProvider sun, final double equatorialRadius) {
- this.sun = sun;
- this.equatorialRadius = equatorialRadius;
+ this.sun = sun;
+ this.equatorialRadius = equatorialRadius;
+ this.otherOccultingBodies = new HashMap<>();
}
/** {@inheritDoc} */
@@ -75,13 +83,34 @@ public abstract class AbstractRadiationForceModel extends AbstractForceModel {
/** {@inheritDoc} */
@Override
public Stream<EventDetector> getEventsDetectors() {
- return Stream.of(new UmbraDetector(), new PenumbraDetector());
+ final EventDetector[] detectors = new EventDetector[2 + 2 * otherOccultingBodies.size()];
+ detectors[0] = new UmbraDetector();
+ detectors[1] = new PenumbraDetector();
+ int i = 2;
+ for (Map.Entry<ExtendedPVCoordinatesProvider, Double> entry : otherOccultingBodies.entrySet()) {
+ detectors[i] = new GeneralUmbraDetector(entry.getKey(), entry.getValue());
+ detectors[i + 1] = new GeneralPenumbraDetector(entry.getKey(), entry.getValue());
+ i = i + 2;
+ }
+ return Stream.of(detectors);
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventsDetectors(final Field<T> field) {
- return Stream.of(new FieldUmbraDetector<>(field), new FieldPenumbraDetector<>(field));
+ final T zero = field.getZero();
+ @SuppressWarnings("unchecked")
+ final FieldEventDetector<T>[] detectors = (FieldEventDetector<T>[]) Array.newInstance(FieldEventDetector.class,
+ 2 + 2 * otherOccultingBodies.size());
+ detectors[0] = new FieldUmbraDetector<>(field);
+ detectors[1] = new FieldPenumbraDetector<>(field);
+ int i = 2;
+ for (Map.Entry<ExtendedPVCoordinatesProvider, Double> entry : otherOccultingBodies.entrySet()) {
+ detectors[i] = new FieldGeneralUmbraDetector<>(field, entry.getKey(), zero.newInstance(entry.getValue()));
+ detectors[i + 1] = new FieldGeneralPenumbraDetector<>(field, entry.getKey(), zero.newInstance(entry.getValue()));
+ i = i + 2;
+ }
+ return Stream.of(detectors);
}
/**
@@ -111,6 +140,35 @@ public abstract class AbstractRadiationForceModel extends AbstractForceModel {
return angle;
}
+
+ /**
+ * Get the useful angles for eclipse computation.
+ * @param position the satellite's position in the selected frame
+ * @param occultingPosition Oculting body position in the selected frame
+ * @param occultingRadius Occulting body mean radius
+ * @param occultedPosition Occulted body position in the selected frame
+ * @param occultedRadius Occulted body mean radius
+ * @return the 3 angles {(satOcculting, satOcculted), Occulting body apparent radius, Occulted body apparent radius}
+ */
+ protected double[] getGeneralEclipseAngles(final Vector3D position, final Vector3D occultingPosition, final double occultingRadius,
+ final Vector3D occultedPosition, final double occultedRadius) {
+ final double[] angle = new double[3];
+
+ final Vector3D satOccultedVector = occultedPosition.subtract(position);
+ final Vector3D satOccultingVector = occultingPosition.subtract(position);
+
+ // Sat-Occulted / Sat-Occulting angle
+ angle[0] = Vector3D.angle(satOccultedVector, satOccultingVector);
+
+ // Occulting body apparent radius
+ angle[1] = FastMath.asin(occultingRadius / satOccultingVector.getNorm());
+
+ // Occulted body apparent radius
+ angle[2] = FastMath.asin(occultedRadius / satOccultedVector.getNorm());
+
+ return angle;
+ }
+
/**
* Get the useful angles for eclipse computation.
* @param sunPosition Sun position in the selected frame
@@ -140,6 +198,62 @@ public abstract class AbstractRadiationForceModel extends AbstractForceModel {
return angle;
}
+ /**
+ * Get the useful angles for eclipse computation.
+ * @param occultingPosition Oculting body position in the selected frame
+ * @param occultingRadius Occulting body mean radius
+ * @param occultedPosition Occulted body position in the selected frame
+ * @param occultedRadius Occulted body mean radius
+ * @param position the satellite's position in the selected frame
+ * @param <T> extends RealFieldElement
+ * @return the 3 angles {(satOcculting, satOcculted), Occulting body apparent radius, Occulted body apparent radius}
+ */
+ protected <T extends CalculusFieldElement<T>> T[] getGeneralEclipseAngles(final FieldVector3D<T> position,
+ final FieldVector3D<T> occultingPosition, final T occultingRadius,
+ final FieldVector3D<T> occultedPosition, final T occultedRadius) {
+ final T[] angle = MathArrays.buildArray(position.getX().getField(), 3);
+
+ final FieldVector3D<T> satOccultedVector = occultedPosition.subtract(position);
+ final FieldVector3D<T> satOccultingVector = occultingPosition.subtract(position);
+
+ // Sat-Occulted / Sat-Occulting angle
+ angle[0] = FieldVector3D.angle(satOccultedVector, satOccultingVector);
+
+ // Occulting body apparent radius
+ angle[1] = satOccultingVector.getNorm().reciprocal().multiply(occultingRadius).asin();
+
+ // Occulted body apparent radius
+ angle[2] = satOccultedVector.getNorm().reciprocal().multiply(occultedRadius).asin();
+
+ return angle;
+ }
+
+ /**
+ * Add a new occulting body.
+ * Central body is already considered, it shall not be added this way.
+ * @param provider body PV provider
+ * @param radius body mean radius
+ */
+ public void addOccultingBody(final ExtendedPVCoordinatesProvider provider, final double radius) {
+ otherOccultingBodies.put(provider, radius);
+ }
+
+ /**
+ * Getter for other occulting bodies to consider.
+ * @return the map of other occulting bodies and corresponding mean radiuses
+ */
+ public Map<ExtendedPVCoordinatesProvider, Double> getOtherOccultingBodies() {
+ return otherOccultingBodies;
+ }
+
+ /**
+ * Getter for equatorial radius.
+ * @return central body equatorial radius
+ */
+ public double getEquatorialRadius() {
+ return equatorialRadius;
+ }
+
/** This class defines the umbra entry/exit detector. */
private class UmbraDetector extends AbstractDetector<UmbraDetector> {
@@ -169,15 +283,15 @@ public abstract class AbstractRadiationForceModel extends AbstractForceModel {
* @param handler event handler to call at event occurrences
* @since 6.1
*/
- private UmbraDetector(final double maxCheck, final double threshold,
- final int maxIter, final EventHandler<? super UmbraDetector> handler) {
+ private UmbraDetector(final double maxCheck, final double threshold, final int maxIter,
+ final EventHandler<? super UmbraDetector> handler) {
super(maxCheck, threshold, maxIter, handler);
}
/** {@inheritDoc} */
@Override
- protected UmbraDetector create(final double newMaxCheck, final double newThreshold,
- final int newMaxIter, final EventHandler<? super UmbraDetector> newHandler) {
+ protected UmbraDetector create(final double newMaxCheck, final double newThreshold, final int newMaxIter,
+ final EventHandler<? super UmbraDetector> newHandler) {
return new UmbraDetector(newMaxCheck, newThreshold, newMaxIter, newHandler);
}
@@ -222,15 +336,15 @@ public abstract class AbstractRadiationForceModel extends AbstractForceModel {
* @param handler event handler to call at event occurrences
* @since 6.1
*/
- private PenumbraDetector(final double maxCheck, final double threshold,
- final int maxIter, final EventHandler<? super PenumbraDetector> handler) {
+ private PenumbraDetector(final double maxCheck, final double threshold, final int maxIter,
+ final EventHandler<? super PenumbraDetector> handler) {
super(maxCheck, threshold, maxIter, handler);
}
/** {@inheritDoc} */
@Override
- protected PenumbraDetector create(final double newMaxCheck, final double newThreshold,
- final int newMaxIter, final EventHandler<? super PenumbraDetector> newHandler) {
+ protected PenumbraDetector create(final double newMaxCheck, final double newThreshold, final int newMaxIter,
+ final EventHandler<? super PenumbraDetector> newHandler) {
return new PenumbraDetector(newMaxCheck, newThreshold, newMaxIter, newHandler);
}
@@ -279,16 +393,14 @@ public abstract class AbstractRadiationForceModel extends AbstractForceModel {
* @param maxIter maximum number of iterations in the event time search
* @param handler event handler to call at event occurrences
*/
- private FieldUmbraDetector(final T maxCheck, final T threshold,
- final int maxIter,
+ private FieldUmbraDetector(final T maxCheck, final T threshold, final int maxIter,
final FieldEventHandler<? super FieldUmbraDetector<T>, T> handler) {
super(maxCheck, threshold, maxIter, handler);
}
/** {@inheritDoc} */
@Override
- protected FieldUmbraDetector<T> create(final T newMaxCheck, final T newThreshold,
- final int newMaxIter,
+ protected FieldUmbraDetector<T> create(final T newMaxCheck, final T newThreshold, final int newMaxIter,
final FieldEventHandler<? super FieldUmbraDetector<T>, T> newHandler) {
return new FieldUmbraDetector<>(newMaxCheck, newThreshold, newMaxIter, newHandler);
}
@@ -338,16 +450,14 @@ public abstract class AbstractRadiationForceModel extends AbstractForceModel {
* @param maxIter maximum number of iterations in the event time search
* @param handler event handler to call at event occurrences
*/
- private FieldPenumbraDetector(final T maxCheck, final T threshold,
- final int maxIter,
+ private FieldPenumbraDetector(final T maxCheck, final T threshold, final int maxIter,
final FieldEventHandler<? super FieldPenumbraDetector<T>, T> handler) {
super(maxCheck, threshold, maxIter, handler);
}
/** {@inheritDoc} */
@Override
- protected FieldPenumbraDetector<T> create(final T newMaxCheck, final T newThreshold,
- final int newMaxIter,
+ protected FieldPenumbraDetector<T> create(final T newMaxCheck, final T newThreshold, final int newMaxIter,
final FieldEventHandler<? super FieldPenumbraDetector<T>, T> newHandler) {
return new FieldPenumbraDetector<>(newMaxCheck, newThreshold, newMaxIter, newHandler);
}
@@ -365,4 +475,274 @@ public abstract class AbstractRadiationForceModel extends AbstractForceModel {
}
+ /** This class defines the umbra entry/exit detector. */
+ private class GeneralUmbraDetector extends AbstractDetector<GeneralUmbraDetector> {
+
+ /** Occulting body PV provider. */
+ private ExtendedPVCoordinatesProvider provider;
+
+ /** Occulting body mean radius. */
+ private double radius;
+
+ /** Build a new instance.
+ * @param provider occulting body PV provider
+ * @param radius occulting body mean radius
+ */
+ GeneralUmbraDetector(final ExtendedPVCoordinatesProvider provider, final double radius) {
+ super(60.0, 1.0e-3, DEFAULT_MAX_ITER, new EventHandler<GeneralUmbraDetector>() {
+
+ /** {@inheritDoc} */
+ public Action eventOccurred(final SpacecraftState s, final GeneralUmbraDetector detector,
+ final boolean increasing) {
+ return Action.RESET_DERIVATIVES;
+ }
+
+ });
+ this.provider = provider;
+ this.radius = radius;
+ }
+
+ /** Private constructor with full parameters.
+ * <p>
+ * This constructor is private as users are expected to use the builder
+ * API with the various {@code withXxx()} methods to set up the instance
+ * in a readable manner without using a huge amount of parameters.
+ * </p>
+ * @param maxCheck maximum checking interval (s)
+ * @param threshold convergence threshold (s)
+ * @param maxIter maximum number of iterations in the event time search
+ * @param handler event handler to call at event occurrences
+ * @since 6.1
+ */
+ private GeneralUmbraDetector(final double maxCheck, final double threshold, final int maxIter,
+ final EventHandler<? super GeneralUmbraDetector> handler) {
+ super(maxCheck, threshold, maxIter, handler);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected GeneralUmbraDetector create(final double newMaxCheck, final double newThreshold, final int newMaxIter,
+ final EventHandler<? super GeneralUmbraDetector> newHandler) {
+ return new GeneralUmbraDetector(newMaxCheck, newThreshold, newMaxIter, newHandler);
+ }
+
+ /** The G-function is the difference between the Sun-Sat-Central-Body angle and
+ * the central body apparent radius.
+ * @param s the current state information : date, kinematics, attitude
+ * @return value of the g function
+ */
+ public double g(final SpacecraftState s) {
+ final double[] angle = getGeneralEclipseAngles(s.getPVCoordinates().getPosition(),
+ provider.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
+ radius, sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
+ Constants.SUN_RADIUS);
+ return angle[0] - angle[1] + angle[2] - ANGULAR_MARGIN;
+ }
+
+ }
+
+ /** This class defines the umbra entry/exit detector. */
+ private class GeneralPenumbraDetector extends AbstractDetector<GeneralPenumbraDetector> {
+
+ /** Occulting body PV provider. */
+ private ExtendedPVCoordinatesProvider provider;
+
+ /** Occulting body mean radius. */
+ private double radius;
+
+ /** Build a new instance.
+ * @param provider occulting body PV provider
+ * @param radius occulting body mean radius
+ */
+ GeneralPenumbraDetector(final ExtendedPVCoordinatesProvider provider, final double radius) {
+ super(60.0, 1.0e-3, DEFAULT_MAX_ITER, new EventHandler<GeneralPenumbraDetector>() {
+
+ /** {@inheritDoc} */
+ public Action eventOccurred(final SpacecraftState s, final GeneralPenumbraDetector detector,
+ final boolean increasing) {
+ return Action.RESET_DERIVATIVES;
+ }
+
+ });
+ this.provider = provider;
+ this.radius = radius;
+ }
+
+ /** Private constructor with full parameters.
+ * <p>
+ * This constructor is private as users are expected to use the builder
+ * API with the various {@code withXxx()} methods to set up the instance
+ * in a readable manner without using a huge amount of parameters.
+ * </p>
+ * @param maxCheck maximum checking interval (s)
+ * @param threshold convergence threshold (s)
+ * @param maxIter maximum number of iterations in the event time search
+ * @param handler event handler to call at event occurrences
+ * @since 6.1
+ */
+ private GeneralPenumbraDetector(final double maxCheck, final double threshold, final int maxIter,
+ final EventHandler<? super GeneralPenumbraDetector> handler) {
+ super(maxCheck, threshold, maxIter, handler);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected GeneralPenumbraDetector create(final double newMaxCheck, final double newThreshold, final int newMaxIter,
+ final EventHandler<? super GeneralPenumbraDetector> newHandler) {
+ return new GeneralPenumbraDetector(newMaxCheck, newThreshold, newMaxIter, newHandler);
+ }
+
+ /** The G-function is the difference between the Sun-Sat-Central-Body angle and
+ * the central body apparent radius.
+ * @param s the current state information : date, kinematics, attitude
+ * @return value of the g function
+ */
+ public double g(final SpacecraftState s) {
+ final double[] angle = getGeneralEclipseAngles(s.getPVCoordinates().getPosition(),
+ provider.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
+ radius, sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
+ Constants.SUN_RADIUS);
+ return angle[0] - angle[1] - angle[2] + ANGULAR_MARGIN;
+ }
+
+ }
+
+ /** This class defines the umbra entry/exit detector. */
+ private class FieldGeneralUmbraDetector<T extends CalculusFieldElement<T>>
+ extends FieldAbstractDetector<FieldGeneralUmbraDetector<T>, T> {
+
+ /** Occulting body PV provider. */
+ private ExtendedPVCoordinatesProvider provider;
+
+ /** Occulting body mean radius. */
+ private T radius;
+
+ /** Build a new instance.
+ * @param field field to which elements belong
+ * @param provider occulting body PV provider
+ * @param radius occulting body mean radius
+ */
+ FieldGeneralUmbraDetector(final Field<T> field, final ExtendedPVCoordinatesProvider provider, final T radius) {
+ super(field.getZero().add(60.0), field.getZero().add(1.0e-3),
+ DEFAULT_MAX_ITER, new FieldEventHandler<FieldGeneralUmbraDetector<T>, T>() {
+
+ /** {@inheritDoc} */
+ public Action eventOccurred(final FieldSpacecraftState<T> s,
+ final FieldGeneralUmbraDetector<T> detector,
+ final boolean increasing) {
+ return Action.RESET_DERIVATIVES;
+ }
+
+ });
+ this.provider = provider;
+ this.radius = radius;
+ }
+
+ /** Private constructor with full parameters.
+ * <p>
+ * This constructor is private as users are expected to use the builder
+ * API with the various {@code withXxx()} methods to set up the instance
+ * in a readable manner without using a huge amount of parameters.
+ * </p>
+ * @param maxCheck maximum checking interval (s)
+ * @param threshold convergence threshold (s)
+ * @param maxIter maximum number of iterations in the event time search
+ * @param handler event handler to call at event occurrences
+ */
+ private FieldGeneralUmbraDetector(final T maxCheck, final T threshold,
+ final int maxIter,
+ final FieldEventHandler<? super FieldGeneralUmbraDetector<T>, T> handler) {
+ super(maxCheck, threshold, maxIter, handler);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected FieldGeneralUmbraDetector<T> create(final T newMaxCheck, final T newThreshold, final int newMaxIter,
+ final FieldEventHandler<? super FieldGeneralUmbraDetector<T>, T> newHandler) {
+ return new FieldGeneralUmbraDetector<>(newMaxCheck, newThreshold, newMaxIter, newHandler);
+ }
+
+ /** The G-function is the difference between the Sun-Sat-Central-Body angle and
+ * the central body apparent radius.
+ * @param s the current state information : date, kinematics, attitude
+ * @return value of the g function
+ */
+ public T g(final FieldSpacecraftState<T> s) {
+ final T[] angle = getGeneralEclipseAngles(s.getPVCoordinates().getPosition(),
+ provider.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
+ radius, sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
+ s.getA().getField().getZero().add(Constants.SUN_RADIUS));
+ return angle[0].subtract(angle[1]).add(angle[2]).subtract(ANGULAR_MARGIN);
+ }
+
+ }
+
+ /** This class defines the umbra entry/exit detector. */
+ private class FieldGeneralPenumbraDetector<T extends CalculusFieldElement<T>>
+ extends FieldAbstractDetector<FieldGeneralPenumbraDetector<T>, T> {
+
+ /** Occulting body PV provider. */
+ private ExtendedPVCoordinatesProvider provider;
+
+ /** Occulting body mean radius. */
+ private T radius;
+
+ /** Build a new instance.
+ * @param field field to which elements belong
+ * @param provider occulting body PV provider
+ * @param radius occulting body mean radius
+ */
+ FieldGeneralPenumbraDetector(final Field<T> field, final ExtendedPVCoordinatesProvider provider, final T radius) {
+ super(field.getZero().add(60.0), field.getZero().add(1.0e-3),
+ DEFAULT_MAX_ITER, new FieldEventHandler<FieldGeneralPenumbraDetector<T>, T>() {
+
+ /** {@inheritDoc} */
+ public Action eventOccurred(final FieldSpacecraftState<T> s,
+ final FieldGeneralPenumbraDetector<T> detector,
+ final boolean increasing) {
+ return Action.RESET_DERIVATIVES;
+ }
+
+ });
+ this.provider = provider;
+ this.radius = radius;
+ }
+
+ /** Private constructor with full parameters.
+ * <p>
+ * This constructor is private as users are expected to use the builder
+ * API with the various {@code withXxx()} methods to set up the instance
+ * in a readable manner without using a huge amount of parameters.
+ * </p>
+ * @param maxCheck maximum checking interval (s)
+ * @param threshold convergence threshold (s)
+ * @param maxIter maximum number of iterations in the event time search
+ * @param handler event handler to call at event occurrences
+ */
+ private FieldGeneralPenumbraDetector(final T maxCheck, final T threshold, final int maxIter,
+ final FieldEventHandler<? super FieldGeneralPenumbraDetector<T>, T> handler) {
+ super(maxCheck, threshold, maxIter, handler);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected FieldGeneralPenumbraDetector<T> create(final T newMaxCheck, final T newThreshold, final int newMaxIter,
+ final FieldEventHandler<? super FieldGeneralPenumbraDetector<T>, T> newHandler) {
+ return new FieldGeneralPenumbraDetector<>(newMaxCheck, newThreshold, newMaxIter, newHandler);
+ }
+
+ /** The G-function is the difference between the Sun-Sat-Central-Body angle and
+ * the central body apparent radius.
+ * @param s the current state information : date, kinematics, attitude
+ * @return value of the g function
+ */
+ public T g(final FieldSpacecraftState<T> s) {
+ final T[] angle = getGeneralEclipseAngles(s.getPVCoordinates().getPosition(),
+ provider.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
+ radius, sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
+ s.getA().getField().getZero().add(Constants.SUN_RADIUS));
+ return angle[0].subtract(angle[1]).subtract(angle[2]).add(ANGULAR_MARGIN);
+ }
+
+ }
}
diff --git a/src/main/java/org/orekit/forces/radiation/SolarRadiationPressure.java b/src/main/java/org/orekit/forces/radiation/SolarRadiationPressure.java
index ff00ce6bac0bceb608b66c8486851bef27a6825f..988c6c3010575ddabe7cb741eeb1b57a12af5b13 100644
--- a/src/main/java/org/orekit/forces/radiation/SolarRadiationPressure.java
+++ b/src/main/java/org/orekit/forces/radiation/SolarRadiationPressure.java
@@ -16,12 +16,15 @@
*/
package org.orekit.forces.radiation;
+import java.util.ArrayList;
import java.util.List;
+import java.util.Map;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
+import org.hipparchus.util.MathArrays;
import org.hipparchus.util.Precision;
import org.orekit.frames.Frame;
import org.orekit.propagation.FieldSpacecraftState;
@@ -109,7 +112,7 @@ public class SolarRadiationPressure extends AbstractRadiationForceModel {
final double r2 = sunSatVector.getNormSq();
// compute flux
- final double ratio = getLightingRatio(position, frame, date);
+ final double ratio = getTotalLightingRatio(position, frame, date);
final double rawP = ratio * kRef / r2;
final Vector3D flux = new Vector3D(rawP / FastMath.sqrt(r2), sunSatVector);
@@ -130,7 +133,7 @@ public class SolarRadiationPressure extends AbstractRadiationForceModel {
final T r2 = sunSatVector.getNormSq();
// compute flux
- final T ratio = getLightingRatio(position, frame, date);
+ final T ratio = getTotalLightingRatio(position, frame, date);
final T rawP = ratio.divide(r2).multiply(kRef);
final FieldVector3D<T> flux = new FieldVector3D<>(rawP.divide(r2.sqrt()), sunSatVector);
@@ -140,6 +143,7 @@ public class SolarRadiationPressure extends AbstractRadiationForceModel {
}
/** Get the lighting ratio ([0-1]).
+ * Considers only central body as occulting body.
* @param position the satellite's position in the selected frame.
* @param frame in which is defined the position
* @param date the date
@@ -201,7 +205,163 @@ public class SolarRadiationPressure extends AbstractRadiationForceModel {
}
+ /** Get eclipse ratio between to bodies seen from a specific object.
+ * Ratio is in [0-1].
+ * @param position the satellite's position in the selected frame
+ * @param occultingPosition the position of the occulting object
+ * @param occultingRadius the mean radius of the occulting object
+ * @param occultedPosition the position of the occulted object
+ * @param occultedRadius the mean radius of the occulted object
+ * @return eclipse ratio
+ */
+ public double getGeneralEclipseRatio(final Vector3D position,
+ final Vector3D occultingPosition,
+ final double occultingRadius,
+ final Vector3D occultedPosition,
+ final double occultedRadius) {
+
+
+ // Compute useful angles
+ final double[] angle = getGeneralEclipseAngles(position, occultingPosition, occultingRadius, occultedPosition, occultedRadius);
+
+ // Sat-Occulted/ Sat-Occulting angle
+ final double occultedSatOcculting = angle[0];
+
+ // Occulting apparent radius
+ final double alphaOcculting = angle[1];
+
+ // Occulted apparent radius
+ final double alphaOcculted = angle[2];
+
+ double result = 1.0;
+
+ // Is the satellite in complete umbra ?
+ if (occultedSatOcculting - alphaOcculting + alphaOcculted <= ANGULAR_MARGIN) {
+ result = 0.0;
+ } else if (occultedSatOcculting - alphaOcculting - alphaOcculted < -ANGULAR_MARGIN) {
+ // Compute an eclipse ratio in penumbra
+ final double sEA2 = occultedSatOcculting * occultedSatOcculting;
+ final double oo2sEA = 1.0 / (2. * occultedSatOcculting);
+ final double aS2 = alphaOcculted * alphaOcculted;
+ final double aE2 = alphaOcculting * alphaOcculting;
+ final double aE2maS2 = aE2 - aS2;
+
+ final double alpha1 = (sEA2 - aE2maS2) * oo2sEA;
+ final double alpha2 = (sEA2 + aE2maS2) * oo2sEA;
+
+ // Protection against numerical inaccuracy at boundaries
+ final double almost0 = Precision.SAFE_MIN;
+ final double almost1 = FastMath.nextDown(1.0);
+ final double a1oaS = FastMath.min(almost1, FastMath.max(-almost1, alpha1 / alphaOcculted));
+ final double aS2ma12 = FastMath.max(almost0, aS2 - alpha1 * alpha1);
+ final double a2oaE = FastMath.min(almost1, FastMath.max(-almost1, alpha2 / alphaOcculting));
+ final double aE2ma22 = FastMath.max(almost0, aE2 - alpha2 * alpha2);
+
+ final double P1 = aS2 * FastMath.acos(a1oaS) - alpha1 * FastMath.sqrt(aS2ma12);
+ final double P2 = aE2 * FastMath.acos(a2oaE) - alpha2 * FastMath.sqrt(aE2ma22);
+
+ result = 1. - (P1 + P2) / (FastMath.PI * aS2);
+ }
+
+ return result;
+ }
+
+ /** Get the total lighting ratio ([0-1]).
+ * This method considers every occulting bodies.
+ * @param position the satellite's position in the selected frame.
+ * @param frame in which is defined the position
+ * @param date the date
+ * @return lighting ratio
+ */
+ public double getTotalLightingRatio(final Vector3D position, final Frame frame, final AbsoluteDate date) {
+
+ double result = 0.0;
+ final Map<ExtendedPVCoordinatesProvider, Double> otherOccultingBodies = getOtherOccultingBodies();
+ final Vector3D sunPosition = sun.getPVCoordinates(date, frame).getPosition();
+ final int n = otherOccultingBodies.size() + 1;
+
+ if (n > 1) {
+
+ final Vector3D[] occultingBodyPositions = new Vector3D[n];
+ final double[] occultingBodyRadiuses = new double[n];
+
+ // Central body
+ occultingBodyPositions[0] = new Vector3D(0.0, 0.0, 0.0);
+ occultingBodyRadiuses[0] = getEquatorialRadius();
+
+ // Other occulting bodies
+ int k = 1;
+ for (ExtendedPVCoordinatesProvider provider: otherOccultingBodies.keySet()) {
+ occultingBodyPositions[k] = provider.getPVCoordinates(date, frame).getPosition();
+ occultingBodyRadiuses[k] = otherOccultingBodies.get(provider);
+ ++k;
+ }
+ for (int i = 0; i < n; ++i) {
+
+ // Lighting ratio computations
+ final double eclipseRatioI = getGeneralEclipseRatio(position,
+ occultingBodyPositions[i],
+ occultingBodyRadiuses[i],
+ sunPosition,
+ Constants.SUN_RADIUS);
+
+ // First body totaly occults Sun, full eclipse is occuring.
+ if (eclipseRatioI == 0.0) {
+ return 0.0;
+ }
+
+
+ result += eclipseRatioI;
+
+
+ // Mutual occulting body eclipse ratio computations between first and secondary bodies
+ for (int j = i + 1; j < n; ++j) {
+
+ final double eclipseRatioJ = getGeneralEclipseRatio(position,
+ occultingBodyPositions[j],
+ occultingBodyRadiuses[j],
+ sunPosition,
+ Constants.SUN_RADIUS);
+ final double eclipseRatioIJ = getGeneralEclipseRatio(position,
+ occultingBodyPositions[i],
+ occultingBodyRadiuses[i],
+ occultingBodyPositions[j],
+ occultingBodyRadiuses[j]);
+
+ final double alphaJ = getGeneralEclipseAngles(position,
+ occultingBodyPositions[i],
+ occultingBodyRadiuses[i],
+ occultingBodyPositions[j],
+ occultingBodyRadiuses[j])[2];
+
+ final double alphaSun = getEclipseAngles(sunPosition, position)[2];
+ final double alphaJSq = alphaJ * alphaJ;
+ final double alphaSunSq = alphaSun * alphaSun;
+ final double mutualEclipseCorrection = (1 - eclipseRatioIJ) * alphaJSq / alphaSunSq;
+
+ // Secondary body totaly occults Sun, no more computations are required, full eclipse is occuring.
+ if (eclipseRatioJ == 0.0 ) {
+ return 0.0;
+ }
+
+ // Secondary body partially occults Sun
+ else if (eclipseRatioJ != 1) {
+ result -= mutualEclipseCorrection;
+ }
+ }
+ }
+ // Final term
+ result -= n - 1;
+ } else {
+ // only central body is considered
+ result = getLightingRatio(position, frame, date);
+ }
+ return result;
+ }
+
+
/** Get the lighting ratio ([0-1]).
+ * Considers only central body as occulting body.
* @param position the satellite's position in the selected frame.
* @param frame in which is defined the position
* @param date the date
@@ -266,6 +426,166 @@ public class SolarRadiationPressure extends AbstractRadiationForceModel {
return result;
}
+ /** Get eclipse ratio between to bodies seen from a specific object.
+ * Ratio is in [0-1].
+ * @param position the satellite's position in the selected frame
+ * @param occultingPosition the position of the occulting object
+ * @param occultingRadius the mean radius of the occulting object
+ * @param occultedPosition the position of the occulted object
+ * @param occultedRadius the mean radius of the occulted object
+ * @param <T> extends RealFieldElement
+ * @return eclipse ratio
+ */
+ public <T extends CalculusFieldElement<T>> T getGeneralEclipseRatio(final FieldVector3D<T> position,
+ final FieldVector3D<T> occultingPosition,
+ final T occultingRadius,
+ final FieldVector3D<T> occultedPosition,
+ final T occultedRadius) {
+
+
+ final T one = occultingRadius.getField().getOne();
+
+ // Compute useful angles
+ final T[] angle = getGeneralEclipseAngles(position, occultingPosition, occultingRadius, occultedPosition, occultedRadius);
+
+ // Sat-Occulted/ Sat-Occulting angle
+ final T occultedSatOcculting = angle[0];
+
+ // Occulting apparent radius
+ final T alphaOcculting = angle[1];
+
+ // Occulted apparent radius
+ final T alphaOcculted = angle[2];
+
+ T result = one;
+
+ // Is the satellite in complete umbra ?
+ if (occultedSatOcculting.getReal() - alphaOcculting.getReal() + alphaOcculted.getReal() <= ANGULAR_MARGIN) {
+ result = occultingRadius.getField().getZero();
+ } else if (occultedSatOcculting.getReal() - alphaOcculting.getReal() - alphaOcculted.getReal() < -ANGULAR_MARGIN) {
+ // Compute a lighting ratio in penumbra
+ final T sEA2 = occultedSatOcculting.multiply(occultedSatOcculting);
+ final T oo2sEA = occultedSatOcculting.multiply(2).reciprocal();
+ final T aS2 = alphaOcculted.multiply(alphaOcculted);
+ final T aE2 = alphaOcculting.multiply(alphaOcculting);
+ final T aE2maS2 = aE2.subtract(aS2);
+
+ final T alpha1 = sEA2.subtract(aE2maS2).multiply(oo2sEA);
+ final T alpha2 = sEA2.add(aE2maS2).multiply(oo2sEA);
+
+ // Protection against numerical inaccuracy at boundaries
+ final double almost0 = Precision.SAFE_MIN;
+ final double almost1 = FastMath.nextDown(1.0);
+ final T a1oaS = min(almost1, max(-almost1, alpha1.divide(alphaOcculted)));
+ final T aS2ma12 = max(almost0, aS2.subtract(alpha1.multiply(alpha1)));
+ final T a2oaE = min(almost1, max(-almost1, alpha2.divide(alphaOcculting)));
+ final T aE2ma22 = max(almost0, aE2.subtract(alpha2.multiply(alpha2)));
+
+ final T P1 = aS2.multiply(a1oaS.acos()).subtract(alpha1.multiply(aS2ma12.sqrt()));
+ final T P2 = aE2.multiply(a2oaE.acos()).subtract(alpha2.multiply(aE2ma22.sqrt()));
+
+ result = one.subtract(P1.add(P2).divide(aS2.multiply(FastMath.PI)));
+ }
+
+ return result;
+ }
+
+ /** Get the total lighting ratio ([0-1]).
+ * This method considers every occulting bodies.
+ * @param position the satellite's position in the selected frame.
+ * @param frame in which is defined the position
+ * @param date the date
+ * @param <T> extends RealFieldElement
+ * @return lighting rati
+ */
+ public <T extends CalculusFieldElement<T>> T getTotalLightingRatio(final FieldVector3D<T> position, final Frame frame, final FieldAbsoluteDate<T> date) {
+
+ final T zero = position.getAlpha().getField().getZero();
+ T result = zero;
+ final Map<ExtendedPVCoordinatesProvider, Double> otherOccultingBodies = getOtherOccultingBodies();
+ final FieldVector3D<T> sunPosition = sun.getPVCoordinates(date, frame).getPosition();
+ final int n = otherOccultingBodies.size() + 1;
+
+ if (n > 1) {
+
+ final List<FieldVector3D<T>> occultingBodyPositions = new ArrayList<FieldVector3D<T>>(n);
+ final T[] occultingBodyRadiuses = MathArrays.buildArray(zero.getField(), n);
+
+ // Central body
+ occultingBodyPositions.add(0, new FieldVector3D<>(zero, zero, zero));
+ occultingBodyRadiuses[0] = zero.add(getEquatorialRadius());
+
+ // Other occulting bodies
+ int k = 1;
+ for (ExtendedPVCoordinatesProvider provider: otherOccultingBodies.keySet()) {
+ occultingBodyPositions.add(k, provider.getPVCoordinates(date, frame).getPosition());
+ occultingBodyRadiuses[k] = zero.add(otherOccultingBodies.get(provider));
+ ++k;
+ }
+ for (int i = 0; i < n; ++i) {
+
+ // Lighting ratio computations
+ final T eclipseRatioI = getGeneralEclipseRatio(position,
+ occultingBodyPositions.get(i),
+ occultingBodyRadiuses[i],
+ sunPosition,
+ zero.add(Constants.SUN_RADIUS));
+
+ // First body totaly occults Sun, full eclipse is occuring.
+ if (eclipseRatioI.getReal() == 0.0) {
+ return zero;
+ }
+
+
+ result = result.add(eclipseRatioI);
+
+
+ // Mutual occulting body eclipse ratio computations between first and secondary bodies
+ for (int j = i + 1; j < n; ++j) {
+
+ final T eclipseRatioJ = getGeneralEclipseRatio(position,
+ occultingBodyPositions.get(i),
+ occultingBodyRadiuses[j],
+ sunPosition,
+ zero.add(Constants.SUN_RADIUS));
+ final T eclipseRatioIJ = getGeneralEclipseRatio(position,
+ occultingBodyPositions.get(i),
+ occultingBodyRadiuses[i],
+ occultingBodyPositions.get(j),
+ occultingBodyRadiuses[j]);
+
+ final T alphaJ = getGeneralEclipseAngles(position,
+ occultingBodyPositions.get(i),
+ occultingBodyRadiuses[i],
+ occultingBodyPositions.get(j),
+ occultingBodyRadiuses[j])[2];
+
+ final T alphaSun = getEclipseAngles(sunPosition, position)[2];
+ final T alphaJSq = alphaJ.multiply(alphaJ);
+ final T alphaSunSq = alphaSun.multiply(alphaSun);
+ final T mutualEclipseCorrection = eclipseRatioIJ.negate().add(1).multiply(alphaJSq).divide(alphaSunSq);
+
+ // Secondary body totaly occults Sun, no more computations are required, full eclipse is occuring.
+ if (eclipseRatioJ.getReal() == 0.0 ) {
+ return zero;
+ }
+
+ // Secondary body partially occults Sun
+ else if (eclipseRatioJ.getReal() != 1) {
+ result = result.subtract(mutualEclipseCorrection);
+ }
+ }
+ }
+ // Final term
+ result = result.subtract(n - 1);
+ } else {
+ // only central body is considered
+ result = getLightingRatio(position, frame, date);
+ }
+ return result;
+ }
+
+
/** {@inheritDoc} */
@Override
public List<ParameterDriver> getParametersDrivers() {
@@ -279,7 +599,7 @@ public class SolarRadiationPressure extends AbstractRadiationForceModel {
* @return min value
*/
private <T extends CalculusFieldElement<T>> T min(final double d, final T f) {
- return (f.getReal() > d) ? f.getField().getZero().add(d) : f;
+ return (f.getReal() > d) ? f.getField().getZero().newInstance(d) : f;
}
/** Compute max of two values, one double and one field element.
@@ -289,7 +609,7 @@ public class SolarRadiationPressure extends AbstractRadiationForceModel {
* @return max value
*/
private <T extends CalculusFieldElement<T>> T max(final double d, final T f) {
- return (f.getReal() <= d) ? f.getField().getZero().add(d) : f;
+ return (f.getReal() <= d) ? f.getField().getZero().newInstance(d) : f;
}
}
diff --git a/src/test/java/org/orekit/forces/radiation/SolarRadiationPressureTest.java b/src/test/java/org/orekit/forces/radiation/SolarRadiationPressureTest.java
index 247b44385b25c8621c707e2a2c78bccb709728a5..166dcddd8d03ce30af1dca19cd640c4fb7ea23bf 100644
--- a/src/test/java/org/orekit/forces/radiation/SolarRadiationPressureTest.java
+++ b/src/test/java/org/orekit/forces/radiation/SolarRadiationPressureTest.java
@@ -17,6 +17,7 @@
package org.orekit.forces.radiation;
+import java.io.File;
import java.io.FileNotFoundException;
import java.lang.reflect.InvocationTargetException;
import java.text.ParseException;
@@ -41,6 +42,9 @@ import org.orekit.Utils;
import org.orekit.attitudes.LofOffset;
import org.orekit.bodies.CelestialBodyFactory;
import org.orekit.bodies.OneAxisEllipsoid;
+import org.orekit.data.DataContext;
+import org.orekit.data.DataProvidersManager;
+import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.forces.AbstractLegacyForceModelTest;
import org.orekit.forces.BoxAndSolarArraySpacecraft;
@@ -175,7 +179,7 @@ public class SolarRadiationPressureTest extends AbstractLegacyForceModelTest {
ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();
SolarRadiationPressure srp =
new SolarRadiationPressure(sun, Constants.SUN_RADIUS,
- (RadiationSensitive) new IsotropicRadiationClassicalConvention(50.0, 0.5, 0.5));
+ new IsotropicRadiationClassicalConvention(50.0, 0.5, 0.5));
Assert.assertFalse(srp.dependsOnPositionOnly());
Vector3D position = orbit.getPVCoordinates().getPosition();
@@ -206,7 +210,7 @@ public class SolarRadiationPressureTest extends AbstractLegacyForceModelTest {
FramesFactory.getITRF(IERSConventions.IERS_2010, true));
SolarRadiationPressure SRP =
new SolarRadiationPressure(sun, earth.getEquatorialRadius(),
- (RadiationSensitive) new IsotropicRadiationCNES95Convention(50.0, 0.5, 0.5));
+ new IsotropicRadiationCNES95Convention(50.0, 0.5, 0.5));
double period = 2*FastMath.PI*FastMath.sqrt(orbit.getA()*orbit.getA()*orbit.getA()/orbit.getMu());
Assert.assertEquals(86164, period, 1);
@@ -603,6 +607,7 @@ public class SolarRadiationPressureTest extends AbstractLegacyForceModelTest {
private static class SolarStepHandler implements OrekitFixedStepHandler {
+ @Override
public void handleStep(SpacecraftState currentState, boolean isLast) {
final double dex = currentState.getEquinoctialEx() - 0.01071166;
final double dey = currentState.getEquinoctialEy() - 0.00654848;
@@ -656,7 +661,7 @@ public class SolarRadiationPressureTest extends AbstractLegacyForceModelTest {
AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator =
new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
integrator.setInitialStepSize(60);
-
+
AdaptiveStepsizeIntegrator RIntegrator =
new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
RIntegrator.setInitialStepSize(60);
@@ -682,7 +687,7 @@ public class SolarRadiationPressureTest extends AbstractLegacyForceModelTest {
FNP.addForceModel(forceModel);
NP.addForceModel(forceModel);
-
+
// Do the test
checkRealFieldPropagation(FKO, PositionAngle.MEAN, 1000., NP, FNP,
1.0e-30, 5.0e-10, 3.0e-11, 3.0e-10,
@@ -761,6 +766,129 @@ public class SolarRadiationPressureTest extends AbstractLegacyForceModelTest {
Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getZ() - finPVC_R.getPosition().getZ()) < FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
}
+ /** Testing if eclipses due to Moon are considered.
+ * Earth is artificially reduced to a single point to only consider Moon effect.
+ * Reference values are presented in "A Study of Solar Radiation Pressure acting on GPS Satellites"
+ * written by Laurent Olivier Froideval in 2009.
+ * Modifications of the step handler and time span able to print lighting ratios other a year and get a reference like graph.
+ */
+ @Test
+ public void testMoonEclipse() {
+
+ // Configure Orekit
+ final File home = new File(System.getProperty("user.home"));
+ final File orekitData = new File(home, "orekit-data");
+ final DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
+ manager.addProvider(new DirectoryCrawler(orekitData));
+
+ final ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();
+ final ExtendedPVCoordinatesProvider moon = CelestialBodyFactory.getMoon();
+ final Frame GCRF = FramesFactory.getGCRF();
+ final AbsoluteDate date = new AbsoluteDate(2007, 01, 01, 6, 0, 0, TimeScalesFactory.getGPS());
+ final double dt = 3600 * 24 * 180;
+ final AbsoluteDate target = date.shiftedBy(dt);
+
+ //PV coordinates from sp3 file esa14081.sp3 PRN 03
+ final PVCoordinates refPV = new PVCoordinates(new Vector3D(14986728.76145754, 14849687.258579938, -16523319.786690142),
+ new Vector3D(-3152.6722260146, 1005.6757515113, -1946.9273038773));
+ final Orbit orbit = new CartesianOrbit(refPV, GCRF, date, Constants.EIGEN5C_EARTH_MU);
+ final SpacecraftState initialState = new SpacecraftState(orbit);
+
+ // Create SRP perturbation with Moon and Earth
+ SolarRadiationPressure srp =
+ new SolarRadiationPressure(sun, Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
+ new IsotropicRadiationClassicalConvention(50.0, 0.5, 0.5));
+ srp.addOccultingBody(moon, Constants.MOON_EQUATORIAL_RADIUS);
+
+ // creation of the propagator
+ double[] absTolerance = {
+ 0.1, 1.0e-9, 1.0e-9, 1.0e-5, 1.0e-5, 1.0e-5, 0.001
+ };
+ double[] relTolerance = {
+ 1.0e-4, 1.0e-4, 1.0e-4, 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-7
+ };
+ AdaptiveStepsizeIntegrator integrator =
+ new DormandPrince853Integrator(600, 1000000, absTolerance, relTolerance);
+ final NumericalPropagator propagator = new NumericalPropagator(integrator);
+ propagator.setInitialState(initialState);
+ propagator.addForceModel(srp);
+ propagator.setMasterMode(600, new MoonEclipseStepHandler(moon, sun, srp));
+ propagator.propagate(target);
+ }
+
+ /** Specialized step handler.
+ * <p>This class extends the step handler in order to print on the output stream at the given step.<p>
+ * @author Thomas Paulet
+ */
+ private static class MoonEclipseStepHandler implements OrekitFixedStepHandler {
+
+ final ExtendedPVCoordinatesProvider moon;
+ final ExtendedPVCoordinatesProvider sun;
+ final SolarRadiationPressure srp;
+
+ /** Simple constructor.
+ */
+ MoonEclipseStepHandler(final ExtendedPVCoordinatesProvider moon, final ExtendedPVCoordinatesProvider sun,
+ final SolarRadiationPressure srp) {
+ this.moon = moon;
+ this.sun = sun;
+ this.srp = srp;
+
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public void init(final SpacecraftState s0, final AbsoluteDate t, final double step) {
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public void handleStep(final SpacecraftState currentState, final boolean isLast) {
+ final AbsoluteDate date = currentState.getDate();
+ final Frame frame = currentState.getFrame();
+
+ final Vector3D moonPos = moon.getPVCoordinates(date, frame).getPosition();
+ final Vector3D sunPos = sun.getPVCoordinates(date, frame).getPosition();
+ final Vector3D statePos = currentState.getPVCoordinates().getPosition();
+
+ // Moon umbra and penumbra conditions
+ final double[] moonAngles = srp.getGeneralEclipseAngles(statePos, moonPos, Constants.MOON_EQUATORIAL_RADIUS,
+ sunPos, Constants.SUN_RADIUS);
+ final double moonUmbra = moonAngles[0] - moonAngles[1] + moonAngles[2];
+ final boolean isInMoonUmbra = (moonUmbra < 1.0e-10);
+
+ final double moonPenumbra = moonAngles[0] - moonAngles[1] - moonAngles[2];
+ final boolean isInMoonPenumbra = (moonPenumbra < -1.0e-10);
+
+ // Earth umbra and penumbra conditions
+ final double[] earthAngles = srp.getEclipseAngles(sunPos, statePos);
+
+ final double earthUmbra = earthAngles[0] - earthAngles[1] + earthAngles[2];
+ final boolean isInEarthUmbra = (earthUmbra < 1.0e-10);
+
+ final double earthPenumbra = earthAngles[0] - earthAngles[1] - earthAngles[2];
+ final boolean isInEarthPenumbra = (earthPenumbra < -1.0e-10);
+
+
+ // Compute lighting ration
+ final double lightingRatio = srp.getTotalLightingRatio(statePos, frame, date);
+
+ // Check behaviour
+ if (isInMoonUmbra || isInEarthUmbra) {
+ Assert.assertEquals(0.0, lightingRatio, 1e-8);
+ }
+ else if (isInMoonPenumbra || isInEarthPenumbra) {
+ Assert.assertEquals(true, (lightingRatio < 1.0));
+ Assert.assertEquals(true, (lightingRatio > 0.0));
+ }
+ else {
+ Assert.assertEquals(1.0, lightingRatio, 1e-8);
+ }
+ }
+ }
+
+
+
@Before
public void setUp() {
Utils.setDataRoot("regular-data");