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.empirical;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.Field;
21  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
22  import org.hipparchus.geometry.euclidean.threed.Rotation;
23  import org.hipparchus.geometry.euclidean.threed.Vector3D;
24  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
25  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
26  import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
27  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
28  import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
29  import org.hipparchus.util.Binary64Field;
30  import org.hipparchus.util.FastMath;
31  import org.junit.jupiter.api.Assertions;
32  import org.junit.jupiter.api.BeforeEach;
33  import org.junit.jupiter.api.Test;
34  import org.orekit.Utils;
35  import org.orekit.attitudes.Attitude;
36  import org.orekit.attitudes.AttitudeProvider;
37  import org.orekit.attitudes.CelestialBodyPointed;
38  import org.orekit.attitudes.FrameAlignedProvider;
39  import org.orekit.attitudes.LofOffset;
40  import org.orekit.bodies.CelestialBodyFactory;
41  import org.orekit.errors.OrekitException;
42  import org.orekit.estimation.leastsquares.BatchLSEstimator;
43  import org.orekit.estimation.measurements.ObservableSatellite;
44  import org.orekit.estimation.measurements.ObservedMeasurement;
45  import org.orekit.estimation.measurements.PV;
46  import org.orekit.forces.AbstractForceModelTest;
47  import org.orekit.forces.maneuvers.ConstantThrustManeuver;
48  import org.orekit.frames.Frame;
49  import org.orekit.frames.FramesFactory;
50  import org.orekit.frames.LOFType;
51  import org.orekit.orbits.*;
52  import org.orekit.propagation.*;
53  import org.orekit.propagation.conversion.DormandPrince853IntegratorBuilder;
54  import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
55  import org.orekit.propagation.numerical.FieldNumericalPropagator;
56  import org.orekit.propagation.numerical.NumericalPropagator;
57  import org.orekit.propagation.sampling.MultiSatStepHandler;
58  import org.orekit.time.AbsoluteDate;
59  import org.orekit.time.DateComponents;
60  import org.orekit.time.FieldAbsoluteDate;
61  import org.orekit.time.TimeComponents;
62  import org.orekit.time.TimeScalesFactory;
63  import org.orekit.utils.Constants;
64  import org.orekit.utils.PVCoordinates;
65  import org.orekit.utils.ParameterDriver;
66  
67  import java.util.ArrayList;
68  import java.util.Arrays;
69  import java.util.List;
70  
71  public class HarmonicAccelerationModelTest extends AbstractForceModelTest {
72  
73      private Orbit initialOrbit;
74  
75      @Test
76      void testEquivalentInertialManeuver() {
77          final double   delta     = FastMath.toRadians(-7.4978);
78          final double   alpha     = FastMath.toRadians(351);
79          final Vector3D direction = new Vector3D(alpha, delta);
80          final double mass        = 2500;
81          final double isp         = Double.POSITIVE_INFINITY;
82          final double duration    = 4000;
83          final double f           = 400;
84  
85          final AttitudeProvider maneuverLaw = new FrameAlignedProvider(new Rotation(direction, Vector3D.PLUS_I));
86          ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
87                                                                       duration, f, isp, Vector3D.PLUS_I);
88          final AttitudeProvider accelerationLaw = new FrameAlignedProvider(new Rotation(direction, Vector3D.PLUS_K));
89          final AccelerationModel accelerationModel = new HarmonicAccelerationModel("", AbsoluteDate.J2000_EPOCH,
90                                                                                            Double.POSITIVE_INFINITY, 1);
91          final ParametricAcceleration inertialAcceleration = new ParametricAcceleration(direction, true, accelerationModel);
92          Assertions.assertTrue(inertialAcceleration.dependsOnPositionOnly());
93          inertialAcceleration.getParametersDrivers().get(0).setValue(f / mass);
94          inertialAcceleration.getParametersDrivers().get(1).setValue(0.5 * FastMath.PI);
95          doTestEquivalentManeuver(mass, maneuverLaw, maneuver, accelerationLaw, inertialAcceleration, 1.0e-15);
96      }
97  
98      @Test
99      void testEquivalentTangentialManeuver() {
100         final double mass     = 2500;
101         final double isp      = Double.POSITIVE_INFINITY;
102         final double duration = 4000;
103         final double f        = 400;
104 
105         final AttitudeProvider commonLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
106         ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
107                                                                      duration, f, isp, Vector3D.PLUS_I);
108         final AccelerationModel accelerationModel = new HarmonicAccelerationModel("", null,
109                                                                                           Double.POSITIVE_INFINITY, 1);
110         final ParametricAcceleration lofAcceleration = new ParametricAcceleration(Vector3D.PLUS_I, false, accelerationModel);
111         Assertions.assertFalse(lofAcceleration.dependsOnPositionOnly());
112         lofAcceleration.getParametersDrivers().get(0).setValue(f / mass);
113         lofAcceleration.getParametersDrivers().get(1).setValue(0.5 * FastMath.PI);
114         doTestEquivalentManeuver(mass, commonLaw, maneuver, commonLaw, lofAcceleration, 1.0e-15);
115     }
116 
117     @Test
118     void testEquivalentTangentialOverriddenManeuver() {
119         final double mass     = 2500;
120         final double isp      = Double.POSITIVE_INFINITY;
121         final double duration = 4000;
122         final double f        = 400;
123 
124         final AttitudeProvider maneuverLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
125         ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
126                                                                      duration, f, isp, Vector3D.PLUS_I);
127         final AttitudeProvider accelerationLaw = new CelestialBodyPointed(initialOrbit.getFrame(),
128                                                                           CelestialBodyFactory.getSun(), Vector3D.PLUS_K,
129                                                                           Vector3D.PLUS_I, Vector3D.PLUS_K);
130         final AccelerationModel accelerationModel = new HarmonicAccelerationModel("prefix", null,
131                                                                                           Double.POSITIVE_INFINITY, 1);
132         final ParametricAcceleration lofAcceleration = new ParametricAcceleration(Vector3D.PLUS_I, maneuverLaw, accelerationModel);
133         lofAcceleration.getParametersDrivers().get(0).setValue(f / mass);
134         lofAcceleration.getParametersDrivers().get(1).setValue(0.5 * FastMath.PI);
135         doTestEquivalentManeuver(mass, maneuverLaw, maneuver, accelerationLaw, lofAcceleration, 1.0e-15);
136     }
137 
138     private void doTestEquivalentManeuver(final double mass,
139                                           final AttitudeProvider maneuverLaw,
140                                           final ConstantThrustManeuver maneuver,
141                                           final AttitudeProvider accelerationLaw,
142                                           final ParametricAcceleration parametricAcceleration,
143                                           final double positionTolerance) {
144 
145         SpacecraftState initialState = new SpacecraftState(initialOrbit,
146                                                            maneuverLaw.getAttitude(initialOrbit,
147                                                                                    initialOrbit.getDate(),
148                                                                                    initialOrbit.getFrame()),
149                                                            mass);
150 
151         double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(initialOrbit, initialOrbit.getType());
152 
153         // propagator 0 uses a maneuver that is so efficient it does not consume any fuel
154         // (hence mass remains constant)
155         AdaptiveStepsizeIntegrator integrator0 =
156             new DormandPrince853Integrator(0.001, 100, tolerance[0], tolerance[1]);
157         integrator0.setInitialStepSize(60);
158         final NumericalPropagator propagator0 = new NumericalPropagator(integrator0);
159         propagator0.setOrbitType(OrbitType.EQUINOCTIAL);
160         propagator0.setPositionAngleType(PositionAngleType.TRUE);
161         propagator0.setInitialState(initialState);
162         propagator0.setAttitudeProvider(maneuverLaw);
163         propagator0.addForceModel(maneuver);
164 
165         // propagator 1 uses a constant acceleration
166         AdaptiveStepsizeIntegrator integrator1 =
167                         new DormandPrince853Integrator(0.001, 100, tolerance[0], tolerance[1]);
168         integrator1.setInitialStepSize(60);
169         final NumericalPropagator propagator1 = new NumericalPropagator(integrator1);
170         propagator1.setOrbitType(propagator0.getOrbitType());
171         propagator1.setPositionAngleType(propagator0.getPositionAngleType());
172         propagator1.setInitialState(initialState);
173         propagator1.setAttitudeProvider(accelerationLaw);
174         propagator1.addForceModel(parametricAcceleration);
175 
176         MultiSatStepHandler handler = interpolators -> {
177             Vector3D p0 = interpolators.get(0).getCurrentState().getPosition();
178             Vector3D p1 = interpolators.get(1).getCurrentState().getPosition();
179             Assertions.assertEquals(0.0, Vector3D.distance(p0, p1), positionTolerance);
180         };
181         PropagatorsParallelizer parallelizer = new PropagatorsParallelizer(Arrays.asList(propagator0, propagator1),
182                                                                            handler);
183 
184         parallelizer.propagate(initialOrbit.getDate(), initialOrbit.getDate().shiftedBy(1000.0));
185 
186     }
187 
188     @Test
189     void testEquivalentInertialManeuverField() {
190         final double   delta     = FastMath.toRadians(-7.4978);
191         final double   alpha     = FastMath.toRadians(351);
192         final Vector3D direction = new Vector3D(alpha, delta);
193         final double mass        = 2500;
194         final double isp         = Double.POSITIVE_INFINITY;
195         final double duration    = 4000;
196         final double f           = 400;
197 
198         final AttitudeProvider maneuverLaw = new FrameAlignedProvider(new Rotation(direction, Vector3D.PLUS_I));
199         ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
200                                                                      duration, f, isp, Vector3D.PLUS_I);
201         final AttitudeProvider accelerationLaw = new FrameAlignedProvider(new Rotation(direction, Vector3D.PLUS_K));
202         final AccelerationModel accelerationModel = new HarmonicAccelerationModel("", AbsoluteDate.J2000_EPOCH,
203                                                                                           Double.POSITIVE_INFINITY, 1);
204         final ParametricAcceleration inertialAcceleration = new ParametricAcceleration(direction, true, accelerationModel);
205         inertialAcceleration.getParametersDrivers().get(0).setValue(f / mass);
206         inertialAcceleration.getParametersDrivers().get(1).setValue(0.5 * FastMath.PI);
207         doTestEquivalentManeuver(Binary64Field.getInstance(),
208                                  mass, maneuverLaw, maneuver, accelerationLaw, inertialAcceleration, 3.0e-9);
209     }
210 
211     @Test
212     void testEquivalentTangentialManeuverField() {
213         final double mass     = 2500;
214         final double isp      = Double.POSITIVE_INFINITY;
215         final double duration = 4000;
216         final double f        = 400;
217 
218         final AttitudeProvider commonLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
219         ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
220                                                                      duration, f, isp, Vector3D.PLUS_I);
221         final HarmonicAccelerationModel accelerationModel = new HarmonicAccelerationModel("", null,
222                                                                                           Double.POSITIVE_INFINITY, 1);
223         final ParametricAcceleration lofAcceleration = new ParametricAcceleration(Vector3D.PLUS_I, false, accelerationModel);
224         lofAcceleration.getParametersDrivers().get(0).setValue(f / mass);
225         lofAcceleration.getParametersDrivers().get(1).setValue(0.5 * FastMath.PI);
226         doTestEquivalentManeuver(Binary64Field.getInstance(),
227                                  mass, commonLaw, maneuver, commonLaw, lofAcceleration, 1.0e-15);
228     }
229 
230     @Test
231     void testEquivalentTangentialOverriddenManeuverField() {
232         final double mass     = 2500;
233         final double isp      = Double.POSITIVE_INFINITY;
234         final double duration = 4000;
235         final double f        = 400;
236 
237         final AttitudeProvider maneuverLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
238         ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
239                                                                      duration, f, isp, Vector3D.PLUS_I);
240         final AttitudeProvider accelerationLaw = new CelestialBodyPointed(initialOrbit.getFrame(),
241                                                                           CelestialBodyFactory.getSun(), Vector3D.PLUS_K,
242                                                                           Vector3D.PLUS_I, Vector3D.PLUS_K);
243         final HarmonicAccelerationModel accelerationModel = new HarmonicAccelerationModel( "prefix", null,
244                                                                                            Double.POSITIVE_INFINITY, 1);
245         final ParametricAcceleration lofAcceleration = new ParametricAcceleration(Vector3D.PLUS_I, maneuverLaw, accelerationModel);
246         lofAcceleration.getParametersDrivers().get(0).setValue(f / mass);
247         lofAcceleration.getParametersDrivers().get(1).setValue(0.5 * FastMath.PI);
248         doTestEquivalentManeuver(Binary64Field.getInstance(),
249                                  mass, maneuverLaw, maneuver, accelerationLaw, lofAcceleration, 1.0e-15);
250     }
251 
252     private <T extends CalculusFieldElement<T>> void doTestEquivalentManeuver(final Field<T> field,
253                                                                           final double mass,
254                                                                           final AttitudeProvider maneuverLaw,
255                                                                           final ConstantThrustManeuver maneuver,
256                                                                           final AttitudeProvider accelerationLaw,
257                                                                           final ParametricAcceleration parametricAcceleration,
258                                                                           final double positionTolerance)
259                                                                                           {
260 
261         FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(field,
262                                                                           new SpacecraftState(initialOrbit,
263                                                                                               maneuverLaw.getAttitude(initialOrbit,
264                                                                                                                       initialOrbit.getDate(),
265                                                                                                                       initialOrbit.getFrame()),
266                                                                                               mass));
267 
268         double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(10).getTolerances(
269                                                                    initialState.getOrbit(),
270                                                                    initialState.getOrbit().getType());
271 
272         // propagator 0 uses a maneuver that is so efficient it does not consume any fuel
273         // (hence mass remains constant)
274         AdaptiveStepsizeFieldIntegrator<T> integrator0 =
275             new DormandPrince853FieldIntegrator<>(field, 0.001, 100, tolerance[0], tolerance[1]);
276         integrator0.setInitialStepSize(60);
277         final FieldNumericalPropagator<T> propagator0 = new FieldNumericalPropagator<>(field, integrator0);
278         propagator0.setOrbitType(OrbitType.EQUINOCTIAL);
279         propagator0.setPositionAngleType(PositionAngleType.TRUE);
280         propagator0.setInitialState(initialState);
281         propagator0.setAttitudeProvider(maneuverLaw);
282         propagator0.addForceModel(maneuver);
283         final FieldEphemerisGenerator<T> generator1 = propagator0.getEphemerisGenerator();
284         propagator0.propagate(initialState.getDate(), initialState.getDate().shiftedBy(1000.0));
285         FieldBoundedPropagator<T> ephemeris0 = generator1.getGeneratedEphemeris();
286 
287         // propagator 1 uses a constant acceleration
288         AdaptiveStepsizeFieldIntegrator<T> integrator1 =
289                         new DormandPrince853FieldIntegrator<>(field, 0.001, 100, tolerance[0], tolerance[1]);
290         integrator1.setInitialStepSize(60);
291         final FieldNumericalPropagator<T> propagator1 = new FieldNumericalPropagator<>(field, integrator1);
292         propagator1.setOrbitType(propagator0.getOrbitType());
293         propagator1.setPositionAngleType(propagator0.getPositionAngleType());
294         propagator1.setInitialState(initialState);
295         propagator1.setAttitudeProvider(accelerationLaw);
296         propagator1.addForceModel(parametricAcceleration);
297         final FieldEphemerisGenerator<T> generator2 = propagator1.getEphemerisGenerator();
298         propagator1.propagate(initialState.getDate(), initialState.getDate().shiftedBy(1000.0));
299         FieldBoundedPropagator<T> ephemeris1 = generator2.getGeneratedEphemeris();
300 
301         for (double dt = 1; dt < 999; dt += 10) {
302             FieldAbsoluteDate<T> t = initialState.getDate().shiftedBy(dt);
303             FieldVector3D<T> p0 = ephemeris0.propagate(t).getPosition();
304             FieldVector3D<T> p1 = ephemeris1.propagate(t).getPosition();
305             Assertions.assertEquals(0, FieldVector3D.distance(p0, p1).getReal(), positionTolerance);
306         }
307 
308     }
309 
310     @Test
311     void testParameterDerivative1T() {
312         doTestParameterDerivative(1, 4.0e-14, 2.0e-11);
313     }
314 
315     @Test
316     void testParameterDerivative2T() {
317         doTestParameterDerivative(2, 3.0e-14, 7.0e-12);
318     }
319 
320     @Test
321     void testParameterDerivative3T() {
322         doTestParameterDerivative(3, 2.0e-14, 2.0e-11);
323     }
324 
325     private void doTestParameterDerivative(final int harmonicMultiplier,
326                                            final double amplitudeDerivativeTolerance,
327                                            final double phaseDerivativeTolerance)
328         {
329 
330         // pos-vel (from a ZOOM ephemeris reference)
331         final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
332         final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
333         final Frame frame = FramesFactory.getGCRF();
334         final AbsoluteDate date = new AbsoluteDate(2005, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI());
335         final CartesianOrbit orbit = new CartesianOrbit(new PVCoordinates(pos, vel), frame,
336                 date, Constants.EIGEN5C_EARTH_MU);
337         final LofOffset lofOffset = new LofOffset(frame, LOFType.LVLH_CCSDS);
338         final Attitude attitude = lofOffset.getAttitude(orbit, date, frame);  // necessary for non-regression
339         final SpacecraftState state = new SpacecraftState(orbit, attitude);
340 
341         final HarmonicAccelerationModel accelerationModel = new HarmonicAccelerationModel("kT",
342                                                                                           state.getDate().shiftedBy(-2.0),
343                                                                                           state.getOrbit().getKeplerianPeriod(), harmonicMultiplier);
344         final ParametricAcceleration hpa = new ParametricAcceleration(Vector3D.PLUS_K, false, accelerationModel);
345         hpa.init(state, state.getDate().shiftedBy(3600.0));
346         hpa.getParametersDrivers().get(0).setValue(0.00001);
347         hpa.getParametersDrivers().get(1).setValue(0.00002);
348         checkParameterDerivative(state, hpa, "kT γ", 1.0e-3, amplitudeDerivativeTolerance);
349         checkParameterDerivative(state, hpa, "kT φ",     1.0e-3, phaseDerivativeTolerance);
350 
351     }
352 
353     @Test
354     void testCoefficientsDetermination() {
355 
356         final double mass = 2500;
357         final Orbit orbit = new CircularOrbit(7500000.0, 1.0e-4, 1.0e-3, 1.7, 0.3, 0.5, PositionAngleType.TRUE,
358                                         FramesFactory.getEME2000(),
359                                         new AbsoluteDate(new DateComponents(2004, 2, 3), TimeComponents.H00,
360                                                          TimeScalesFactory.getUTC()),
361                                         Constants.EIGEN5C_EARTH_MU);
362         final double period = orbit.getKeplerianPeriod();
363         AttitudeProvider maneuverLaw = new LofOffset(orbit.getFrame(), LOFType.VNC);
364         SpacecraftState initialState = new SpacecraftState(orbit,
365                                                            maneuverLaw.getAttitude(orbit,
366                                                                                    orbit.getDate(),
367                                                                                    orbit.getFrame()),
368                                                            mass);
369 
370         double dP      = 10.0;
371         double minStep = 0.001;
372         double maxStep = 100;
373         double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(dP).getTolerances(orbit, orbit.getType());
374 
375         // generate PV measurements corresponding to a tangential maneuver
376         AdaptiveStepsizeIntegrator integrator0 =
377             new DormandPrince853Integrator(minStep, maxStep, tolerance[0], tolerance[1]);
378         integrator0.setInitialStepSize(60);
379         final NumericalPropagator propagator0 = new NumericalPropagator(integrator0);
380         propagator0.setOrbitType(OrbitType.EQUINOCTIAL);
381         propagator0.setPositionAngleType(PositionAngleType.TRUE);
382         propagator0.setInitialState(initialState);
383         propagator0.setAttitudeProvider(maneuverLaw);
384         final ParametricAcceleration hpaRefX1 = new ParametricAcceleration(Vector3D.PLUS_I, true,
385                                                                            new HarmonicAccelerationModel("refX1", null, period, 1));
386         final ParametricAcceleration hpaRefY1 = new ParametricAcceleration(Vector3D.PLUS_J, true,
387                                                                            new HarmonicAccelerationModel("refY1", null, period, 1));
388         final ParametricAcceleration hpaRefZ2 = new ParametricAcceleration(Vector3D.PLUS_K, true,
389                                                                            new HarmonicAccelerationModel("refZ2", null, period, 2));
390         hpaRefX1.getParametersDrivers().get(0).setValue(2.4e-2);
391         hpaRefX1.getParametersDrivers().get(1).setValue(3.1);
392         hpaRefY1.getParametersDrivers().get(0).setValue(4.0e-2);
393         hpaRefY1.getParametersDrivers().get(1).setValue(0.3);
394         hpaRefZ2.getParametersDrivers().get(0).setValue(1.0e-2);
395         hpaRefZ2.getParametersDrivers().get(1).setValue(1.8);
396         propagator0.addForceModel(hpaRefX1);
397         propagator0.addForceModel(hpaRefY1);
398         propagator0.addForceModel(hpaRefZ2);
399         ObservableSatellite sat0 = new ObservableSatellite(0);
400         final List<ObservedMeasurement<?>> measurements = new ArrayList<>();
401         propagator0.setStepHandler(10.0,
402                                    state ->
403                                    measurements.add(new PV(state.getDate(),
404                                                            state.getPosition(), state.getPVCoordinates().getVelocity(),
405                                                            1.0e-3, 1.0e-6, 1.0, sat0)));
406         propagator0.propagate(orbit.getDate().shiftedBy(900));
407 
408         // set up an estimator to retrieve the maneuver as several harmonic accelerations in inertial frame
409         final NumericalPropagatorBuilder propagatorBuilder =
410                         new NumericalPropagatorBuilder(orbit,
411                                                        new DormandPrince853IntegratorBuilder(minStep, maxStep, dP),
412                                                        PositionAngleType.TRUE, dP);
413         propagatorBuilder.addForceModel(new ParametricAcceleration(Vector3D.PLUS_I, true,
414                                                                    new HarmonicAccelerationModel("X1", null, period, 1)));
415         propagatorBuilder.addForceModel(new ParametricAcceleration(Vector3D.PLUS_J, true,
416                                                                    new HarmonicAccelerationModel("Y1", null, period, 1)));
417         propagatorBuilder.addForceModel(new ParametricAcceleration(Vector3D.PLUS_K, true,
418                                                                    new HarmonicAccelerationModel("Z2", null, period, 2)));
419         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), propagatorBuilder);
420         estimator.setParametersConvergenceThreshold(1.0e-2);
421         estimator.setMaxIterations(20);
422         estimator.setMaxEvaluations(100);
423         for (final ObservedMeasurement<?> measurement : measurements) {
424             estimator.addMeasurement(measurement);
425         }
426 
427         // we will estimate only the force model parameters, not the orbit
428         for (final ParameterDriver d : estimator.getOrbitalParametersDrivers(false).getDrivers()) {
429             d.setSelected(false);
430         }
431         setParameter(estimator, "X1 γ", 1.0e-2);
432         setParameter(estimator, "X1 φ", 4.0);
433         setParameter(estimator, "Y1 γ", 1.0e-2);
434         setParameter(estimator, "Y1 φ", 0.0);
435         setParameter(estimator, "Z2 γ", 1.0e-2);
436         setParameter(estimator, "Z2 φ", 1.0);
437 
438         estimator.estimate();
439         Assertions.assertTrue(estimator.getIterationsCount()  < 15);
440         Assertions.assertTrue(estimator.getEvaluationsCount() < 15);
441         Assertions.assertEquals(0.0, estimator.getOptimum().getRMS(), 1.0e-5);
442 
443         Assertions.assertEquals(hpaRefX1.getParametersDrivers().get(0).getValue(), getParameter(estimator, "X1 γ"), 1.e-12);
444         Assertions.assertEquals(hpaRefX1.getParametersDrivers().get(1).getValue(), getParameter(estimator, "X1 φ"), 1.e-12);
445         Assertions.assertEquals(hpaRefY1.getParametersDrivers().get(0).getValue(), getParameter(estimator, "Y1 γ"), 1.e-12);
446         Assertions.assertEquals(hpaRefY1.getParametersDrivers().get(1).getValue(), getParameter(estimator, "Y1 φ"), 1.e-12);
447         Assertions.assertEquals(hpaRefZ2.getParametersDrivers().get(0).getValue(), getParameter(estimator, "Z2 γ"), 1.e-12);
448         Assertions.assertEquals(hpaRefZ2.getParametersDrivers().get(1).getValue(), getParameter(estimator, "Z2 φ"), 1.e-12);
449 
450     }
451 
452     private void setParameter(BatchLSEstimator estimator, String name, double value)
453         {
454         for (final ParameterDriver driver : estimator.getPropagatorParametersDrivers(false).getDrivers()) {
455             if (driver.getName().equals(name)) {
456                 driver.setSelected(true);
457                 driver.setValue(value);
458                 return;
459             }
460         }
461         Assertions.fail("unknown parameter " + name);
462     }
463 
464     // if Pdriver has only 1 value driven
465     private double getParameter(BatchLSEstimator estimator, String name)
466     {
467     for (final ParameterDriver driver : estimator.getPropagatorParametersDrivers(false).getDrivers()) {
468         if (driver.getName().equals(name)) {
469             return driver.getValue();
470         }
471     }
472     Assertions.fail("unknown parameter " + name);
473     return Double.NaN;
474     }
475 
476     @BeforeEach
477     public void setUp() {
478         try {
479             Utils.setDataRoot("regular-data");
480 
481             final double a = 24396159;
482             final double e = 0.72831215;
483             final double i = FastMath.toRadians(7);
484             final double omega = FastMath.toRadians(180);
485             final double OMEGA = FastMath.toRadians(261);
486             final double lv = 0;
487 
488             final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 1, 1),
489                                                            new TimeComponents(23, 30, 00.000),
490                                                            TimeScalesFactory.getUTC());
491             initialOrbit =
492                             new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
493                                                FramesFactory.getEME2000(), initDate, Constants.EIGEN5C_EARTH_MU);
494         } catch (OrekitException oe) {
495             Assertions.fail(oe.getLocalizedMessage());
496         }
497 
498     }
499 
500 }