1   /* Contributed in the public domain.
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.propagation.conversion;
18  
19  import org.hipparchus.Field;
20  import org.hipparchus.exception.MathIllegalArgumentException;
21  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
22  import org.hipparchus.ode.FieldODEIntegrator;
23  import org.hipparchus.ode.nonstiff.AdamsBashforthFieldIntegrator;
24  import org.hipparchus.ode.nonstiff.AdamsMoultonFieldIntegrator;
25  import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaFieldIntegrator;
26  import org.hipparchus.ode.nonstiff.DormandPrince54FieldIntegrator;
27  import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
28  import org.hipparchus.ode.nonstiff.EulerFieldIntegrator;
29  import org.hipparchus.ode.nonstiff.GillFieldIntegrator;
30  import org.hipparchus.ode.nonstiff.HighamHall54FieldIntegrator;
31  import org.hipparchus.ode.nonstiff.LutherFieldIntegrator;
32  import org.hipparchus.ode.nonstiff.MidpointFieldIntegrator;
33  import org.hipparchus.ode.nonstiff.ThreeEighthesFieldIntegrator;
34  import org.hipparchus.util.Binary64;
35  import org.hipparchus.util.Binary64Field;
36  import org.junit.jupiter.api.Assertions;
37  import org.junit.jupiter.api.DisplayName;
38  import org.junit.jupiter.api.Test;
39  import org.orekit.Utils;
40  import org.orekit.bodies.CelestialBodyFactory;
41  import org.orekit.forces.ForceModel;
42  import org.orekit.forces.gravity.ThirdBodyAttraction;
43  import org.orekit.frames.Frame;
44  import org.orekit.frames.FramesFactory;
45  import org.orekit.orbits.FieldCartesianOrbit;
46  import org.orekit.orbits.FieldOrbit;
47  import org.orekit.orbits.OrbitType;
48  import org.orekit.propagation.FieldSpacecraftState;
49  import org.orekit.propagation.ToleranceProvider;
50  import org.orekit.propagation.numerical.FieldNumericalPropagator;
51  import org.orekit.time.FieldAbsoluteDate;
52  import org.orekit.utils.FieldPVCoordinates;
53  
54  public class FieldODEIntegratorTest {
55  
56      @Test
57      @DisplayName("Test FieldODEIntegrator builders")
58      void testFieldODEIntegratorBuilders() {
59  
60          // Given
61  
62          // Load Orekit data
63          Utils.setDataRoot("regular-data");
64  
65          final Field<Binary64>      field = Binary64Field.getInstance();
66          final FieldOrbit<Binary64> orbit = getReferenceOrbit();
67  
68          final double   step      = 100;
69          final Binary64 fieldStep = new Binary64(step);
70  
71          final int    nSteps  = 2;
72          final double minStep = 1;
73          final double maxStep = 300;
74          final double dP      = 1;
75  
76          final AdamsBashforthFieldIntegratorBuilder<Binary64> integratorBuilder01 =
77                  new AdamsBashforthFieldIntegratorBuilder<>(nSteps, minStep, maxStep, dP);
78  
79          final AdamsMoultonFieldIntegratorBuilder<Binary64> integratorBuilder02 =
80                  new AdamsMoultonFieldIntegratorBuilder<>(nSteps, minStep, maxStep, dP);
81  
82          final ClassicalRungeKuttaFieldIntegratorBuilder<Binary64> integratorBuilder03 =
83                  new ClassicalRungeKuttaFieldIntegratorBuilder<>(fieldStep);
84  
85          final ClassicalRungeKuttaFieldIntegratorBuilder<Binary64> integratorBuilder03Bis =
86                  new ClassicalRungeKuttaFieldIntegratorBuilder<>(step);
87  
88          final DormandPrince54FieldIntegratorBuilder<Binary64> integratorBuilder04 =
89                  new DormandPrince54FieldIntegratorBuilder<>(minStep, maxStep, dP);
90  
91          final DormandPrince853FieldIntegratorBuilder<Binary64> integratorBuilder05 =
92                  new DormandPrince853FieldIntegratorBuilder<>(minStep, maxStep, dP);
93  
94          final EulerFieldIntegratorBuilder<Binary64> integratorBuilder06 =
95                  new EulerFieldIntegratorBuilder<>(fieldStep);
96  
97          final EulerFieldIntegratorBuilder<Binary64> integratorBuilder06Bis =
98                  new EulerFieldIntegratorBuilder<>(step);
99  
100         final GillFieldIntegratorBuilder<Binary64> integratorBuilder07 =
101                 new GillFieldIntegratorBuilder<>(fieldStep);
102 
103         final GillFieldIntegratorBuilder<Binary64> integratorBuilder07Bis =
104                 new GillFieldIntegratorBuilder<>(step);
105 
106         final HighamHall54FieldIntegratorBuilder<Binary64> integratorBuilder08 =
107                 new HighamHall54FieldIntegratorBuilder<>(minStep, maxStep, dP);
108 
109         final LutherFieldIntegratorBuilder<Binary64> integratorBuilder09 =
110                 new LutherFieldIntegratorBuilder<>(fieldStep);
111 
112         final LutherFieldIntegratorBuilder<Binary64> integratorBuilder09Bis =
113                 new LutherFieldIntegratorBuilder<>(step);
114 
115         final MidpointFieldIntegratorBuilder<Binary64> integratorBuilder10 =
116                 new MidpointFieldIntegratorBuilder<>(fieldStep);
117 
118         final MidpointFieldIntegratorBuilder<Binary64> integratorBuilder10Bis =
119                 new MidpointFieldIntegratorBuilder<>(step);
120 
121         final ThreeEighthesFieldIntegratorBuilder<Binary64> integratorBuilder11 =
122                 new ThreeEighthesFieldIntegratorBuilder<>(fieldStep);
123 
124         final ThreeEighthesFieldIntegratorBuilder<Binary64> integratorBuilder11Bis =
125                 new ThreeEighthesFieldIntegratorBuilder<>(step);
126 
127         // When
128         final FieldODEIntegrator<Binary64> builtIntegrator01 =
129                 integratorBuilder01.buildIntegrator(orbit, OrbitType.CARTESIAN);
130 
131         final FieldODEIntegrator<Binary64> builtIntegrator02 =
132                 integratorBuilder02.buildIntegrator(orbit, OrbitType.CARTESIAN);
133 
134         final FieldODEIntegrator<Binary64> builtIntegrator03 =
135                 integratorBuilder03.buildIntegrator(orbit, OrbitType.CARTESIAN);
136 
137         final FieldODEIntegrator<Binary64> builtIntegrator03Bis =
138                 integratorBuilder03Bis.buildIntegrator(orbit, OrbitType.CARTESIAN);
139 
140         final FieldODEIntegrator<Binary64> builtIntegrator04 =
141                 integratorBuilder04.buildIntegrator(orbit, OrbitType.CARTESIAN);
142 
143         final FieldODEIntegrator<Binary64> builtIntegrator05 =
144                 integratorBuilder05.buildIntegrator(orbit, OrbitType.CARTESIAN);
145 
146         final FieldODEIntegrator<Binary64> builtIntegrator06 =
147                 integratorBuilder06.buildIntegrator(orbit, OrbitType.CARTESIAN);
148 
149         final FieldODEIntegrator<Binary64> builtIntegrator06Bis =
150                 integratorBuilder06Bis.buildIntegrator(orbit, OrbitType.CARTESIAN);
151 
152         final FieldODEIntegrator<Binary64> builtIntegrator07 =
153                 integratorBuilder07.buildIntegrator(orbit, OrbitType.CARTESIAN);
154 
155         final FieldODEIntegrator<Binary64> builtIntegrator07Bis =
156                 integratorBuilder07Bis.buildIntegrator(orbit, OrbitType.CARTESIAN);
157 
158         final FieldODEIntegrator<Binary64> builtIntegrator08 =
159                 integratorBuilder08.buildIntegrator(orbit, OrbitType.CARTESIAN);
160 
161         final FieldODEIntegrator<Binary64> builtIntegrator09 =
162                 integratorBuilder09.buildIntegrator(orbit, OrbitType.CARTESIAN);
163 
164         final FieldODEIntegrator<Binary64> builtIntegrator09Bis =
165                 integratorBuilder09Bis.buildIntegrator(orbit, OrbitType.CARTESIAN);
166 
167         final FieldODEIntegrator<Binary64> builtIntegrator10 =
168                 integratorBuilder10.buildIntegrator(orbit, OrbitType.CARTESIAN);
169 
170         final FieldODEIntegrator<Binary64> builtIntegrator10Bis =
171                 integratorBuilder10Bis.buildIntegrator(orbit, OrbitType.CARTESIAN);
172 
173         final FieldODEIntegrator<Binary64> builtIntegrator11 =
174                 integratorBuilder11.buildIntegrator(orbit, OrbitType.CARTESIAN);
175 
176         final FieldODEIntegrator<Binary64> builtIntegrator11Bis =
177                 integratorBuilder11Bis.buildIntegrator(orbit, OrbitType.CARTESIAN);
178 
179         // Then
180 
181         // Creating reference integrators
182         final double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(dP).getTolerances(orbit.toOrbit(), OrbitType.CARTESIAN);
183 
184         final FieldODEIntegrator<Binary64> referenceIntegrator01 =
185                 new AdamsBashforthFieldIntegrator<>(field, nSteps, minStep, maxStep, tolerances[0], tolerances[1]);
186 
187         final FieldODEIntegrator<Binary64> referenceIntegrator02 =
188                 new AdamsMoultonFieldIntegrator<>(field, nSteps, minStep, maxStep, tolerances[0], tolerances[1]);
189 
190         final FieldODEIntegrator<Binary64> referenceIntegrator03 =
191                 new ClassicalRungeKuttaFieldIntegrator<>(field, fieldStep);
192 
193         final FieldODEIntegrator<Binary64> referenceIntegrator04 =
194                 new DormandPrince54FieldIntegrator<>(field, minStep, maxStep, tolerances[0], tolerances[1]);
195 
196         final FieldODEIntegrator<Binary64> referenceIntegrator05 =
197                 new DormandPrince853FieldIntegrator<>(field, minStep, maxStep, tolerances[0], tolerances[1]);
198 
199         final FieldODEIntegrator<Binary64> referenceIntegrator06 =
200                 new EulerFieldIntegrator<>(field, fieldStep);
201 
202         final FieldODEIntegrator<Binary64> referenceIntegrator07 =
203                 new GillFieldIntegrator<>(field, fieldStep);
204 
205         final FieldODEIntegrator<Binary64> referenceIntegrator08 =
206                 new HighamHall54FieldIntegrator<>(field, minStep, maxStep, tolerances[0], tolerances[1]);
207 
208         final FieldODEIntegrator<Binary64> referenceIntegrator09 =
209                 new LutherFieldIntegrator<>(field, fieldStep);
210 
211         final FieldODEIntegrator<Binary64> referenceIntegrator10 =
212                 new MidpointFieldIntegrator<>(field, fieldStep);
213 
214         final FieldODEIntegrator<Binary64> referenceIntegrator11 =
215                 new ThreeEighthesFieldIntegrator<>(field, fieldStep);
216 
217         assertBuiltIntegrator(referenceIntegrator01, builtIntegrator01, orbit);
218         assertBuiltIntegrator(referenceIntegrator02, builtIntegrator02, orbit);
219         assertBuiltIntegrator(referenceIntegrator03, builtIntegrator03, orbit);
220         assertBuiltIntegrator(referenceIntegrator03, builtIntegrator03Bis, orbit);
221         assertBuiltIntegrator(referenceIntegrator04, builtIntegrator04, orbit);
222         assertBuiltIntegrator(referenceIntegrator05, builtIntegrator05, orbit);
223         assertBuiltIntegrator(referenceIntegrator06, builtIntegrator06, orbit);
224         assertBuiltIntegrator(referenceIntegrator06, builtIntegrator06Bis, orbit);
225         assertBuiltIntegrator(referenceIntegrator07, builtIntegrator07, orbit);
226         assertBuiltIntegrator(referenceIntegrator07, builtIntegrator07Bis, orbit);
227         assertBuiltIntegrator(referenceIntegrator08, builtIntegrator08, orbit);
228         assertBuiltIntegrator(referenceIntegrator09, builtIntegrator09, orbit);
229         assertBuiltIntegrator(referenceIntegrator09, builtIntegrator09Bis, orbit);
230         assertBuiltIntegrator(referenceIntegrator10, builtIntegrator10, orbit);
231         assertBuiltIntegrator(referenceIntegrator10, builtIntegrator10Bis, orbit);
232         assertBuiltIntegrator(referenceIntegrator11, builtIntegrator11, orbit);
233         assertBuiltIntegrator(referenceIntegrator11, builtIntegrator11Bis, orbit);
234     }
235 
236     @Test
237     @DisplayName("Test that an error is thrown if given zero as step size")
238     void testErrorThrownWhenGivenZeroStepSize() {
239         // Given
240         final double step = 0;
241 
242         // When & Then
243         Assertions.assertThrows(MathIllegalArgumentException.class, () -> new MidpointFieldIntegratorBuilder<>(step));
244     }
245 
246     private void assertBuiltIntegrator(final FieldODEIntegrator<Binary64> referenceIntegrator,
247                                        final FieldODEIntegrator<Binary64> builtIntegrator,
248                                        final FieldOrbit<Binary64> initialOrbit) {
249 
250         final Field<Binary64> field               = Binary64Field.getInstance();
251         final double          propagationDuration = 1200;
252         final ForceModel      moonAttraction      = new ThirdBodyAttraction(CelestialBodyFactory.getMoon());
253         final ForceModel      sunAttraction       = new ThirdBodyAttraction(CelestialBodyFactory.getSun());
254 
255         // Create initial state from given initial orbit
256         final FieldSpacecraftState<Binary64> initialState = new FieldSpacecraftState<>(initialOrbit);
257 
258         // Create propagator with reference integrator and initialize its state
259         final FieldNumericalPropagator<Binary64> referencePropagator =
260                 new FieldNumericalPropagator<>(field, referenceIntegrator);
261         referencePropagator.setInitialState(initialState);
262         referencePropagator.addForceModel(sunAttraction);
263         referencePropagator.addForceModel(moonAttraction);
264 
265         // Create propagator with built integrator and initialize its state
266         final FieldNumericalPropagator<Binary64> builtPropagator =
267                 new FieldNumericalPropagator<>(field, builtIntegrator);
268         builtPropagator.setInitialState(initialState);
269         builtPropagator.addForceModel(sunAttraction);
270         builtPropagator.addForceModel(moonAttraction);
271 
272         // Propagate
273         final FieldSpacecraftState<Binary64> referencePropagatedState =
274                 referencePropagator.propagate(initialOrbit.getDate().shiftedBy(propagationDuration));
275         final FieldOrbit<Binary64> referencePropagatedOrbit = referencePropagatedState.getOrbit();
276 
277         final FieldSpacecraftState<Binary64> builtPropagatedState =
278                 builtPropagator.propagate(initialOrbit.getDate().shiftedBy(propagationDuration));
279         final FieldOrbit<Binary64> builtPropagatedOrbit = builtPropagatedState.getOrbit();
280 
281         // Assert that results given with reference and built integrators are the same
282         final double assertionTolerance = 1e-15;
283         Assertions.assertEquals(referencePropagatedOrbit.getA().getReal(), builtPropagatedOrbit.getA().getReal(),
284                                 assertionTolerance);
285         Assertions.assertEquals(referencePropagatedOrbit.getEquinoctialEx().getReal(),
286                                 builtPropagatedOrbit.getEquinoctialEx().getReal(), assertionTolerance);
287         Assertions.assertEquals(referencePropagatedOrbit.getEquinoctialEy().getReal(),
288                                 builtPropagatedOrbit.getEquinoctialEy().getReal(), assertionTolerance);
289         Assertions.assertEquals(referencePropagatedOrbit.getHx().getReal(), builtPropagatedOrbit.getHx().getReal(),
290                                 assertionTolerance);
291         Assertions.assertEquals(referencePropagatedOrbit.getHy().getReal(), builtPropagatedOrbit.getHy().getReal(),
292                                 assertionTolerance);
293         Assertions.assertEquals(referencePropagatedOrbit.getLM().getReal(), builtPropagatedOrbit.getLM().getReal(),
294                                 assertionTolerance);
295     }
296 
297     private FieldOrbit<Binary64> getReferenceOrbit() {
298 
299         final Field<Binary64> field = Binary64Field.getInstance();
300 
301         final FieldAbsoluteDate<Binary64> date  = new FieldAbsoluteDate<>(field);
302         final Frame                       frame = FramesFactory.getGCRF();
303         final Binary64                    mu    = new Binary64(398600e9);
304 
305         final FieldVector3D<Binary64> position = new FieldVector3D<>(new Binary64(6378000 + 400000),
306                                                                      new Binary64(0),
307                                                                      new Binary64(0));
308         final FieldVector3D<Binary64> velocity = new FieldVector3D<>(new Binary64(5500),
309                                                                      new Binary64(5500),
310                                                                      new Binary64(0));
311         final FieldPVCoordinates<Binary64> pv = new FieldPVCoordinates<>(position, velocity);
312 
313         return new FieldCartesianOrbit<>(pv, frame, date, mu);
314     }
315 }