1   package org.orekit.forces.radiation;
2   
3   import org.hamcrest.MatcherAssert;
4   import org.hamcrest.Matchers;
5   import org.hamcrest.comparator.ComparatorMatcherBuilder;
6   import org.hipparchus.Field;
7   import org.hipparchus.complex.Complex;
8   import org.hipparchus.complex.ComplexField;
9   import org.hipparchus.dfp.Dfp;
10  import org.hipparchus.dfp.DfpField;
11  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
12  import org.hipparchus.geometry.euclidean.threed.Vector3D;
13  import org.hipparchus.util.FastMath;
14  import org.junit.jupiter.api.Assertions;
15  import org.junit.jupiter.api.BeforeEach;
16  import org.junit.jupiter.api.Test;
17  import org.orekit.OrekitMatchers;
18  import org.orekit.Utils;
19  import org.orekit.bodies.AnalyticalSolarPositionProvider;
20  import org.orekit.bodies.OneAxisEllipsoid;
21  import org.orekit.errors.OrekitException;
22  import org.orekit.frames.Frame;
23  import org.orekit.frames.FramesFactory;
24  import org.orekit.orbits.CartesianOrbit;
25  import org.orekit.propagation.FieldSpacecraftState;
26  import org.orekit.propagation.SpacecraftState;
27  import org.orekit.propagation.events.EventDetectionSettings;
28  import org.orekit.propagation.events.EventDetector;
29  import org.orekit.propagation.events.FieldEventDetector;
30  import org.orekit.propagation.events.handlers.FieldResetDerivativesOnEvent;
31  import org.orekit.propagation.events.handlers.ResetDerivativesOnEvent;
32  import org.orekit.time.AbsoluteDate;
33  import org.orekit.utils.Constants;
34  import org.orekit.utils.ExtendedPositionProvider;
35  import org.orekit.utils.PVCoordinates;
36  
37  import java.util.Comparator;
38  import java.util.List;
39  
40  class ConicallyShadowedLightFluxModelTest {
41  
42      /** Compares Dfp values. */
43      private static final Comparator<Dfp> dfpComparator = (o1, o2) -> {
44          if (o1.lessThan(o2)) {
45              return -1;
46          }
47          if (o1.greaterThan(o2)) {
48              return 1;
49          }
50          return 0;
51      };
52  
53      @BeforeEach
54      public void setUp() throws OrekitException {
55          Utils.setDataRoot("potential/shm-format:regular-data");
56      }
57  
58      @Test
59      void testConstructor() {
60          // GIVEN
61          final ExtendedPositionProvider positionProvider = new AnalyticalSolarPositionProvider();
62          final ConicallyShadowedLightFluxModel fluxModel = new ConicallyShadowedLightFluxModel(1., 1.,
63                  positionProvider, 1.);
64          // WHEN
65          final List<EventDetector> detectors = fluxModel.getEclipseConditionsDetector();
66          // THEN
67          final EventDetectionSettings detectionSettings = ConicallyShadowedLightFluxModel.getDefaultEclipseDetectionSettings();
68          Assertions.assertEquals(detectionSettings.getMaxIterationCount(), detectors.get(0).getMaxIterationCount());
69          Assertions.assertEquals(detectionSettings.getThreshold(), detectors.get(0).getThreshold());
70      }
71  
72      @Test
73      void testGetEclipseConditionsDetector() {
74          // GIVEN
75          final ExtendedPositionProvider positionProvider = new AnalyticalSolarPositionProvider();
76          final ConicallyShadowedLightFluxModel fluxModel = new ConicallyShadowedLightFluxModel(1., positionProvider,
77                  1.);
78          // WHEN
79          final List<EventDetector> detectors = fluxModel.getEclipseConditionsDetector();
80          // THEN
81          Assertions.assertEquals(2, detectors.size());
82          for (EventDetector detector : detectors) {
83              Assertions.assertInstanceOf(ResetDerivativesOnEvent.class, detector.getHandler());
84          }
85      }
86  
87      @Test
88      void testGetFieldEclipseConditionsDetector() {
89          // GIVEN
90          final ExtendedPositionProvider positionProvider = new AnalyticalSolarPositionProvider();
91          final ConicallyShadowedLightFluxModel fluxModel = new ConicallyShadowedLightFluxModel(1., positionProvider,
92                  1.);
93          final ComplexField field = ComplexField.getInstance();
94          // WHEN
95          final List<FieldEventDetector<Complex>> fieldDetectors = fluxModel.getFieldEclipseConditionsDetector(field);
96          // THEN
97          final List<EventDetector> detectors = fluxModel.getEclipseConditionsDetector();
98          Assertions.assertEquals(detectors.size(), fieldDetectors.size());
99          for (int i = 0; i < detectors.size(); i++) {
100             Assertions.assertInstanceOf(FieldResetDerivativesOnEvent.class, fieldDetectors.get(i).getHandler());
101             Assertions.assertEquals(detectors.get(i).getThreshold(), fieldDetectors.get(i).getThreshold().getReal());
102             Assertions.assertEquals(detectors.get(i).getMaxIterationCount(), fieldDetectors.get(i).getMaxIterationCount());
103         }
104         final SpacecraftState state = new SpacecraftState(new CartesianOrbit(new PVCoordinates(Vector3D.PLUS_I, Vector3D.MINUS_J),
105                 FramesFactory.getEME2000(), AbsoluteDate.ARBITRARY_EPOCH, 1.));
106         final FieldSpacecraftState<Complex> fieldState = new FieldSpacecraftState<>(field, state);
107         fluxModel.init(fieldState, null);
108         for (int i = 0; i < detectors.size(); i++) {
109             Assertions.assertEquals(detectors.get(i).g(state), fieldDetectors.get(i).g(fieldState).getReal(), 1e-8);
110             Assertions.assertEquals(detectors.get(i).getMaxCheckInterval().currentInterval(state, true),
111                     fieldDetectors.get(i).getMaxCheckInterval().currentInterval(fieldState, true), 1e-14);
112         }
113     }
114 
115     @Test
116     void testFieldGetLightingRatio() {
117         // GIVEN
118         final ExtendedPositionProvider positionProvider = new AnalyticalSolarPositionProvider();
119         final ConicallyShadowedLightFluxModel fluxModel = new ConicallyShadowedLightFluxModel(Constants.SUN_RADIUS, positionProvider,
120                 Constants.EGM96_EARTH_EQUATORIAL_RADIUS);
121         final ComplexField field = ComplexField.getInstance();
122         // WHEN & THEN
123         for (int i = 0; i < 100000; i++) {
124             final SpacecraftState state = new SpacecraftState(new CartesianOrbit(
125                     new PVCoordinates(new Vector3D(fluxModel.getOccultingBodyRadius() + 500e3, 0., 1e3), new Vector3D(1e1, 7e3, 1e2)),
126                     FramesFactory.getEME2000(), AbsoluteDate.ARBITRARY_EPOCH.shiftedBy(i * 3600.), Constants.EGM96_EARTH_MU));
127             final FieldSpacecraftState<Complex> fieldState = new FieldSpacecraftState<>(field, state);
128             final double expectedRatio = fluxModel.getLightingRatio(state);
129             Assertions.assertEquals(expectedRatio, fluxModel.getLightingRatio(fieldState).getReal(), 1e-8);
130         }
131     }
132 
133     @Test
134     void testGetLightingRatio() {
135         // GIVEN
136         final ExtendedPositionProvider positionProvider = new AnalyticalSolarPositionProvider();
137         final ConicallyShadowedLightFluxModel fluxModel = new ConicallyShadowedLightFluxModel(Constants.SUN_RADIUS, positionProvider,
138                 Constants.EGM96_EARTH_EQUATORIAL_RADIUS);
139         final Frame frame = FramesFactory.getGCRF();
140         final SolarRadiationPressure radiationPressure = new SolarRadiationPressure(positionProvider,
141                 new OneAxisEllipsoid(fluxModel.getOccultingBodyRadius(), 0., FramesFactory.getGCRF()), null);
142         // WHEN & THEN
143         for (int i = 0; i < 10000; i++) {
144             final SpacecraftState state = new SpacecraftState(new CartesianOrbit(
145                     new PVCoordinates(new Vector3D(fluxModel.getOccultingBodyRadius() + 500e3, 0., 1e3), new Vector3D(1e1, 7e3, 1e2)),
146                     frame, AbsoluteDate.ARBITRARY_EPOCH.shiftedBy(i * 3600.), Constants.EGM96_EARTH_MU));
147             final double expectedRatio = radiationPressure.getLightingRatio(state);
148             Assertions.assertEquals(expectedRatio, fluxModel.getLightingRatio(state), 1e-6);
149         }
150     }
151 
152     @Test
153     public void testGetLightingRatioConditions() {
154         // setup
155         final double au = Constants.JPL_SSD_ASTRONOMICAL_UNIT;
156         Vector3D sunP = new Vector3D(-au, 0, 0);
157         final double rS = 6.957E8;
158         final double rE = 6378136.3;
159         ConicallyShadowedLightFluxModel model =
160                 new ConicallyShadowedLightFluxModel(0, rS, null, rE);
161         final double f2 = FastMath.asin((rS - rE) / au);
162         final double c2 = au - rE / FastMath.sin(f2);
163         final double l2 = c2 * FastMath.tan(f2);
164         // small angle approximation for ratio of disk area
165         final double expected = 1.0 - 4 * rE * rE / (rS * rS);
166 
167         // action
168         MatcherAssert.assertThat(
169                 model.getLightingRatio(new Vector3D(au, 0, 0), sunP),
170                 Matchers.allOf(
171                         Matchers.greaterThanOrEqualTo(0.0),
172                         Matchers.lessThanOrEqualTo(1.0),
173                         Matchers.closeTo(expected, 1e-8)));
174         // on the edge
175         MatcherAssert.assertThat(
176                 model.getLightingRatio(new Vector3D(au, l2, 0), sunP),
177                 Matchers.allOf(
178                         Matchers.greaterThanOrEqualTo(0.0),
179                         Matchers.lessThanOrEqualTo(1.0),
180                         Matchers.greaterThanOrEqualTo(expected),
181                         Matchers.closeTo(expected, 1e-8)));
182         // disks intersect
183         MatcherAssert.assertThat(
184                 model.getLightingRatio(new Vector3D(au, l2 + 1e6, 0), sunP),
185                 Matchers.allOf(
186                         Matchers.greaterThanOrEqualTo(0.0),
187                         Matchers.lessThanOrEqualTo(1.0),
188                         Matchers.greaterThanOrEqualTo(expected),
189                         Matchers.closeTo(expected, 1e-5)));
190     }
191 
192     @Test
193     public void testGetLightingRatioConditionsField() {
194         // setup
195         Field<Dfp> field = new DfpField(20);
196         final Dfp zero = field.getZero();
197         final Dfp one = field.getOne();
198         final Dfp au = zero.add(Constants.JPL_SSD_ASTRONOMICAL_UNIT);
199         FieldVector3D<Dfp> sunP = new FieldVector3D<>(au.negate(), zero, zero);
200         final Dfp rS = zero.add(6.957E8);
201         final Dfp rE = zero.add(6378136.3);
202         ConicallyShadowedLightFluxModel model =
203                 new ConicallyShadowedLightFluxModel(0, rS.getReal(), null, rE.getReal());
204         final Dfp f2 = FastMath.asin((rS.subtract(rE)).divide(au));
205         final Dfp c2 = au.subtract(rE.divide(FastMath.sin(f2)));
206         final Dfp l2 = c2.multiply(FastMath.tan(f2));
207         // small angle approximation for ratio of disk area
208         final Dfp expected = one.subtract(rE.square().divide(rS.square()).multiply(4.0));
209 
210         // action
211         final ComparatorMatcherBuilder<Dfp> builder =
212                 ComparatorMatcherBuilder.comparedBy(dfpComparator);
213         MatcherAssert.assertThat(
214                 model.getLightingRatio(new FieldVector3D<>(au, zero, zero), sunP),
215                 Matchers.allOf(
216                         builder.greaterThanOrEqualTo(zero),
217                         builder.lessThanOrEqualTo(one),
218                         OrekitMatchers.closeTo(expected, 1e-9)));
219         // on the edge
220         MatcherAssert.assertThat(
221                 model.getLightingRatio(new FieldVector3D<>(au, l2, zero), sunP),
222                 Matchers.allOf(
223                         builder.greaterThanOrEqualTo(zero),
224                         builder.lessThanOrEqualTo(one),
225                         builder.greaterThanOrEqualTo(expected),
226                         OrekitMatchers.closeTo(expected, 1e-8)));
227         // disks intersect
228         MatcherAssert.assertThat(
229                 model.getLightingRatio(new FieldVector3D<>(au, l2.add(1e6), zero), sunP),
230                 Matchers.allOf(
231                         builder.greaterThanOrEqualTo(zero),
232                         builder.lessThanOrEqualTo(one),
233                         builder.greaterThanOrEqualTo(expected),
234                         OrekitMatchers.closeTo(expected, 1e-5)));
235     }
236 
237     /**
238      * Check for NaN in case discovered when running RadiationPressureModelTest.
239      * This case causes numerical issues because the Sun is almost entirely
240      * occulted. The remaining area is computed as 1 - (occulted area) which
241      * causes catastrophic cancellation in this case. Computing the occulted
242      * area is further problematic because it sums two quantities that are
243      * almost the same magnitude, but with opposite signs. Even in extended
244      * precision, changing nominally equivalent expressions can produce a
245      * negative (infeasible) result. For example, computing an angle with acos
246      * instead of atan2.
247      */
248     @Test
249     public void testGetLightingRatioNan() {
250         // setup
251         final ConicallyShadowedLightFluxModel model = new ConicallyShadowedLightFluxModel(
252                 1.0205062355092827E17,
253                 6.957E8,
254                 null,
255                 6378136.3);
256         final Vector3D position = new Vector3D(
257                 4634793.393351071,
258                 6979388.527107593,
259                 0.1946021760969457);
260         final Vector3D occultedBodyPosition = new Vector3D(
261                 2.6535624002170452E10,
262                 -1.3275125037266681E11,
263                 -5.755404493076553E10);
264         // value was computed with extended precision below.
265         // Though I don't trust it, making some mathematically allowable
266         // transformations gives different, even negative results.
267         DfpField field = new DfpField(50);
268         final Dfp zero = field.getZero();
269         final Dfp expected = zero.newInstance("7.0542775362244865347040611195042e-21");
270 
271         // action
272         double actual = model.getLightingRatio(position, occultedBodyPosition);
273 
274         // verify
275         MatcherAssert.assertThat(actual, Matchers.greaterThanOrEqualTo(0.0));
276         MatcherAssert.assertThat(actual,
277                 Matchers.closeTo(expected.getReal(), 1e-5));
278 
279         // now for the field version
280         FieldVector3D<Dfp> fieldPosition = new FieldVector3D<>(field, position);
281         FieldVector3D<Dfp> fieldOccultedBodyPosition =
282                 new FieldVector3D<>(field, occultedBodyPosition);
283         Dfp fieldActual = model.getLightingRatio(
284                 fieldPosition,
285                 fieldOccultedBodyPosition);
286         final ComparatorMatcherBuilder<Dfp> builder =
287                 ComparatorMatcherBuilder.comparedBy(dfpComparator);
288         MatcherAssert.assertThat(
289                 fieldActual,
290                 builder.greaterThanOrEqualTo(zero));
291         MatcherAssert.assertThat(
292                 fieldActual,
293                 OrekitMatchers.closeTo(expected, 1e-25));
294     }
295 
296     /**
297      * Check another case that produced NaNs. Satellite is on the line between
298      * penumbra and umbra. #1892.
299      */
300     @Test
301     public void testGetLightingRatioNanAgain() {
302         // setup
303         final ConicallyShadowedLightFluxModel model = new ConicallyShadowedLightFluxModel(
304                 1.0205062355092827E17,
305                 6.957E8,
306                 null,
307                 6378136.3);
308         final Vector3D position = new Vector3D(
309                 5416037.98611839,
310                 4557004.944798493,
311                 2.09861286566582);
312         final Vector3D occultedBodyPosition = new Vector3D(
313                 2.7402005760327873E10,
314                 -1.3260269763721999E11,
315                 -5.748964657098426E10);
316         // point is right on the line between umbra and penumbra, so depending
317         // on the check used it should either be zero, or a very small positive
318         // value.
319         final double expected = 0.0;
320 
321         // action
322         double actual = model.getLightingRatio(position, occultedBodyPosition);
323 
324         // verify
325         MatcherAssert.assertThat(actual, Matchers.greaterThanOrEqualTo(0.0));
326         MatcherAssert.assertThat(actual, Matchers.closeTo(expected, 0.0));
327 
328         // now for the field version
329         Field<Dfp> field = new DfpField(50);
330         FieldVector3D<Dfp> fieldPosition = new FieldVector3D<>(field, position);
331         FieldVector3D<Dfp> fieldOccultedBodyPosition =
332                 new FieldVector3D<>(field, occultedBodyPosition);
333         Dfp fieldActual = model.getLightingRatio(
334                 fieldPosition,
335                 fieldOccultedBodyPosition);
336         final Dfp zero = field.getZero();
337         final ComparatorMatcherBuilder<Dfp> builder =
338                 ComparatorMatcherBuilder.comparedBy(dfpComparator);
339         MatcherAssert.assertThat(
340                 fieldActual,
341                 builder.greaterThanOrEqualTo(zero));
342         MatcherAssert.assertThat(
343                 fieldActual,
344                 OrekitMatchers.closeTo(zero.add(expected), 0.0));
345     }
346 
347 }