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