1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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.util.Binary64Field;
29 import org.hipparchus.util.FastMath;
30 import org.junit.jupiter.api.Assertions;
31 import org.junit.jupiter.api.BeforeEach;
32 import org.junit.jupiter.api.Test;
33 import org.orekit.Utils;
34 import org.orekit.attitudes.AttitudeProvider;
35 import org.orekit.attitudes.CelestialBodyPointed;
36 import org.orekit.attitudes.FrameAlignedProvider;
37 import org.orekit.attitudes.LofOffset;
38 import org.orekit.bodies.CelestialBodyFactory;
39 import org.orekit.errors.OrekitException;
40 import org.orekit.forces.AbstractForceModelTest;
41 import org.orekit.forces.maneuvers.ConstantThrustManeuver;
42 import org.orekit.frames.FramesFactory;
43 import org.orekit.frames.LOFType;
44 import org.orekit.orbits.*;
45 import org.orekit.propagation.*;
46 import org.orekit.propagation.numerical.FieldNumericalPropagator;
47 import org.orekit.propagation.numerical.NumericalPropagator;
48 import org.orekit.propagation.sampling.MultiSatStepHandler;
49 import org.orekit.time.AbsoluteDate;
50 import org.orekit.time.DateComponents;
51 import org.orekit.time.FieldAbsoluteDate;
52 import org.orekit.time.TimeComponents;
53 import org.orekit.time.TimeScalesFactory;
54 import org.orekit.utils.Constants;
55 import org.orekit.utils.PVCoordinates;
56
57 import java.util.Arrays;
58
59 public class PolynomialAccelerationModelTest extends AbstractForceModelTest {
60
61 private Orbit initialOrbit;
62
63 @Test
64 void testEquivalentInertialManeuver() {
65 final double delta = FastMath.toRadians(-7.4978);
66 final double alpha = FastMath.toRadians(351);
67 final Vector3D direction = new Vector3D(alpha, delta);
68 final double mass = 2500;
69 final double isp = Double.POSITIVE_INFINITY;
70 final double duration = 4000;
71 final double f = 400;
72
73 final AttitudeProvider maneuverLaw = new FrameAlignedProvider(new Rotation(direction, Vector3D.PLUS_I));
74 ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
75 duration, f, isp, Vector3D.PLUS_I);
76 final AttitudeProvider accelerationLaw = new FrameAlignedProvider(new Rotation(direction, Vector3D.PLUS_K));
77 final AccelerationModel accelerationModel = new PolynomialAccelerationModel("", AbsoluteDate.J2000_EPOCH, 0);
78 final ParametricAcceleration inertialAcceleration = new ParametricAcceleration(direction, true, accelerationModel);
79 Assertions.assertTrue(inertialAcceleration.dependsOnPositionOnly());
80 inertialAcceleration.getParametersDrivers().get(0).setValue(f / mass);
81 doTestEquivalentManeuver(mass, maneuverLaw, maneuver, accelerationLaw, inertialAcceleration, 1.0e-15);
82 }
83
84 @Test
85 void testEquivalentTangentialManeuver() {
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 commonLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
92 ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
93 duration, f, isp, Vector3D.PLUS_I);
94 final AccelerationModel accelerationModel = new PolynomialAccelerationModel("", null, 0);
95 final ParametricAcceleration lofAcceleration = new ParametricAcceleration(Vector3D.PLUS_I, false, accelerationModel);
96 Assertions.assertFalse(lofAcceleration.dependsOnPositionOnly());
97 lofAcceleration.getParametersDrivers().get(0).setValue(f / mass);
98 doTestEquivalentManeuver(mass, commonLaw, maneuver, commonLaw, lofAcceleration, 1.0e-15);
99 }
100
101 @Test
102 void testEquivalentTangentialOverriddenManeuver() {
103 final double mass = 2500;
104 final double isp = Double.POSITIVE_INFINITY;
105 final double duration = 4000;
106 final double f = 400;
107
108 final AttitudeProvider maneuverLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
109 ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
110 duration, f, isp, Vector3D.PLUS_I);
111 final AttitudeProvider accelerationLaw = new CelestialBodyPointed(initialOrbit.getFrame(),
112 CelestialBodyFactory.getSun(), Vector3D.PLUS_K,
113 Vector3D.PLUS_I, Vector3D.PLUS_K);
114 final AccelerationModel accelerationModel = new PolynomialAccelerationModel("prefix", null, 0);
115 final ParametricAcceleration lofAcceleration = new ParametricAcceleration(Vector3D.PLUS_I, maneuverLaw, accelerationModel);
116 lofAcceleration.getParametersDrivers().get(0).setValue(f / mass);
117 doTestEquivalentManeuver(mass, maneuverLaw, maneuver, accelerationLaw, lofAcceleration, 1.0e-15);
118 }
119
120 private void doTestEquivalentManeuver(final double mass,
121 final AttitudeProvider maneuverLaw,
122 final ConstantThrustManeuver maneuver,
123 final AttitudeProvider accelerationLaw,
124 final ParametricAcceleration parametricAcceleration,
125 final double positionTolerance)
126 {
127
128 SpacecraftState initialState = new SpacecraftState(initialOrbit,
129 maneuverLaw.getAttitude(initialOrbit,
130 initialOrbit.getDate(),
131 initialOrbit.getFrame())).withMass(mass);
132
133 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(initialOrbit, initialOrbit.getType());
134
135
136
137 AdaptiveStepsizeIntegrator integrator0 =
138 new DormandPrince853Integrator(0.001, 100, tolerance[0], tolerance[1]);
139 integrator0.setInitialStepSize(60);
140 final NumericalPropagator propagator0 = new NumericalPropagator(integrator0);
141 propagator0.setOrbitType(OrbitType.EQUINOCTIAL);
142 propagator0.setPositionAngleType(PositionAngleType.TRUE);
143 propagator0.setInitialState(initialState);
144 propagator0.setAttitudeProvider(maneuverLaw);
145 propagator0.addForceModel(maneuver);
146
147
148 AdaptiveStepsizeIntegrator integrator1 =
149 new DormandPrince853Integrator(0.001, 100, tolerance[0], tolerance[1]);
150 integrator1.setInitialStepSize(60);
151 final NumericalPropagator propagator1 = new NumericalPropagator(integrator1);
152 propagator1.setOrbitType(propagator0.getOrbitType());
153 propagator1.setPositionAngleType(propagator0.getPositionAngleType());
154 propagator1.setInitialState(initialState);
155 propagator1.setAttitudeProvider(accelerationLaw);
156 propagator1.addForceModel(parametricAcceleration);
157
158 MultiSatStepHandler handler = interpolators -> {
159 Vector3D p0 = interpolators.get(0).getCurrentState().getPosition();
160 Vector3D p1 = interpolators.get(1).getCurrentState().getPosition();
161 Assertions.assertEquals(0.0, Vector3D.distance(p0, p1), positionTolerance);
162 };
163 PropagatorsParallelizer parallelizer = new PropagatorsParallelizer(Arrays.asList(propagator0, propagator1),
164 handler);
165
166 parallelizer.propagate(initialOrbit.getDate(), initialOrbit.getDate().shiftedBy(1000.0));
167
168 }
169
170 @Test
171 void testEquivalentInertialManeuverField() {
172 final double delta = FastMath.toRadians(-7.4978);
173 final double alpha = FastMath.toRadians(351);
174 final Vector3D direction = new Vector3D(alpha, delta);
175 final double mass = 2500;
176 final double isp = Double.POSITIVE_INFINITY;
177 final double duration = 4000;
178 final double f = 400;
179
180 final AttitudeProvider maneuverLaw = new FrameAlignedProvider(new Rotation(direction, Vector3D.PLUS_I));
181 ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
182 duration, f, isp, Vector3D.PLUS_I);
183 final AttitudeProvider accelerationLaw = new FrameAlignedProvider(new Rotation(direction, Vector3D.PLUS_K));
184 final AccelerationModel accelerationModel = new PolynomialAccelerationModel("", AbsoluteDate.J2000_EPOCH, 0);
185 final ParametricAcceleration inertialAcceleration = new ParametricAcceleration(direction, true, accelerationModel);
186 inertialAcceleration.getParametersDrivers().get(0).setValue(f / mass);
187 doTestEquivalentManeuver(Binary64Field.getInstance(),
188 mass, maneuverLaw, maneuver, accelerationLaw, inertialAcceleration, 3.0e-9);
189 }
190
191 @Test
192 void testEquivalentTangentialManeuverField() {
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 commonLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
199 ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
200 duration, f, isp, Vector3D.PLUS_I);
201 final AccelerationModel accelerationModel = new PolynomialAccelerationModel("", null, 0);
202 final ParametricAcceleration lofAcceleration = new ParametricAcceleration(Vector3D.PLUS_I, false, accelerationModel);
203 lofAcceleration.getParametersDrivers().get(0).setValue(f / mass);
204 doTestEquivalentManeuver(Binary64Field.getInstance(),
205 mass, commonLaw, maneuver, commonLaw, lofAcceleration, 1.0e-15);
206 }
207
208 @Test
209 void testEquivalentTangentialOverriddenManeuverField() {
210 final double mass = 2500;
211 final double isp = Double.POSITIVE_INFINITY;
212 final double duration = 4000;
213 final double f = 400;
214
215 final AttitudeProvider maneuverLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
216 ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
217 duration, f, isp, Vector3D.PLUS_I);
218 final AttitudeProvider accelerationLaw = new CelestialBodyPointed(initialOrbit.getFrame(),
219 CelestialBodyFactory.getSun(), Vector3D.PLUS_K,
220 Vector3D.PLUS_I, Vector3D.PLUS_K);
221 final AccelerationModel accelerationModel = new PolynomialAccelerationModel("prefix", null, 0);
222 final ParametricAcceleration lofAcceleration = new ParametricAcceleration(Vector3D.PLUS_I, maneuverLaw, accelerationModel);
223 lofAcceleration.getParametersDrivers().get(0).setValue(f / mass);
224 doTestEquivalentManeuver(Binary64Field.getInstance(),
225 mass, maneuverLaw, maneuver, accelerationLaw, lofAcceleration, 1.0e-15);
226 }
227
228 private <T extends CalculusFieldElement<T>> void doTestEquivalentManeuver(final Field<T> field,
229 final double mass,
230 final AttitudeProvider maneuverLaw,
231 final ConstantThrustManeuver maneuver,
232 final AttitudeProvider accelerationLaw,
233 final ParametricAcceleration parametricAcceleration,
234 final double positionTolerance)
235 {
236
237 FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(field,
238 new SpacecraftState(initialOrbit,
239 maneuverLaw.getAttitude(initialOrbit,
240 initialOrbit.getDate(),
241 initialOrbit.getFrame())).withMass(mass));
242
243 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(10).getTolerances(
244 initialState.getOrbit(),
245 initialState.getOrbit().getType());
246
247
248
249 AdaptiveStepsizeFieldIntegrator<T> integrator0 =
250 new DormandPrince853FieldIntegrator<>(field, 0.001, 100, tolerance[0], tolerance[1]);
251 integrator0.setInitialStepSize(60);
252 final FieldNumericalPropagator<T> propagator0 = new FieldNumericalPropagator<>(field, integrator0);
253 propagator0.setOrbitType(OrbitType.EQUINOCTIAL);
254 propagator0.setPositionAngleType(PositionAngleType.TRUE);
255 propagator0.setInitialState(initialState);
256 propagator0.setAttitudeProvider(maneuverLaw);
257 propagator0.addForceModel(maneuver);
258 final FieldEphemerisGenerator<T> generator0 = propagator0.getEphemerisGenerator();
259 propagator0.propagate(initialState.getDate(), initialState.getDate().shiftedBy(1000.0));
260 FieldBoundedPropagator<T> ephemeris0 = generator0.getGeneratedEphemeris();
261
262
263 AdaptiveStepsizeFieldIntegrator<T> integrator1 =
264 new DormandPrince853FieldIntegrator<>(field, 0.001, 100, tolerance[0], tolerance[1]);
265 integrator1.setInitialStepSize(60);
266 final FieldNumericalPropagator<T> propagator1 = new FieldNumericalPropagator<>(field, integrator1);
267 propagator1.setOrbitType(propagator0.getOrbitType());
268 propagator1.setPositionAngleType(propagator0.getPositionAngleType());
269 propagator1.setInitialState(initialState);
270 propagator1.setAttitudeProvider(accelerationLaw);
271 propagator1.addForceModel(parametricAcceleration);
272 final FieldEphemerisGenerator<T> generator1 = propagator1.getEphemerisGenerator();
273 propagator1.propagate(initialState.getDate(), initialState.getDate().shiftedBy(1000.0));
274 FieldBoundedPropagator<T> ephemeris1 = generator1.getGeneratedEphemeris();
275
276 for (double dt = 1; dt < 999; dt += 10) {
277 FieldAbsoluteDate<T> t = initialState.getDate().shiftedBy(dt);
278 FieldVector3D<T> p0 = ephemeris0.propagate(t).getPosition();
279 FieldVector3D<T> p1 = ephemeris1.propagate(t).getPosition();
280 Assertions.assertEquals(0, FieldVector3D.distance(p0, p1).getReal(), positionTolerance);
281 }
282
283 }
284
285 @Test
286 void testParameterDerivative() {
287
288
289 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
290 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
291 final SpacecraftState state =
292 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
293 FramesFactory.getGCRF(),
294 new AbsoluteDate(2005, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
295 Constants.EIGEN5C_EARTH_MU));
296
297 final AccelerationModel accelerationModel = new PolynomialAccelerationModel("ppa",
298 state.getDate().shiftedBy(-2.0),
299 3);
300 final ParametricAcceleration ppa = new ParametricAcceleration(Vector3D.PLUS_K, false, accelerationModel);
301
302 ppa.init(state, state.getDate().shiftedBy(3600.0));
303 ppa.getParametersDrivers().get(0).setValue(0.00001);
304 ppa.getParametersDrivers().get(1).setValue(0.00002);
305 ppa.getParametersDrivers().get(2).setValue(0.00003);
306 ppa.getParametersDrivers().get(3).setValue(0.00004);
307
308 checkParameterDerivative(state, ppa, "ppa[0]", 1.0e-3, 7.0e-12);
309 checkParameterDerivative(state, ppa, "ppa[1]", 1.0e-3, 9.0e-13);
310 checkParameterDerivative(state, ppa, "ppa[2]", 1.0e-3, 2.0e-12);
311 checkParameterDerivative(state, ppa, "ppa[3]", 1.0e-3, 3.0e-12);
312
313 }
314
315 @BeforeEach
316 public void setUp() {
317 try {
318 Utils.setDataRoot("regular-data");
319
320 final double a = 24396159;
321 final double e = 0.72831215;
322 final double i = FastMath.toRadians(7);
323 final double omega = FastMath.toRadians(180);
324 final double OMEGA = FastMath.toRadians(261);
325 final double lv = 0;
326
327 final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 01, 01),
328 new TimeComponents(23, 30, 00.000),
329 TimeScalesFactory.getUTC());
330 initialOrbit =
331 new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
332 FramesFactory.getEME2000(), initDate, Constants.EIGEN5C_EARTH_MU);
333 } catch (OrekitException oe) {
334 Assertions.fail(oe.getLocalizedMessage());
335 }
336
337 }
338
339 }