1   /* Copyright 2002-2025 CS GROUP
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.forces.maneuvers;
18  
19  import java.util.List;
20  
21  import org.hamcrest.MatcherAssert;
22  import org.hipparchus.Field;
23  import org.hipparchus.analysis.differentiation.DSFactory;
24  import org.hipparchus.analysis.differentiation.DerivativeStructure;
25  import org.hipparchus.analysis.differentiation.Gradient;
26  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
27  import org.hipparchus.geometry.euclidean.threed.Rotation;
28  import org.hipparchus.geometry.euclidean.threed.Vector3D;
29  import org.hipparchus.ode.ODEIntegrator;
30  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
31  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
32  import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
33  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
34  import org.hipparchus.util.FastMath;
35  import org.hipparchus.util.MathUtils;
36  import org.junit.jupiter.api.Assertions;
37  import org.junit.jupiter.api.BeforeEach;
38  import org.junit.jupiter.api.Test;
39  import org.orekit.OrekitMatchers;
40  import org.orekit.Utils;
41  import org.orekit.attitudes.AttitudeProvider;
42  import org.orekit.attitudes.FrameAlignedProvider;
43  import org.orekit.attitudes.LofOffset;
44  import org.orekit.forces.AbstractLegacyForceModelTest;
45  import org.orekit.forces.ForceModel;
46  import org.orekit.frames.Frame;
47  import org.orekit.frames.FramesFactory;
48  import org.orekit.frames.LOFType;
49  import org.orekit.orbits.CartesianOrbit;
50  import org.orekit.orbits.CircularOrbit;
51  import org.orekit.orbits.FieldKeplerianOrbit;
52  import org.orekit.orbits.KeplerianOrbit;
53  import org.orekit.orbits.Orbit;
54  import org.orekit.orbits.OrbitType;
55  import org.orekit.orbits.PositionAngleType;
56  import org.orekit.propagation.FieldSpacecraftState;
57  import org.orekit.propagation.SpacecraftState;
58  import org.orekit.propagation.ToleranceProvider;
59  import org.orekit.propagation.events.EventDetector;
60  import org.orekit.propagation.numerical.FieldNumericalPropagator;
61  import org.orekit.propagation.numerical.NumericalPropagator;
62  import org.orekit.time.AbsoluteDate;
63  import org.orekit.time.DateComponents;
64  import org.orekit.time.FieldAbsoluteDate;
65  import org.orekit.time.TimeComponents;
66  import org.orekit.time.TimeScalesFactory;
67  import org.orekit.utils.Constants;
68  import org.orekit.utils.FieldPVCoordinates;
69  import org.orekit.utils.PVCoordinates;
70  import org.orekit.utils.ParameterDriver;
71  import org.orekit.utils.TimeStampedPVCoordinates;
72  
73  class ConstantThrustManeuverTest extends AbstractLegacyForceModelTest {
74  
75      // Body mu
76      private double mu;
77  
78      @Override
79      protected FieldVector3D<DerivativeStructure> accelerationDerivatives(final ForceModel forceModel,
80                                                                           final FieldSpacecraftState<DerivativeStructure> state) {
81          try {
82              final boolean firing = ((ConstantThrustManeuver) forceModel).isFiring(state);
83              
84              final Vector3D thrustVector = ((ConstantThrustManeuver) forceModel).getThrustVector();
85              final double thrust = thrustVector.getNorm();
86              final Vector3D direction = thrustVector.normalize();
87  
88              if (firing) {
89                  return new FieldVector3D<>(state.getMass().reciprocal().multiply(thrust),
90                                             state.getAttitude().getRotation().applyInverseTo(direction));
91              } else {
92                  // constant (and null) acceleration when not firing
93                  return FieldVector3D.getZero(state.getMass().getField());
94              }
95          } catch (IllegalArgumentException | SecurityException e) {
96              Assertions.fail(e.getLocalizedMessage());
97              return null;
98          }
99      }
100 
101     @Override
102     protected FieldVector3D<Gradient> accelerationDerivativesGradient(final ForceModel forceModel,
103                                                                       final FieldSpacecraftState<Gradient> state) {
104         try {
105             final boolean firing = ((ConstantThrustManeuver) forceModel).isFiring(state);
106             
107             final Vector3D thrustVector = ((ConstantThrustManeuver) forceModel).getThrustVector();
108             final double thrust = thrustVector.getNorm();
109             final Vector3D direction = thrustVector.normalize();
110 
111             if (firing) {
112                 return new FieldVector3D<>(state.getMass().reciprocal().multiply(thrust),
113                                            state.getAttitude().getRotation().applyInverseTo(direction));
114             } else {
115                 // constant (and null) acceleration when not firing
116                 return FieldVector3D.getZero(state.getMass().getField());
117             }
118         } catch (IllegalArgumentException | SecurityException e) {
119             Assertions.fail(e.getLocalizedMessage());
120             return null;
121         }
122     }
123 
124     private CircularOrbit dummyOrbit(AbsoluteDate date) {
125         return new CircularOrbit(new PVCoordinates(Vector3D.PLUS_I, Vector3D.PLUS_J),
126                                  FramesFactory.getEME2000(), date, mu);
127     }
128 
129     @Test
130     void testJacobianVs80Implementation() {
131         final double isp = 318;
132         final double mass = 2500;
133         final double a = 24396159;
134         final double e = 0.72831215;
135         final double i = FastMath.toRadians(7);
136         final double omega = FastMath.toRadians(180);
137         final double OMEGA = FastMath.toRadians(261);
138         final double lv = 0;
139 
140         final double duration = 3653.99;
141         final double f = 420;
142 
143         final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 01, 01),
144                                                        new TimeComponents(23, 30, 00.000),
145                                                        TimeScalesFactory.getUTC());
146         final Orbit orbit =
147             new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
148                                FramesFactory.getEME2000(), initDate, mu);
149         final AttitudeProvider law = new LofOffset(orbit.getFrame(), LOFType.LVLH_CCSDS);
150        final SpacecraftState initialState =
151             new SpacecraftState(orbit, law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
152 
153         final AbsoluteDate fireDate = new AbsoluteDate(new DateComponents(2004, 01, 02),
154                                                        new TimeComponents(04, 15, 34.080),
155                                                        TimeScalesFactory.getUTC());
156         final ConstantThrustManeuver maneuver =
157             new ConstantThrustManeuver(fireDate, duration, f, isp, Vector3D.PLUS_I);
158 
159         // before maneuver (Jacobian wrt. state is zero)
160         checkStateJacobianVs80Implementation(initialState, maneuver, law, 1.0e-50, false);
161 
162         // in maneuver
163         SpacecraftState startState = initialState.shiftedBy(fireDate.durationFrom(initDate));
164         EventDetector d = maneuver.getEventDetectors().findFirst().get();
165         d.getHandler().eventOccurred(startState, d, true);
166         SpacecraftState midState = startState.shiftedBy(duration / 2.0);
167         checkStateJacobianVs80Implementation(midState, maneuver, law, 1.0e-20, false);
168 
169     }
170 
171     @Test
172     void testJacobianVs80ImplementationGradient() {
173         final double isp = 318;
174         final double mass = 2500;
175         final double a = 24396159;
176         final double e = 0.72831215;
177         final double i = FastMath.toRadians(7);
178         final double omega = FastMath.toRadians(180);
179         final double OMEGA = FastMath.toRadians(261);
180         final double lv = 0;
181 
182         final double duration = 3653.99;
183         final double f = 420;
184 
185         final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 01, 01),
186                                                        new TimeComponents(23, 30, 00.000),
187                                                        TimeScalesFactory.getUTC());
188         final Orbit orbit =
189             new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
190                                FramesFactory.getEME2000(), initDate, mu);
191         final AttitudeProvider law = new LofOffset(orbit.getFrame(), LOFType.LVLH_CCSDS);
192        final SpacecraftState initialState =
193             new SpacecraftState(orbit, law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
194 
195         final AbsoluteDate fireDate = new AbsoluteDate(new DateComponents(2004, 01, 02),
196                                                        new TimeComponents(04, 15, 34.080),
197                                                        TimeScalesFactory.getUTC());
198         final ConstantThrustManeuver maneuver =
199             new ConstantThrustManeuver(fireDate, duration, f, isp, Vector3D.PLUS_I);
200 
201         // before maneuver (Jacobian wrt. state is zero)
202         checkStateJacobianVs80ImplementationGradient(initialState, maneuver, law, 1.0e-50, false);
203 
204         // in maneuver
205         SpacecraftState startState = initialState.shiftedBy(fireDate.durationFrom(initDate));
206         EventDetector d = maneuver.getEventDetectors().findFirst().get();
207         d.getHandler().eventOccurred(startState, d, true);
208         SpacecraftState midState = startState.shiftedBy(duration / 2.0);
209         checkStateJacobianVs80ImplementationGradient(midState, maneuver, law, 1.0e-20, false);
210 
211     }
212 
213     @Test
214     void testPositiveDuration() {
215         AbsoluteDate date = new AbsoluteDate(new DateComponents(2004, 01, 01),
216                                              new TimeComponents(23, 30, 00.000),
217                                              TimeScalesFactory.getUTC());
218         ConstantThrustManeuver maneuver =
219             new ConstantThrustManeuver(date, 10.0, 400.0, 300.0, Vector3D.PLUS_K);
220         Assertions.assertFalse(maneuver.dependsOnPositionOnly());
221         Assertions.assertNull(maneuver.getAttitudeOverride());
222         Assertions.assertEquals(0.0, Vector3D.distance(maneuver.getDirection(), Vector3D.PLUS_K), 1.0e-15);
223         Assertions.assertEquals(10.0, maneuver.getDuration(), 1.0e-15);
224         Assertions.assertEquals(0.0, date.durationFrom(maneuver.getStartDate()), 1.0e-15);
225         Assertions.assertEquals(0.0, date.shiftedBy(10.0).durationFrom(maneuver.getEndDate()), 1.0e-15);
226         Assertions.assertEquals("", maneuver.getName());
227         List<ParameterDriver> drivers = maneuver.getParametersDrivers();
228         Assertions.assertEquals(6, drivers.size());
229         Assertions.assertEquals("thrust", drivers.get(0).getName());
230         Assertions.assertEquals("flow rate", drivers.get(1).getName());
231         EventDetector detector = maneuver.getEventDetectors().findFirst().get();
232 
233         Orbit o1 = dummyOrbit(date.shiftedBy(- 1.0));
234         Assertions.assertTrue(detector.g(new SpacecraftState(o1)) < 0);
235         Orbit o2 = dummyOrbit(date.shiftedBy(  1.0));
236         Assertions.assertTrue(detector.g(new SpacecraftState(o2)) > 0);
237         Orbit o3 = dummyOrbit(date.shiftedBy(  9.0));
238         Assertions.assertTrue(detector.g(new SpacecraftState(o3)) > 0);
239         Orbit o4 = dummyOrbit(date.shiftedBy( 11.0));
240         Assertions.assertTrue(detector.g(new SpacecraftState(o4)) < 0);
241     }
242 
243     @Test
244     void testNegativeDuration() {
245         AbsoluteDate date = new AbsoluteDate(new DateComponents(2004, 01, 01),
246                                              new TimeComponents(23, 30, 00.000),
247                                              TimeScalesFactory.getUTC());
248         ConstantThrustManeuver maneuver =
249             new ConstantThrustManeuver(date, -10.0, 400.0, 300.0, Vector3D.PLUS_K,
250                                        "1A-");
251         List<ParameterDriver> drivers = maneuver.getParametersDrivers();
252         Assertions.assertEquals(6, drivers.size());
253         Assertions.assertEquals("1A-thrust", drivers.get(0).getName());
254         Assertions.assertEquals("1A-flow rate", drivers.get(1).getName());
255         EventDetector detector = maneuver.getEventDetectors().findFirst().get();
256 
257         Orbit o1 = dummyOrbit(date.shiftedBy(-11.0));
258         Assertions.assertTrue(detector.g(new SpacecraftState(o1)) < 0);
259         Orbit o2 = dummyOrbit(date.shiftedBy( -9.0));
260         Assertions.assertTrue(detector.g(new SpacecraftState(o2)) > 0);
261         Orbit o3 = dummyOrbit(date.shiftedBy( -1.0));
262         Assertions.assertTrue(detector.g(new SpacecraftState(o3)) > 0);
263         Orbit o4 = dummyOrbit(date.shiftedBy(  1.0));
264         Assertions.assertTrue(detector.g(new SpacecraftState(o4)) < 0);
265     }
266 
267     @Test
268     void testRoughBehaviour() {
269         final double isp = 318;
270         final double mass = 2500;
271         final double a = 24396159;
272         final double e = 0.72831215;
273         final double i = FastMath.toRadians(7);
274         final double omega = FastMath.toRadians(180);
275         final double OMEGA = FastMath.toRadians(261);
276         final double lv = 0;
277 
278         final double duration = 3653.99;
279         final double f = 420;
280         final double delta = FastMath.toRadians(-7.4978);
281         final double alpha = FastMath.toRadians(351);
282         final AttitudeProvider law = new FrameAlignedProvider(new Rotation(new Vector3D(alpha, delta),
283                                                                            Vector3D.PLUS_I));
284 
285         final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 01, 01),
286                                                        new TimeComponents(23, 30, 00.000),
287                                                        TimeScalesFactory.getUTC());
288         final Orbit orbit =
289             new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
290                                FramesFactory.getEME2000(), initDate, mu);
291         final SpacecraftState initialState =
292             new SpacecraftState(orbit, law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
293 
294         final AbsoluteDate fireDate = new AbsoluteDate(new DateComponents(2004, 01, 02),
295                                                        new TimeComponents(04, 15, 34.080),
296                                                        TimeScalesFactory.getUTC());
297         final ConstantThrustManeuver maneuver =
298             new ConstantThrustManeuver(fireDate, duration, f, isp, Vector3D.PLUS_I);
299         Assertions.assertEquals(f,   maneuver.getThrustMagnitude(), 1.0e-10);
300         Assertions.assertEquals(isp, maneuver.getIsp(),    1.0e-10);
301 
302         double[] absTolerance = {
303             0.001, 1.0e-9, 1.0e-9, 1.0e-6, 1.0e-6, 1.0e-6, 0.001
304         };
305         double[] relTolerance = {
306             1.0e-7, 1.0e-4, 1.0e-4, 1.0e-7, 1.0e-7, 1.0e-7, 1.0e-7
307         };
308         AdaptiveStepsizeIntegrator integrator =
309             new DormandPrince853Integrator(0.001, 1000, absTolerance, relTolerance);
310         integrator.setInitialStepSize(60);
311         final NumericalPropagator propagator = new NumericalPropagator(integrator);
312         propagator.setInitialState(initialState);
313         propagator.setAttitudeProvider(law);
314         propagator.addForceModel(maneuver);
315         final SpacecraftState finalorb = propagator.propagate(fireDate.shiftedBy(3800));
316 
317         final double massTolerance =
318                 FastMath.abs(maneuver.getFlowRate()) * maneuver.getEventDetectors().findFirst().get().getThreshold();
319         Assertions.assertEquals(2007.8824544261233, finalorb.getMass(), massTolerance);
320         Assertions.assertEquals(2.6872, FastMath.toDegrees(MathUtils.normalizeAngle(finalorb.getOrbit().getI(), FastMath.PI)), 1e-4);
321         Assertions.assertEquals(28970, finalorb.getOrbit().getA()/1000, 1);
322 
323     }
324 
325     /**Testing if the propagation between the FieldPropagation and the propagation
326      * is equivalent.
327      * Also testing if propagating X+dX with the propagation is equivalent to
328      * propagation X with the FieldPropagation and then applying the taylor
329      * expansion of dX to the result.*/
330     @Test
331     void testRealField() {
332         DSFactory factory = new DSFactory(6, 5);
333         DerivativeStructure a_0 = factory.variable(0, 7e7);
334         DerivativeStructure e_0 = factory.variable(1, 0.4);
335         DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
336         DerivativeStructure R_0 = factory.variable(3, 0.7);
337         DerivativeStructure O_0 = factory.variable(4, 0.5);
338         DerivativeStructure n_0 = factory.variable(5, 0.1);
339 
340         Field<DerivativeStructure> field = a_0.getField();
341         DerivativeStructure zero = field.getZero();
342 
343         FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
344 
345         Frame EME = FramesFactory.getEME2000();
346 
347         FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
348                                                                                  PositionAngleType.MEAN,
349                                                                                  EME,
350                                                                                  J2000,
351                                                                                  zero.add(Constants.EIGEN5C_EARTH_MU));
352 
353         FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
354 
355         SpacecraftState iSR = initialState.toSpacecraftState();
356 
357         final OrbitType type = OrbitType.KEPLERIAN;
358         double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(FKO.toOrbit(), type);
359 
360 
361         AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator =
362                         new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
363         integrator.setInitialStepSize(60);
364         AdaptiveStepsizeIntegrator RIntegrator =
365                         new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
366         RIntegrator.setInitialStepSize(60);
367 
368         FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
369         FNP.setOrbitType(type);
370         FNP.setInitialState(initialState);
371 
372         NumericalPropagator NP = new NumericalPropagator(RIntegrator);
373         NP.setOrbitType(type);
374         NP.setInitialState(iSR);
375 
376         final ConstantThrustManeuver forceModel =  new ConstantThrustManeuver(J2000.toAbsoluteDate().shiftedBy(100), 100.0, 400.0, 300.0, Vector3D.PLUS_K);
377 
378         FNP.addForceModel(forceModel);
379         NP.addForceModel(forceModel);
380 
381         // Do the test
382         checkRealFieldPropagation(FKO, PositionAngleType.MEAN, 1005., NP, FNP,
383                                   1.0e-15, 5.0e-10, 3.0e-11, 3.0e-10,
384                                   1, false);
385     }
386 
387     /**Testing if the propagation between the FieldPropagation and the propagation
388      * is equivalent.
389      * Also testing if propagating X+dX with the propagation is equivalent to
390      * propagation X with the FieldPropagation and then applying the taylor
391      * expansion of dX to the result.*/
392     @Test
393     void testRealFieldGradient() {
394         int freeParameters = 6;
395         Gradient a_0 = Gradient.variable(freeParameters, 0, 7e7);
396         Gradient e_0 = Gradient.variable(freeParameters, 1, 0.4);
397         Gradient i_0 = Gradient.variable(freeParameters, 2, 85 * FastMath.PI / 180);
398         Gradient R_0 = Gradient.variable(freeParameters, 3, 0.7);
399         Gradient O_0 = Gradient.variable(freeParameters, 4, 0.5);
400         Gradient n_0 = Gradient.variable(freeParameters, 5, 0.1);
401 
402         Field<Gradient> field = a_0.getField();
403         Gradient zero = field.getZero();
404 
405         FieldAbsoluteDate<Gradient> J2000 = new FieldAbsoluteDate<>(field);
406 
407         Frame EME = FramesFactory.getEME2000();
408 
409         FieldKeplerianOrbit<Gradient> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
410                                                                       PositionAngleType.MEAN,
411                                                                       EME,
412                                                                       J2000,
413                                                                       zero.add(Constants.EIGEN5C_EARTH_MU));
414 
415         FieldSpacecraftState<Gradient> initialState = new FieldSpacecraftState<>(FKO);
416 
417         SpacecraftState iSR = initialState.toSpacecraftState();
418 
419         final OrbitType type = OrbitType.KEPLERIAN;
420         double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(FKO.toOrbit(), type);
421 
422 
423         AdaptiveStepsizeFieldIntegrator<Gradient> integrator =
424                         new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
425         integrator.setInitialStepSize(60);
426         AdaptiveStepsizeIntegrator RIntegrator =
427                         new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
428         RIntegrator.setInitialStepSize(60);
429 
430         FieldNumericalPropagator<Gradient> FNP = new FieldNumericalPropagator<>(field, integrator);
431         FNP.setOrbitType(type);
432         FNP.setInitialState(initialState);
433 
434         NumericalPropagator NP = new NumericalPropagator(RIntegrator);
435         NP.setOrbitType(type);
436         NP.setInitialState(iSR);
437 
438         final ConstantThrustManeuver forceModel =  new ConstantThrustManeuver(J2000.toAbsoluteDate().shiftedBy(100), 100.0, 400.0, 300.0, Vector3D.PLUS_K);
439 
440         FNP.addForceModel(forceModel);
441         NP.addForceModel(forceModel);
442 
443         // Do the test
444         checkRealFieldPropagationGradient(FKO, PositionAngleType.MEAN, 1005., NP, FNP,
445                                   1.0e-15, 1.3e-02, 2.9e-04, 2.4e-3,
446                                   1, false);
447     }
448 
449     /**Same test as the previous one but not adding the ForceModel to the NumericalPropagator
450     it is a test to validate the previous test.
451     (to test if the ForceModel it's actually
452     doing something in the Propagator and the FieldPropagator)*/
453     @Test
454     public void RealFieldExpectErrorTest() {
455         DSFactory factory = new DSFactory(6, 0);
456         DerivativeStructure a_0 = factory.variable(0, 7e7);
457         DerivativeStructure e_0 = factory.variable(1, 0.4);
458         DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
459         DerivativeStructure R_0 = factory.variable(3, 0.7);
460         DerivativeStructure O_0 = factory.variable(4, 0.5);
461         DerivativeStructure n_0 = factory.variable(5, 0.1);
462 
463         Field<DerivativeStructure> field = a_0.getField();
464         DerivativeStructure zero = field.getZero();
465 
466         FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
467 
468         Frame EME = FramesFactory.getEME2000();
469 
470         FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
471                                                                                  PositionAngleType.MEAN,
472                                                                                  EME,
473                                                                                  J2000,
474                                                                                  zero.add(Constants.EIGEN5C_EARTH_MU));
475 
476         FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
477 
478         SpacecraftState iSR = initialState.toSpacecraftState();
479 
480         final OrbitType type = OrbitType.KEPLERIAN;
481         double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(FKO.toOrbit(), type);
482 
483 
484         AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator =
485                         new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
486         integrator.setInitialStepSize(60);
487         AdaptiveStepsizeIntegrator RIntegrator =
488                         new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
489         RIntegrator.setInitialStepSize(60);
490 
491         FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
492         FNP.setOrbitType(type);
493         FNP.setInitialState(initialState);
494 
495         NumericalPropagator NP = new NumericalPropagator(RIntegrator);
496         NP.setInitialState(iSR);
497 
498         final ConstantThrustManeuver forceModel =  new ConstantThrustManeuver(J2000.toAbsoluteDate().shiftedBy(100), 100.0, 400.0, 300.0, Vector3D.PLUS_K);
499 
500         FNP.addForceModel(forceModel);
501      //NOT ADDING THE FORCE MODEL TO THE NUMERICAL PROPAGATOR   NP.addForceModel(forceModel);
502 
503         FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(1000.);
504         FieldSpacecraftState<DerivativeStructure> finalState_DS = FNP.propagate(target);
505         SpacecraftState finalState_R = NP.propagate(target.toAbsoluteDate());
506         FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
507         PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
508 
509         Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getX() - finPVC_R.getPosition().getX()) < FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
510         Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getY() - finPVC_R.getPosition().getY()) < FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
511         Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getZ() - finPVC_R.getPosition().getZ()) < FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
512     }
513 
514     @Test
515     void testForwardAndBackward() {
516         final double isp = 318;
517         final double mass = 2500;
518         final double a = 24396159;
519         final double e = 0.72831215;
520         final double i = FastMath.toRadians(7);
521         final double omega = FastMath.toRadians(180);
522         final double OMEGA = FastMath.toRadians(261);
523         final double lv = 0;
524 
525         final double duration = 3653.99;
526         final double f = 420;
527         final double delta = FastMath.toRadians(-7.4978);
528         final double alpha = FastMath.toRadians(351);
529         final AttitudeProvider law = new FrameAlignedProvider(new Rotation(new Vector3D(alpha, delta),
530                                                                            Vector3D.PLUS_I));
531 
532         final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 01, 01),
533                                                        new TimeComponents(23, 30, 00.000),
534                                                        TimeScalesFactory.getUTC());
535         final Orbit orbit =
536             new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
537                                FramesFactory.getEME2000(), initDate, mu);
538         final SpacecraftState initialState =
539             new SpacecraftState(orbit, law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
540 
541         final AbsoluteDate fireDate = new AbsoluteDate(new DateComponents(2004, 01, 02),
542                                                        new TimeComponents(04, 15, 34.080),
543                                                        TimeScalesFactory.getUTC());
544         final ConstantThrustManeuver maneuver =
545             new ConstantThrustManeuver(fireDate, duration, f, isp, Vector3D.PLUS_I);
546         Assertions.assertEquals(f,   maneuver.getThrustMagnitude(), 1.0e-10);
547         Assertions.assertEquals(isp, maneuver.getIsp(),    1.0e-10);
548 
549         final OrbitType orbitType = OrbitType.KEPLERIAN;
550         double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1e-5).getTolerances(orbit, orbitType);
551         AdaptiveStepsizeIntegrator integrator1 =
552             new DormandPrince853Integrator(1.0e-5, 1000, tol[0], tol[1]);
553         integrator1.setInitialStepSize(60);
554         final NumericalPropagator propagator1 = new NumericalPropagator(integrator1);
555         propagator1.setOrbitType(orbitType);
556         propagator1.setPositionAngleType(PositionAngleType.TRUE);
557         propagator1.setInitialState(initialState);
558         propagator1.setAttitudeProvider(law);
559         propagator1.addForceModel(maneuver);
560         final SpacecraftState finalState = propagator1.propagate(fireDate.shiftedBy(3800));
561         Assertions.assertFalse(maneuver.isFiring(fireDate.shiftedBy(-1.0e-6)));
562         Assertions.assertTrue(maneuver.isFiring(fireDate.shiftedBy(+1.0e-6)));
563         Assertions.assertTrue(maneuver.isFiring(fireDate.shiftedBy(0.5 * duration)));
564         Assertions.assertTrue(maneuver.isFiring(fireDate.shiftedBy(duration - 1.0e-6)));
565         Assertions.assertFalse(maneuver.isFiring(fireDate.shiftedBy(duration + 1.0e-6)));
566 
567         AdaptiveStepsizeIntegrator integrator2 =
568                         new DormandPrince853Integrator(1.0e-5, 1000, tol[0], tol[1]);
569         integrator2.setInitialStepSize(60);
570         final NumericalPropagator propagator2 = new NumericalPropagator(integrator2);
571         propagator2.setOrbitType(orbitType);
572         propagator2.setPositionAngleType(propagator1.getPositionAngleType());
573         propagator2.setInitialState(finalState);
574         propagator2.setAttitudeProvider(law);
575         propagator2.addForceModel(maneuver);
576         final SpacecraftState recoveredState = propagator2.propagate(orbit.getDate());
577         final Vector3D refPosition = initialState.getPosition();
578         final Vector3D recoveredPosition = recoveredState.getPosition();
579         Assertions.assertEquals(0.0, Vector3D.distance(refPosition, recoveredPosition), 1.1e-3);
580         Assertions.assertEquals(initialState.getMass(), recoveredState.getMass(), 1.4e-8);
581         Assertions.assertFalse(maneuver.isFiring(fireDate.shiftedBy(-1.0e-6)));
582         Assertions.assertTrue(maneuver.isFiring(fireDate.shiftedBy(+1.0e-6)));
583         Assertions.assertTrue(maneuver.isFiring(fireDate.shiftedBy(0.5 * duration)));
584         Assertions.assertTrue(maneuver.isFiring(fireDate.shiftedBy(duration - 1.0e-6)));
585         Assertions.assertFalse(maneuver.isFiring(fireDate.shiftedBy(duration + 1.0e-6)));
586 
587     }
588 
589     @Test
590     void testParameterDerivative() {
591 
592         // pos-vel (from a ZOOM ephemeris reference)
593         final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
594         final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
595         final SpacecraftState state =
596                 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
597                                                        FramesFactory.getGCRF(),
598                                                        new AbsoluteDate(2005, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
599                                                        Constants.EIGEN5C_EARTH_MU));
600 
601         final ConstantThrustManeuver maneuver =
602                 new ConstantThrustManeuver(state.getDate().shiftedBy(-10), 100.0, 20.0, 350.0,
603                                            Vector3D.PLUS_I, "along-X-");
604         maneuver.init(state, state.getDate().shiftedBy(3600.0));
605 
606         checkParameterDerivative(state, maneuver, "along-X-thrust",    1.0e-3, 7.1e-14);
607         checkParameterDerivative(state, maneuver, "along-X-flow rate", 1.0e-3, 1.0e-15);
608 
609     }
610 
611     @Test
612     void testParameterDerivativeGradient() {
613 
614         // pos-vel (from a ZOOM ephemeris reference)
615         final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
616         final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
617         final SpacecraftState state =
618                 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
619                                                        FramesFactory.getGCRF(),
620                                                        new AbsoluteDate(2005, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
621                                                        Constants.EIGEN5C_EARTH_MU));
622 
623         final ConstantThrustManeuver maneuver =
624                 new ConstantThrustManeuver(state.getDate().shiftedBy(-10), 100.0, 20.0, 350.0,
625                                            Vector3D.PLUS_I, "along-X-");
626         maneuver.init(state, state.getDate().shiftedBy(3600.0));
627 
628         checkParameterDerivativeGradient(state, maneuver, "along-X-thrust",    1.0e-3, 3.0e-13);
629         checkParameterDerivativeGradient(state, maneuver, "along-X-flow rate", 1.0e-3, 1.0e-15);
630 
631     }
632 
633     @Test
634     void testInertialManeuver() {
635         final double isp = 318;
636         final double mass = 2500;
637         final double a = 24396159;
638         final double e = 0.72831215;
639         final double i = FastMath.toRadians(7);
640         final double omega = FastMath.toRadians(180);
641         final double OMEGA = FastMath.toRadians(261);
642         final double lv = 0;
643 
644         final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 01, 01),
645                                                        new TimeComponents(23, 30, 00.000),
646                                                        TimeScalesFactory.getUTC());
647         final Orbit orbit =
648             new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
649                                FramesFactory.getEME2000(), initDate, mu);
650 
651         final double duration = 3653.99;
652         final double f = 420;
653         final double delta = FastMath.toRadians(-7.4978);
654         final double alpha = FastMath.toRadians(351);
655         final AttitudeProvider inertialLaw = new FrameAlignedProvider(new Rotation(new Vector3D(alpha, delta),
656                                                                                    Vector3D.PLUS_I));
657         final AttitudeProvider lofLaw = new LofOffset(orbit.getFrame(), LOFType.VNC);
658 
659         final SpacecraftState initialState =
660             new SpacecraftState(orbit, inertialLaw.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
661 
662         final AbsoluteDate fireDate = new AbsoluteDate(new DateComponents(2004, 01, 02),
663                                                        new TimeComponents(04, 15, 34.080),
664                                                        TimeScalesFactory.getUTC());
665         final ConstantThrustManeuver maneuverWithoutOverride =
666             new ConstantThrustManeuver(fireDate, duration, f, isp, Vector3D.PLUS_I);
667         Assertions.assertEquals(f,   maneuverWithoutOverride.getThrustMagnitude(), 1.0e-10);
668         Assertions.assertEquals(isp, maneuverWithoutOverride.getIsp(),    1.0e-10);
669 
670         // reference propagation:
671         // propagator already uses inertial law
672         // maneuver does not need to override it to get an inertial maneuver
673         double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1.).getTolerances(orbit, OrbitType.KEPLERIAN);
674         AdaptiveStepsizeIntegrator integrator1 =
675             new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
676         integrator1.setInitialStepSize(60);
677         final NumericalPropagator propagator1 = new NumericalPropagator(integrator1);
678         propagator1.setInitialState(initialState);
679         propagator1.setAttitudeProvider(inertialLaw);
680         propagator1.addForceModel(maneuverWithoutOverride);
681         final SpacecraftState finalState1 = propagator1.propagate(fireDate.shiftedBy(3800));
682 
683 
684         // test propagation:
685         // propagator uses a LOF-aligned law
686         // maneuver needs to override it to get an inertial maneuver
687         final ConstantThrustManeuver maneuverWithOverride =
688                         new ConstantThrustManeuver(fireDate, duration, f, isp,
689                                                    inertialLaw, Vector3D.PLUS_I);
690         Assertions.assertEquals(f,   maneuverWithoutOverride.getThrustMagnitude(), 1.0e-10);
691         Assertions.assertEquals(isp, maneuverWithoutOverride.getIsp(),    1.0e-10);
692 
693         AdaptiveStepsizeIntegrator integrator2 =
694                         new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
695         integrator2.setInitialStepSize(60);
696         final NumericalPropagator propagator2 = new NumericalPropagator(integrator2);
697         propagator2.setInitialState(initialState);
698         propagator2.setAttitudeProvider(lofLaw);
699         propagator2.addForceModel(maneuverWithOverride);
700         final SpacecraftState finalState2 = propagator2.propagate(finalState1.getDate());
701         MatcherAssert.assertThat(finalState2.getPVCoordinates(),
702                                  OrekitMatchers.pvCloseTo(finalState1.getPVCoordinates(),
703                                                           1.0e-10));
704         Assertions.assertFalse(maneuverWithoutOverride.isFiring(fireDate.shiftedBy(-1.0e-6)));
705         Assertions.assertTrue(maneuverWithoutOverride.isFiring(fireDate.shiftedBy(+1.0e-6)));
706         Assertions.assertTrue(maneuverWithoutOverride.isFiring(fireDate.shiftedBy(0.5 * duration)));
707         Assertions.assertTrue(maneuverWithoutOverride.isFiring(fireDate.shiftedBy(duration - 1.0e-6)));
708         Assertions.assertFalse(maneuverWithoutOverride.isFiring(fireDate.shiftedBy(duration + 1.0e-6)));
709 
710         // intentionally wrong propagation, that will produce a very different state
711         // propagator uses LOF attitude,
712         // maneuver forget to override it, so maneuver will be LOF-aligned in this case
713         AdaptiveStepsizeIntegrator integrator3 =
714                         new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
715         integrator3.setInitialStepSize(60);
716         final NumericalPropagator propagator3 = new NumericalPropagator(integrator3);
717         propagator3.setInitialState(initialState);
718         propagator3.setAttitudeProvider(lofLaw);
719         propagator3.addForceModel(maneuverWithoutOverride);
720         final SpacecraftState finalState3 = propagator3.propagate(finalState1.getDate());
721         Assertions.assertEquals(345859.0,
722                            Vector3D.distance(finalState1.getPosition(),
723                                              finalState3.getPosition()),
724                            1.0);
725     }
726 
727     @Test
728     void testStopInMiddle() {
729         final double isp = 318;
730         final double mass = 2500;
731         final double a = 24396159;
732         final double e = 0.72831215;
733         final double i = FastMath.toRadians(7);
734         final double omega = FastMath.toRadians(180);
735         final double OMEGA = FastMath.toRadians(261);
736         final double lv = 0;
737 
738         final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 01, 01),
739                                                        new TimeComponents(23, 30, 00.000),
740                                                        TimeScalesFactory.getUTC());
741         final Orbit orbit =
742             new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
743                                FramesFactory.getEME2000(), initDate, mu);
744 
745         final double duration = 3653.99;
746         final double f = 420;
747         final double delta = FastMath.toRadians(-7.4978);
748         final double alpha = FastMath.toRadians(351);
749         final AttitudeProvider inertialLaw = new FrameAlignedProvider(new Rotation(new Vector3D(alpha, delta),
750                                                                                    Vector3D.PLUS_I));
751 
752         final SpacecraftState initialState =
753             new SpacecraftState(orbit, inertialLaw.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
754 
755         final AbsoluteDate fireDate = new AbsoluteDate(new DateComponents(2004, 01, 02),
756                                                        new TimeComponents(04, 15, 34.080),
757                                                        TimeScalesFactory.getUTC());
758         final ConstantThrustManeuver maneuver =
759             new ConstantThrustManeuver(fireDate, duration, f, isp, Vector3D.PLUS_I);
760         Assertions.assertEquals(f,   maneuver.getThrustMagnitude(), 1.0e-10);
761         Assertions.assertEquals(isp, maneuver.getIsp(),    1.0e-10);
762 
763         // before events have been encountered, the maneuver is not yet allowed to generate non zero acceleration
764         for (double dt = 0 ; dt < fireDate.durationFrom(initDate) + duration + 100; dt += 0.1) {
765             Assertions.assertFalse(maneuver.isFiring(initDate.shiftedBy(dt)));
766             Assertions.assertEquals(0.0,
767                                 maneuver.acceleration(initialState.shiftedBy(dt), maneuver.getParameters()).getNorm(),
768                                 1.0e-15);
769         }
770 
771         // reference propagation:
772         // propagator already uses inertial law
773         // maneuver does not need to override it to get an inertial maneuver
774         double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1.).getTolerances(orbit, OrbitType.KEPLERIAN);
775         AdaptiveStepsizeIntegrator integrator1 =
776             new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
777         integrator1.setInitialStepSize(60);
778         final NumericalPropagator propagator1 = new NumericalPropagator(integrator1);
779         propagator1.setInitialState(initialState);
780         propagator1.setAttitudeProvider(inertialLaw);
781         propagator1.addForceModel(maneuver);
782         final SpacecraftState middleState = propagator1.propagate(fireDate.shiftedBy(0.5 * duration));
783 
784         Assertions.assertFalse(maneuver.isFiring(initialState));
785         Assertions.assertEquals(0.0,
786                             maneuver.acceleration(initialState, maneuver.getParameters()).getNorm(),
787                             1.0e-15);
788         Assertions.assertTrue(maneuver.isFiring(middleState));
789         Assertions.assertEquals(0.186340263,
790                             maneuver.acceleration(middleState, maneuver.getParameters()).getNorm(),
791                             1.0e-9);
792         Assertions.assertTrue(maneuver.isFiring(middleState.shiftedBy(3 * duration)));
793         Assertions.assertEquals(0.186340263,
794                             maneuver.acceleration(middleState.shiftedBy(3 * duration), maneuver.getParameters()).getNorm(),
795                             1.0e-9);
796 
797 
798     }
799 
800     @Test
801     void testNullDuration() {
802 
803         // Defining initial state
804         final Frame                    eme2000      = FramesFactory.getEME2000();
805         final AbsoluteDate             initialDate  = new AbsoluteDate();
806         final TimeStampedPVCoordinates initialPV    = new TimeStampedPVCoordinates(initialDate,
807                                                                                    new Vector3D(6378e3 + 400e3, 0, 0),
808                                                                                    new Vector3D(0, 7669, 0));
809         final double         initialMass    = 1000;
810         final CartesianOrbit cartesianOrbit = new CartesianOrbit(initialPV, eme2000, Constants.EIGEN5C_EARTH_MU);
811 
812         // Defining ConstantThrustManeuver with null duration
813         final AbsoluteDate     startManeuverDate = initialDate.shiftedBy(30);
814         final double           duration          = 0;                                  // in s
815         final double           thrust            = 100;                                // default value
816         final double           isp               = 300;                                // default value
817         final Vector3D         direction         = new Vector3D(1, 0, 0);
818         final AttitudeProvider attitudeProvider  = new LofOffset(eme2000, LOFType.TNW);
819 
820         final ConstantThrustManeuver nullDurationManeuver = new ConstantThrustManeuver(startManeuverDate, duration,
821                 thrust, isp, attitudeProvider, direction);
822 
823         // Defining propagator
824         // Default Values
825         final double dP      = 0.001;
826         final double minStep = 0.1;
827         final double maxStep = 3600;
828 
829         // Defining integrator
830         final double[][]          tolerances = ToleranceProvider.getDefaultToleranceProvider(dP).getTolerances(cartesianOrbit, OrbitType.CARTESIAN);
831         final ODEIntegrator       integrator = new DormandPrince853Integrator(minStep, maxStep, tolerances[0],
832                 tolerances[1]);
833         final NumericalPropagator numProp = new NumericalPropagator(integrator);
834 
835         // Configuring propagator
836         numProp.setOrbitType(OrbitType.CARTESIAN);
837         numProp.setInitialState(new SpacecraftState(cartesianOrbit, initialMass));
838         numProp.addForceModel(nullDurationManeuver);
839 
840         // Propagation
841         final SpacecraftState finalState = numProp.propagate(initialDate.shiftedBy(60));
842         Assertions.assertEquals(cartesianOrbit.getA(), finalState.getOrbit().getA(), 1.0e-15);
843 
844     }
845 
846     @BeforeEach
847     public void setUp() {
848         Utils.setDataRoot("regular-data");
849 
850         // Body mu
851         mu = 3.9860047e14;
852 
853     }
854 
855 }