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