1   package org.orekit.forces.radiation;
2   
3   import org.hipparchus.CalculusFieldElement;
4   import org.hipparchus.Field;
5   import org.hipparchus.complex.Complex;
6   import org.hipparchus.complex.ComplexField;
7   import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
8   import org.hipparchus.geometry.euclidean.threed.Vector3D;
9   import org.hipparchus.ode.FieldODEIntegrator;
10  import org.hipparchus.ode.ODEIntegrator;
11  import org.hipparchus.util.Binary64;
12  import org.hipparchus.util.Binary64Field;
13  import org.junit.jupiter.api.Assertions;
14  import org.junit.jupiter.api.Test;
15  import org.junit.jupiter.params.ParameterizedTest;
16  import org.junit.jupiter.params.provider.ValueSource;
17  import org.mockito.Mockito;
18  import org.orekit.Utils;
19  import org.orekit.bodies.CelestialBodyFactory;
20  import org.orekit.bodies.OneAxisEllipsoid;
21  import org.orekit.forces.BoxAndSolarArraySpacecraft;
22  import org.orekit.forces.ForceModel;
23  import org.orekit.frames.FramesFactory;
24  import org.orekit.orbits.EquinoctialOrbit;
25  import org.orekit.orbits.Orbit;
26  import org.orekit.orbits.OrbitType;
27  import org.orekit.orbits.PositionAngleType;
28  import org.orekit.propagation.FieldSpacecraftState;
29  import org.orekit.propagation.SpacecraftState;
30  import org.orekit.propagation.conversion.DormandPrince54FieldIntegratorBuilder;
31  import org.orekit.propagation.conversion.DormandPrince54IntegratorBuilder;
32  import org.orekit.propagation.conversion.FieldODEIntegratorBuilder;
33  import org.orekit.propagation.conversion.ODEIntegratorBuilder;
34  import org.orekit.propagation.events.EventDetector;
35  import org.orekit.propagation.events.FieldEventDetector;
36  import org.orekit.propagation.numerical.FieldNumericalPropagator;
37  import org.orekit.propagation.numerical.NumericalPropagator;
38  import org.orekit.time.AbsoluteDate;
39  import org.orekit.time.FieldAbsoluteDate;
40  import org.orekit.utils.Constants;
41  import org.orekit.utils.ExtendedPositionProvider;
42  
43  import java.util.ArrayList;
44  import java.util.List;
45  import java.util.stream.Stream;
46  
47  class RadiationPressureModelTest {
48  
49      @Test
50      void testAcceleration() {
51          // GIVEN
52          final SpacecraftState state = createState(42000e3);
53          final ComplexField field = ComplexField.getInstance();
54          final FieldSpacecraftState<Complex> fieldState = new FieldSpacecraftState<>(field, state);
55          final IsotropicRadiationSingleCoefficient radiationSingleCoefficient = new IsotropicRadiationSingleCoefficient(1., 2.);
56          final LightFluxModel mockedFluxModel = Mockito.mock(LightFluxModel.class);
57          final Vector3D fluxVector = new Vector3D(1., 2., 3.);
58          Mockito.when(mockedFluxModel.getLightFluxVector(state)).thenReturn(fluxVector);
59          Mockito.when(mockedFluxModel.getLightFluxVector(fieldState)).thenReturn(new FieldVector3D<>(field, fluxVector));
60          final RadiationPressureModel forceModel = new RadiationPressureModel(mockedFluxModel,
61                  radiationSingleCoefficient);
62          // WHEN
63          final FieldVector3D<Complex> fieldAcceleration = forceModel.acceleration(fieldState, forceModel.getParameters(field));
64          // THEN
65          final Vector3D expectedAcceleration = forceModel.acceleration(state, forceModel.getParameters());
66          Assertions.assertEquals(expectedAcceleration, fieldAcceleration.toVector3D());
67      }
68  
69      @Test
70      void testDependsOnPositionOnlyTrue() {
71          // GIVEN
72          final IsotropicRadiationSingleCoefficient mockedIsotropicRadiationSingleCoefficient = Mockito.mock(IsotropicRadiationSingleCoefficient.class);
73          final LightFluxModel mockedFluxModel = Mockito.mock(LightFluxModel.class);
74          final RadiationPressureModel forceModel = new RadiationPressureModel(mockedFluxModel,
75                  mockedIsotropicRadiationSingleCoefficient);
76          // WHEN
77          final boolean dependsOnPositionOnly = forceModel.dependsOnPositionOnly();
78          // THEN
79          Assertions.assertTrue(dependsOnPositionOnly);
80      }
81  
82      @Test
83      void testDependsOnPositionOnlyFalse() {
84          // GIVEN
85          final BoxAndSolarArraySpacecraft mockedBoxAndSolarArraySpacecraft = Mockito.mock(BoxAndSolarArraySpacecraft.class);
86          final LightFluxModel mockedFluxModel = Mockito.mock(LightFluxModel.class);
87          final RadiationPressureModel forceModel = new RadiationPressureModel(mockedFluxModel,
88                  mockedBoxAndSolarArraySpacecraft);
89          // WHEN
90          final boolean dependsOnPositionOnly = forceModel.dependsOnPositionOnly();
91          // THEN
92          Assertions.assertFalse(dependsOnPositionOnly);
93      }
94  
95      @Test
96      void testGetEventDetectors() {
97          // GIVEN
98          final RadiationSensitive mockedRadiationSensitive = Mockito.mock(RadiationSensitive.class);
99          final LightFluxModel mockedFluxModel = Mockito.mock(LightFluxModel.class);
100         final List<EventDetector> eclipseDetectors = new ArrayList<>();
101         eclipseDetectors.add(Mockito.mock(EventDetector.class));
102         Mockito.when(mockedFluxModel.getEclipseConditionsDetector()).thenReturn(eclipseDetectors);
103         final RadiationPressureModel forceModel = new RadiationPressureModel(mockedFluxModel,
104                 mockedRadiationSensitive);
105         // WHEN
106         final Stream<?> detectors = forceModel.getEventDetectors();
107         // THEN
108         Assertions.assertEquals(eclipseDetectors.size(), detectors.toArray().length);
109     }
110 
111     @SuppressWarnings("unchecked")
112     @Test
113     void testGetFieldEventDetectors() {
114         // GIVEN
115         final RadiationSensitive mockedRadiationSensitive = Mockito.mock(RadiationSensitive.class);
116         final LightFluxModel mockedFluxModel = Mockito.mock(LightFluxModel.class);
117         final List<FieldEventDetector<Complex>> eclipseDetectors = new ArrayList<>();
118         eclipseDetectors.add(Mockito.mock(FieldEventDetector.class));
119         final ComplexField field = ComplexField.getInstance();
120         Mockito.when(mockedFluxModel.getFieldEclipseConditionsDetector(field)).thenReturn(eclipseDetectors);
121         final RadiationPressureModel forceModel = new RadiationPressureModel(mockedFluxModel,
122                 mockedRadiationSensitive);
123         // WHEN
124         final Stream<?> detectors = forceModel.getFieldEventDetectors(field);
125         // THEN
126         Assertions.assertEquals(eclipseDetectors.size(), detectors.toArray().length);
127     }
128 
129     @ParameterizedTest
130     @ValueSource(doubles = {400e3, 500e3, 800e3, 1000e3, 2000e3})
131     void testLeoPropagation(final double altitude) {
132         // GIVEN
133         Utils.setDataRoot("regular-data");
134         final IsotropicRadiationSingleCoefficient isotropicRadiationSingleCoefficient = new IsotropicRadiationSingleCoefficient(10., 1.6);
135         final ConicallyShadowedLightFluxModel lightFluxModel = new ConicallyShadowedLightFluxModel(Constants.SUN_RADIUS,
136                 CelestialBodyFactory.getSun(), Constants.EGM96_EARTH_EQUATORIAL_RADIUS);
137         final RadiationPressureModel forceModel = new RadiationPressureModel(lightFluxModel,
138                 isotropicRadiationSingleCoefficient);
139         final double radius = lightFluxModel.getOccultingBodyRadius() + altitude;
140         final NumericalPropagator propagator = createPropagator(radius);
141         propagator.addForceModel(forceModel);
142         final Orbit orbit = propagator.getInitialState().getOrbit();
143         final AbsoluteDate epoch = orbit.getDate();
144         // WHEN
145         final AbsoluteDate terminalDate = epoch.shiftedBy(orbit.getKeplerianPeriod() * 5);
146         final SpacecraftState propagateState = propagator.propagate(terminalDate);
147         // THEN
148         final SpacecraftState comparableState = computeComparableState(forceModel, radius, terminalDate);
149         final Vector3D relativePosition = comparableState.getPosition().subtract(propagateState.getPosition(comparableState.getFrame()));
150         Assertions.assertEquals(0., relativePosition.getNorm(), 1e0);
151     }
152 
153     @Test
154     void testLeoFieldPropagation() {
155         // GIVEN
156         Utils.setDataRoot("regular-data");
157         final IsotropicRadiationSingleCoefficient isotropicRadiationSingleCoefficient =
158               new IsotropicRadiationSingleCoefficient(10., 1.6);
159         final ConicallyShadowedLightFluxModel lightFluxModel = new ConicallyShadowedLightFluxModel(Constants.SUN_RADIUS,
160                                                                                                    CelestialBodyFactory.getSun(),
161                                                                                                    Constants.EGM96_EARTH_EQUATORIAL_RADIUS);
162         final RadiationPressureModel forceModel = new RadiationPressureModel(lightFluxModel,
163                                                                              isotropicRadiationSingleCoefficient);
164         final double                             radius          = lightFluxModel.getOccultingBodyRadius() + 700e3;
165         final Binary64Field                      field           = Binary64Field.getInstance();
166         final SpacecraftState                    initialState    = createState(radius);
167         final FieldNumericalPropagator<Binary64> fieldPropagator =
168               doBuildFieldPropagator(field, forceModel, initialState);
169         final AbsoluteDate                       epoch           = initialState.getDate();
170         final AbsoluteDate terminalDate =
171               epoch.shiftedBy(initialState.getOrbit().getKeplerianPeriod() * 10);
172         // WHEN
173         final FieldSpacecraftState<Binary64> propagatedState =
174               fieldPropagator.propagate(new FieldAbsoluteDate<>(field, terminalDate));
175         // THEN
176         final NumericalPropagator propagator = createPropagator(radius);
177         propagator.setOrbitType(fieldPropagator.getOrbitType());
178         propagator.setPositionAngleType(fieldPropagator.getPositionAngleType());
179         propagator.addForceModel(forceModel);
180         final SpacecraftState comparableState = propagator.propagate(terminalDate);
181         final Vector3D relativePosition = comparableState.getPosition()
182                                                          .subtract(propagatedState.getPosition(
183                                                                comparableState.getFrame()).toVector3D());
184         Assertions.assertEquals(0., relativePosition.getNorm(), 1e-3);
185     }
186 
187     private static <T extends CalculusFieldElement<T>> FieldNumericalPropagator<T>
188     doBuildFieldPropagator(final Field<T> field,
189                            final ForceModel forceModel,
190                            final SpacecraftState initialState) {
191         final FieldODEIntegratorBuilder<T> fieldIntegratoBuilder =
192               new DormandPrince54FieldIntegratorBuilder<>(1e-3, 1e2, 1e-3);
193 
194         final OrbitType propagationType = OrbitType.EQUINOCTIAL;
195         final FieldODEIntegrator<T> fieldIntegrator =
196               fieldIntegratoBuilder.buildIntegrator(field, initialState.getOrbit(), propagationType);
197 
198         final FieldNumericalPropagator<T> fieldPropagator =
199               new FieldNumericalPropagator<>(field, fieldIntegrator);
200 
201         fieldPropagator.addForceModel(forceModel);
202         fieldPropagator.setOrbitType(propagationType);
203         fieldPropagator.setInitialState(new FieldSpacecraftState<>(field, initialState));
204 
205         return fieldPropagator;
206     }
207 
208     @Test
209     void testGeoPropagation() {
210         // GIVEN
211         Utils.setDataRoot("regular-data");
212         final IsotropicRadiationSingleCoefficient isotropicRadiationSingleCoefficient = new IsotropicRadiationSingleCoefficient(10., 1.6);
213         final CylindricallyShadowedLightFluxModel lightFluxModel = new CylindricallyShadowedLightFluxModel(CelestialBodyFactory.getSun(),
214                 Constants.EGM96_EARTH_EQUATORIAL_RADIUS);
215         final RadiationPressureModel forceModel = new RadiationPressureModel(lightFluxModel,
216                 isotropicRadiationSingleCoefficient);
217         final double radius = 42000e3;
218         final NumericalPropagator propagator = createPropagator(radius);
219         propagator.addForceModel(forceModel);
220         final AbsoluteDate epoch = propagator.getInitialState().getDate();
221         // WHEN
222         final AbsoluteDate terminalDate = epoch.shiftedBy(Constants.JULIAN_DAY * 10.);
223         final SpacecraftState propagateState = propagator.propagate(terminalDate);
224         // THEN
225         final SpacecraftState comparableState = computeComparableState(forceModel, radius, terminalDate);
226         final Vector3D relativePosition = comparableState.getPosition().subtract(propagateState.getPosition(comparableState.getFrame()));
227         Assertions.assertEquals(0., relativePosition.getNorm(), 1e-6);
228     }
229 
230     private SpacecraftState computeComparableState(final RadiationPressureModel radiationPressureModel,
231                                                    final double radius, final AbsoluteDate terminalDate) {
232         final NumericalPropagator propagator = createPropagator(radius);
233         final AbstractSolarLightFluxModel lightFluxModel = (AbstractSolarLightFluxModel) radiationPressureModel.getLightFluxModel();
234         final SolarRadiationPressure solarRadiationPressure = new SolarRadiationPressure((ExtendedPositionProvider) lightFluxModel.getOccultedBody(),
235                 new OneAxisEllipsoid(lightFluxModel.getOccultingBodyRadius(), 0., FramesFactory.getGTOD(false)),
236                 radiationPressureModel.getRadiationSensitive());
237         propagator.addForceModel(solarRadiationPressure);
238         return propagator.propagate(terminalDate);
239     }
240 
241     private SpacecraftState createState(final double radius) {
242         return new SpacecraftState(createOrbit(radius), 100);
243     }
244 
245     private Orbit createOrbit(final double radius) {
246         return new EquinoctialOrbit(radius, 0., 0., 0., 0., 0., PositionAngleType.ECCENTRIC, FramesFactory.getGCRF(),
247                 AbsoluteDate.ARBITRARY_EPOCH, Constants.EGM96_EARTH_MU);
248     }
249 
250     private NumericalPropagator createPropagator(final double radius) {
251         final SpacecraftState initialState = createState(radius);
252         final Orbit initialOrbit = initialState.getOrbit();
253         final ODEIntegratorBuilder integratorBuilder = new DormandPrince54IntegratorBuilder(1e-3, 1e2, 1e-3);
254         final ODEIntegrator integrator = integratorBuilder.buildIntegrator(initialOrbit, initialOrbit.getType());
255         final NumericalPropagator propagator = new NumericalPropagator(integrator);
256         propagator.setOrbitType(initialOrbit.getType());
257         propagator.setInitialState(initialState);
258         return propagator;
259     }
260 
261 }