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())).withMass(mass);
149 
150         double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(initialOrbit, initialOrbit.getType());
151 
152         // propagator 0 uses a maneuver that is so efficient it does not consume any fuel
153         // (hence mass remains constant)
154         AdaptiveStepsizeIntegrator integrator0 =
155             new DormandPrince853Integrator(0.001, 100, tolerance[0], tolerance[1]);
156         integrator0.setInitialStepSize(60);
157         final NumericalPropagator propagator0 = new NumericalPropagator(integrator0);
158         propagator0.setOrbitType(OrbitType.EQUINOCTIAL);
159         propagator0.setPositionAngleType(PositionAngleType.TRUE);
160         propagator0.setInitialState(initialState);
161         propagator0.setAttitudeProvider(maneuverLaw);
162         propagator0.addForceModel(maneuver);
163 
164         // propagator 1 uses a constant acceleration
165         AdaptiveStepsizeIntegrator integrator1 =
166                         new DormandPrince853Integrator(0.001, 100, tolerance[0], tolerance[1]);
167         integrator1.setInitialStepSize(60);
168         final NumericalPropagator propagator1 = new NumericalPropagator(integrator1);
169         propagator1.setOrbitType(propagator0.getOrbitType());
170         propagator1.setPositionAngleType(propagator0.getPositionAngleType());
171         propagator1.setInitialState(initialState);
172         propagator1.setAttitudeProvider(accelerationLaw);
173         propagator1.addForceModel(parametricAcceleration);
174 
175         MultiSatStepHandler handler = interpolators -> {
176             Vector3D p0 = interpolators.get(0).getCurrentState().getPosition();
177             Vector3D p1 = interpolators.get(1).getCurrentState().getPosition();
178             Assertions.assertEquals(0.0, Vector3D.distance(p0, p1), positionTolerance);
179         };
180         PropagatorsParallelizer parallelizer = new PropagatorsParallelizer(Arrays.asList(propagator0, propagator1),
181                                                                            handler);
182 
183         parallelizer.propagate(initialOrbit.getDate(), initialOrbit.getDate().shiftedBy(1000.0));
184 
185     }
186 
187     @Test
188     void testEquivalentInertialManeuverField() {
189         final double   delta     = FastMath.toRadians(-7.4978);
190         final double   alpha     = FastMath.toRadians(351);
191         final Vector3D direction = new Vector3D(alpha, delta);
192         final double mass        = 2500;
193         final double isp         = Double.POSITIVE_INFINITY;
194         final double duration    = 4000;
195         final double f           = 400;
196 
197         final AttitudeProvider maneuverLaw = new FrameAlignedProvider(new Rotation(direction, Vector3D.PLUS_I));
198         ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
199                                                                      duration, f, isp, Vector3D.PLUS_I);
200         final AttitudeProvider accelerationLaw = new FrameAlignedProvider(new Rotation(direction, Vector3D.PLUS_K));
201         final AccelerationModel accelerationModel = new HarmonicAccelerationModel("", AbsoluteDate.J2000_EPOCH,
202                                                                                           Double.POSITIVE_INFINITY, 1);
203         final ParametricAcceleration inertialAcceleration = new ParametricAcceleration(direction, true, accelerationModel);
204         inertialAcceleration.getParametersDrivers().get(0).setValue(f / mass);
205         inertialAcceleration.getParametersDrivers().get(1).setValue(0.5 * FastMath.PI);
206         doTestEquivalentManeuver(Binary64Field.getInstance(),
207                                  mass, maneuverLaw, maneuver, accelerationLaw, inertialAcceleration, 3.0e-9);
208     }
209 
210     @Test
211     void testEquivalentTangentialManeuverField() {
212         final double mass     = 2500;
213         final double isp      = Double.POSITIVE_INFINITY;
214         final double duration = 4000;
215         final double f        = 400;
216 
217         final AttitudeProvider commonLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
218         ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
219                                                                      duration, f, isp, Vector3D.PLUS_I);
220         final HarmonicAccelerationModel accelerationModel = new HarmonicAccelerationModel("", null,
221                                                                                           Double.POSITIVE_INFINITY, 1);
222         final ParametricAcceleration lofAcceleration = new ParametricAcceleration(Vector3D.PLUS_I, false, accelerationModel);
223         lofAcceleration.getParametersDrivers().get(0).setValue(f / mass);
224         lofAcceleration.getParametersDrivers().get(1).setValue(0.5 * FastMath.PI);
225         doTestEquivalentManeuver(Binary64Field.getInstance(),
226                                  mass, commonLaw, maneuver, commonLaw, lofAcceleration, 1.0e-15);
227     }
228 
229     @Test
230     void testEquivalentTangentialOverriddenManeuverField() {
231         final double mass     = 2500;
232         final double isp      = Double.POSITIVE_INFINITY;
233         final double duration = 4000;
234         final double f        = 400;
235 
236         final AttitudeProvider maneuverLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
237         ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
238                                                                      duration, f, isp, Vector3D.PLUS_I);
239         final AttitudeProvider accelerationLaw = new CelestialBodyPointed(initialOrbit.getFrame(),
240                                                                           CelestialBodyFactory.getSun(), Vector3D.PLUS_K,
241                                                                           Vector3D.PLUS_I, Vector3D.PLUS_K);
242         final HarmonicAccelerationModel accelerationModel = new HarmonicAccelerationModel( "prefix", null,
243                                                                                            Double.POSITIVE_INFINITY, 1);
244         final ParametricAcceleration lofAcceleration = new ParametricAcceleration(Vector3D.PLUS_I, maneuverLaw, accelerationModel);
245         lofAcceleration.getParametersDrivers().get(0).setValue(f / mass);
246         lofAcceleration.getParametersDrivers().get(1).setValue(0.5 * FastMath.PI);
247         doTestEquivalentManeuver(Binary64Field.getInstance(),
248                                  mass, maneuverLaw, maneuver, accelerationLaw, lofAcceleration, 1.0e-15);
249     }
250 
251     private <T extends CalculusFieldElement<T>> void doTestEquivalentManeuver(final Field<T> field,
252                                                                           final double mass,
253                                                                           final AttitudeProvider maneuverLaw,
254                                                                           final ConstantThrustManeuver maneuver,
255                                                                           final AttitudeProvider accelerationLaw,
256                                                                           final ParametricAcceleration parametricAcceleration,
257                                                                           final double positionTolerance)
258                                                                                           {
259 
260         FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(field,
261                                                                           new SpacecraftState(initialOrbit,
262                                                                                               maneuverLaw.getAttitude(initialOrbit,
263                                                                                                                       initialOrbit.getDate(),
264                                                                                                                       initialOrbit.getFrame())).withMass(mass));
265 
266         double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(10).getTolerances(
267                                                                    initialState.getOrbit(),
268                                                                    initialState.getOrbit().getType());
269 
270         // propagator 0 uses a maneuver that is so efficient it does not consume any fuel
271         // (hence mass remains constant)
272         AdaptiveStepsizeFieldIntegrator<T> integrator0 =
273             new DormandPrince853FieldIntegrator<>(field, 0.001, 100, tolerance[0], tolerance[1]);
274         integrator0.setInitialStepSize(60);
275         final FieldNumericalPropagator<T> propagator0 = new FieldNumericalPropagator<>(field, integrator0);
276         propagator0.setOrbitType(OrbitType.EQUINOCTIAL);
277         propagator0.setPositionAngleType(PositionAngleType.TRUE);
278         propagator0.setInitialState(initialState);
279         propagator0.setAttitudeProvider(maneuverLaw);
280         propagator0.addForceModel(maneuver);
281         final FieldEphemerisGenerator<T> generator1 = propagator0.getEphemerisGenerator();
282         propagator0.propagate(initialState.getDate(), initialState.getDate().shiftedBy(1000.0));
283         FieldBoundedPropagator<T> ephemeris0 = generator1.getGeneratedEphemeris();
284 
285         // propagator 1 uses a constant acceleration
286         AdaptiveStepsizeFieldIntegrator<T> integrator1 =
287                         new DormandPrince853FieldIntegrator<>(field, 0.001, 100, tolerance[0], tolerance[1]);
288         integrator1.setInitialStepSize(60);
289         final FieldNumericalPropagator<T> propagator1 = new FieldNumericalPropagator<>(field, integrator1);
290         propagator1.setOrbitType(propagator0.getOrbitType());
291         propagator1.setPositionAngleType(propagator0.getPositionAngleType());
292         propagator1.setInitialState(initialState);
293         propagator1.setAttitudeProvider(accelerationLaw);
294         propagator1.addForceModel(parametricAcceleration);
295         final FieldEphemerisGenerator<T> generator2 = propagator1.getEphemerisGenerator();
296         propagator1.propagate(initialState.getDate(), initialState.getDate().shiftedBy(1000.0));
297         FieldBoundedPropagator<T> ephemeris1 = generator2.getGeneratedEphemeris();
298 
299         for (double dt = 1; dt < 999; dt += 10) {
300             FieldAbsoluteDate<T> t = initialState.getDate().shiftedBy(dt);
301             FieldVector3D<T> p0 = ephemeris0.propagate(t).getPosition();
302             FieldVector3D<T> p1 = ephemeris1.propagate(t).getPosition();
303             Assertions.assertEquals(0, FieldVector3D.distance(p0, p1).getReal(), positionTolerance);
304         }
305 
306     }
307 
308     @Test
309     void testParameterDerivative1T() {
310         doTestParameterDerivative(1, 4.0e-14, 2.0e-11);
311     }
312 
313     @Test
314     void testParameterDerivative2T() {
315         doTestParameterDerivative(2, 3.0e-14, 7.0e-12);
316     }
317 
318     @Test
319     void testParameterDerivative3T() {
320         doTestParameterDerivative(3, 2.0e-14, 2.0e-11);
321     }
322 
323     private void doTestParameterDerivative(final int harmonicMultiplier,
324                                            final double amplitudeDerivativeTolerance,
325                                            final double phaseDerivativeTolerance)
326         {
327 
328         // pos-vel (from a ZOOM ephemeris reference)
329         final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
330         final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
331         final Frame frame = FramesFactory.getGCRF();
332         final AbsoluteDate date = new AbsoluteDate(2005, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI());
333         final CartesianOrbit orbit = new CartesianOrbit(new PVCoordinates(pos, vel), frame,
334                 date, Constants.EIGEN5C_EARTH_MU);
335         final LofOffset lofOffset = new LofOffset(frame, LOFType.LVLH_CCSDS);
336         final Attitude attitude = lofOffset.getAttitude(orbit, date, frame);  // necessary for non-regression
337         final SpacecraftState state = new SpacecraftState(orbit, attitude);
338 
339         final HarmonicAccelerationModel accelerationModel = new HarmonicAccelerationModel("kT",
340                                                                                           state.getDate().shiftedBy(-2.0),
341                                                                                           state.getOrbit().getKeplerianPeriod(), harmonicMultiplier);
342         final ParametricAcceleration hpa = new ParametricAcceleration(Vector3D.PLUS_K, false, accelerationModel);
343         hpa.init(state, state.getDate().shiftedBy(3600.0));
344         hpa.getParametersDrivers().get(0).setValue(0.00001);
345         hpa.getParametersDrivers().get(1).setValue(0.00002);
346         checkParameterDerivative(state, hpa, "kT γ", 1.0e-3, amplitudeDerivativeTolerance);
347         checkParameterDerivative(state, hpa, "kT φ",     1.0e-3, phaseDerivativeTolerance);
348 
349     }
350 
351     @Test
352     void testCoefficientsDetermination() {
353 
354         final double mass = 2500;
355         final Orbit orbit = new CircularOrbit(7500000.0, 1.0e-4, 1.0e-3, 1.7, 0.3, 0.5, PositionAngleType.TRUE,
356                                         FramesFactory.getEME2000(),
357                                         new AbsoluteDate(new DateComponents(2004, 2, 3), TimeComponents.H00,
358                                                          TimeScalesFactory.getUTC()),
359                                         Constants.EIGEN5C_EARTH_MU);
360         final double period = orbit.getKeplerianPeriod();
361         AttitudeProvider maneuverLaw = new LofOffset(orbit.getFrame(), LOFType.VNC);
362         SpacecraftState initialState = new SpacecraftState(orbit,
363                                                            maneuverLaw.getAttitude(orbit,
364                                                                                    orbit.getDate(),
365                                                                                    orbit.getFrame())).withMass(mass);
366 
367         double dP      = 10.0;
368         double minStep = 0.001;
369         double maxStep = 100;
370         double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(dP).getTolerances(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.setOrbitType(OrbitType.EQUINOCTIAL);
378         propagator0.setPositionAngleType(PositionAngleType.TRUE);
379         propagator0.setInitialState(initialState);
380         propagator0.setAttitudeProvider(maneuverLaw);
381         final ParametricAcceleration hpaRefX1 = new ParametricAcceleration(Vector3D.PLUS_I, true,
382                                                                            new HarmonicAccelerationModel("refX1", null, period, 1));
383         final ParametricAcceleration hpaRefY1 = new ParametricAcceleration(Vector3D.PLUS_J, true,
384                                                                            new HarmonicAccelerationModel("refY1", null, period, 1));
385         final ParametricAcceleration hpaRefZ2 = new ParametricAcceleration(Vector3D.PLUS_K, true,
386                                                                            new HarmonicAccelerationModel("refZ2", null, period, 2));
387         hpaRefX1.getParametersDrivers().get(0).setValue(2.4e-2);
388         hpaRefX1.getParametersDrivers().get(1).setValue(3.1);
389         hpaRefY1.getParametersDrivers().get(0).setValue(4.0e-2);
390         hpaRefY1.getParametersDrivers().get(1).setValue(0.3);
391         hpaRefZ2.getParametersDrivers().get(0).setValue(1.0e-2);
392         hpaRefZ2.getParametersDrivers().get(1).setValue(1.8);
393         propagator0.addForceModel(hpaRefX1);
394         propagator0.addForceModel(hpaRefY1);
395         propagator0.addForceModel(hpaRefZ2);
396         ObservableSatellite sat0 = new ObservableSatellite(0);
397         final List<ObservedMeasurement<?>> measurements = new ArrayList<>();
398         propagator0.setStepHandler(10.0,
399                                    state ->
400                                    measurements.add(new PV(state.getDate(),
401                                                            state.getPosition(), state.getVelocity(),
402                                                            1.0e-3, 1.0e-6, 1.0, sat0)));
403         propagator0.propagate(orbit.getDate().shiftedBy(900));
404 
405         // set up an estimator to retrieve the maneuver as several harmonic accelerations in inertial frame
406         final NumericalPropagatorBuilder propagatorBuilder =
407                         new NumericalPropagatorBuilder(orbit,
408                                                        new DormandPrince853IntegratorBuilder(minStep, maxStep, dP),
409                                                        PositionAngleType.TRUE, dP);
410         propagatorBuilder.addForceModel(new ParametricAcceleration(Vector3D.PLUS_I, true,
411                                                                    new HarmonicAccelerationModel("X1", null, period, 1)));
412         propagatorBuilder.addForceModel(new ParametricAcceleration(Vector3D.PLUS_J, true,
413                                                                    new HarmonicAccelerationModel("Y1", null, period, 1)));
414         propagatorBuilder.addForceModel(new ParametricAcceleration(Vector3D.PLUS_K, true,
415                                                                    new HarmonicAccelerationModel("Z2", null, period, 2)));
416         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), propagatorBuilder);
417         estimator.setParametersConvergenceThreshold(1.0e-2);
418         estimator.setMaxIterations(20);
419         estimator.setMaxEvaluations(100);
420         for (final ObservedMeasurement<?> measurement : measurements) {
421             estimator.addMeasurement(measurement);
422         }
423 
424         // we will estimate only the force model parameters, not the orbit
425         for (final ParameterDriver d : estimator.getOrbitalParametersDrivers(false).getDrivers()) {
426             d.setSelected(false);
427         }
428         setParameter(estimator, "X1 γ", 1.0e-2);
429         setParameter(estimator, "X1 φ", 4.0);
430         setParameter(estimator, "Y1 γ", 1.0e-2);
431         setParameter(estimator, "Y1 φ", 0.0);
432         setParameter(estimator, "Z2 γ", 1.0e-2);
433         setParameter(estimator, "Z2 φ", 1.0);
434 
435         estimator.estimate();
436         Assertions.assertTrue(estimator.getIterationsCount()  < 15);
437         Assertions.assertTrue(estimator.getEvaluationsCount() < 15);
438         Assertions.assertEquals(0.0, estimator.getOptimum().getRMS(), 1.6e-5);
439 
440         Assertions.assertEquals(hpaRefX1.getParametersDrivers().get(0).getValue(), getParameter(estimator, "X1 γ"), 1.e-12);
441         Assertions.assertEquals(hpaRefX1.getParametersDrivers().get(1).getValue(), getParameter(estimator, "X1 φ"), 1.e-12);
442         Assertions.assertEquals(hpaRefY1.getParametersDrivers().get(0).getValue(), getParameter(estimator, "Y1 γ"), 1.e-12);
443         Assertions.assertEquals(hpaRefY1.getParametersDrivers().get(1).getValue(), getParameter(estimator, "Y1 φ"), 1.e-12);
444         Assertions.assertEquals(hpaRefZ2.getParametersDrivers().get(0).getValue(), getParameter(estimator, "Z2 γ"), 1.e-12);
445         Assertions.assertEquals(hpaRefZ2.getParametersDrivers().get(1).getValue(), getParameter(estimator, "Z2 φ"), 1.e-12);
446 
447     }
448 
449     private void setParameter(BatchLSEstimator estimator, String name, double value)
450         {
451         for (final ParameterDriver driver : estimator.getPropagatorParametersDrivers(false).getDrivers()) {
452             if (driver.getName().equals(name)) {
453                 driver.setSelected(true);
454                 driver.setValue(value);
455                 return;
456             }
457         }
458         Assertions.fail("unknown parameter " + name);
459     }
460 
461     // if Pdriver has only 1 value driven
462     private double getParameter(BatchLSEstimator estimator, String name)
463     {
464     for (final ParameterDriver driver : estimator.getPropagatorParametersDrivers(false).getDrivers()) {
465         if (driver.getName().equals(name)) {
466             return driver.getValue();
467         }
468     }
469     Assertions.fail("unknown parameter " + name);
470     return Double.NaN;
471     }
472 
473     @BeforeEach
474     public void setUp() {
475         try {
476             Utils.setDataRoot("regular-data");
477 
478             final double a = 24396159;
479             final double e = 0.72831215;
480             final double i = FastMath.toRadians(7);
481             final double omega = FastMath.toRadians(180);
482             final double OMEGA = FastMath.toRadians(261);
483             final double lv = 0;
484 
485             final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 1, 1),
486                                                            new TimeComponents(23, 30, 00.000),
487                                                            TimeScalesFactory.getUTC());
488             initialOrbit =
489                             new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
490                                                FramesFactory.getEME2000(), initDate, Constants.EIGEN5C_EARTH_MU);
491         } catch (OrekitException oe) {
492             Assertions.fail(oe.getLocalizedMessage());
493         }
494 
495     }
496 
497 }