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()),
132 mass);
133
134 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(initialOrbit, initialOrbit.getType());
135
136
137
138 AdaptiveStepsizeIntegrator integrator0 =
139 new DormandPrince853Integrator(0.001, 100, tolerance[0], tolerance[1]);
140 integrator0.setInitialStepSize(60);
141 final NumericalPropagator propagator0 = new NumericalPropagator(integrator0);
142 propagator0.setOrbitType(OrbitType.EQUINOCTIAL);
143 propagator0.setPositionAngleType(PositionAngleType.TRUE);
144 propagator0.setInitialState(initialState);
145 propagator0.setAttitudeProvider(maneuverLaw);
146 propagator0.addForceModel(maneuver);
147
148
149 AdaptiveStepsizeIntegrator integrator1 =
150 new DormandPrince853Integrator(0.001, 100, tolerance[0], tolerance[1]);
151 integrator1.setInitialStepSize(60);
152 final NumericalPropagator propagator1 = new NumericalPropagator(integrator1);
153 propagator1.setOrbitType(propagator0.getOrbitType());
154 propagator1.setPositionAngleType(propagator0.getPositionAngleType());
155 propagator1.setInitialState(initialState);
156 propagator1.setAttitudeProvider(accelerationLaw);
157 propagator1.addForceModel(parametricAcceleration);
158
159 MultiSatStepHandler handler = interpolators -> {
160 Vector3D p0 = interpolators.get(0).getCurrentState().getPosition();
161 Vector3D p1 = interpolators.get(1).getCurrentState().getPosition();
162 Assertions.assertEquals(0.0, Vector3D.distance(p0, p1), positionTolerance);
163 };
164 PropagatorsParallelizer parallelizer = new PropagatorsParallelizer(Arrays.asList(propagator0, propagator1),
165 handler);
166
167 parallelizer.propagate(initialOrbit.getDate(), initialOrbit.getDate().shiftedBy(1000.0));
168
169 }
170
171 @Test
172 void testEquivalentInertialManeuverField() {
173 final double delta = FastMath.toRadians(-7.4978);
174 final double alpha = FastMath.toRadians(351);
175 final Vector3D direction = new Vector3D(alpha, delta);
176 final double mass = 2500;
177 final double isp = Double.POSITIVE_INFINITY;
178 final double duration = 4000;
179 final double f = 400;
180
181 final AttitudeProvider maneuverLaw = new FrameAlignedProvider(new Rotation(direction, Vector3D.PLUS_I));
182 ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
183 duration, f, isp, Vector3D.PLUS_I);
184 final AttitudeProvider accelerationLaw = new FrameAlignedProvider(new Rotation(direction, Vector3D.PLUS_K));
185 final AccelerationModel accelerationModel = new PolynomialAccelerationModel("", AbsoluteDate.J2000_EPOCH, 0);
186 final ParametricAcceleration inertialAcceleration = new ParametricAcceleration(direction, true, accelerationModel);
187 inertialAcceleration.getParametersDrivers().get(0).setValue(f / mass);
188 doTestEquivalentManeuver(Binary64Field.getInstance(),
189 mass, maneuverLaw, maneuver, accelerationLaw, inertialAcceleration, 3.0e-9);
190 }
191
192 @Test
193 void testEquivalentTangentialManeuverField() {
194 final double mass = 2500;
195 final double isp = Double.POSITIVE_INFINITY;
196 final double duration = 4000;
197 final double f = 400;
198
199 final AttitudeProvider commonLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
200 ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
201 duration, f, isp, Vector3D.PLUS_I);
202 final AccelerationModel accelerationModel = new PolynomialAccelerationModel("", null, 0);
203 final ParametricAcceleration lofAcceleration = new ParametricAcceleration(Vector3D.PLUS_I, false, accelerationModel);
204 lofAcceleration.getParametersDrivers().get(0).setValue(f / mass);
205 doTestEquivalentManeuver(Binary64Field.getInstance(),
206 mass, commonLaw, maneuver, commonLaw, lofAcceleration, 1.0e-15);
207 }
208
209 @Test
210 void testEquivalentTangentialOverriddenManeuverField() {
211 final double mass = 2500;
212 final double isp = Double.POSITIVE_INFINITY;
213 final double duration = 4000;
214 final double f = 400;
215
216 final AttitudeProvider maneuverLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
217 ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialOrbit.getDate().shiftedBy(-10.0),
218 duration, f, isp, Vector3D.PLUS_I);
219 final AttitudeProvider accelerationLaw = new CelestialBodyPointed(initialOrbit.getFrame(),
220 CelestialBodyFactory.getSun(), Vector3D.PLUS_K,
221 Vector3D.PLUS_I, Vector3D.PLUS_K);
222 final AccelerationModel accelerationModel = new PolynomialAccelerationModel("prefix", null, 0);
223 final ParametricAcceleration lofAcceleration = new ParametricAcceleration(Vector3D.PLUS_I, maneuverLaw, accelerationModel);
224 lofAcceleration.getParametersDrivers().get(0).setValue(f / mass);
225 doTestEquivalentManeuver(Binary64Field.getInstance(),
226 mass, maneuverLaw, maneuver, accelerationLaw, lofAcceleration, 1.0e-15);
227 }
228
229 private <T extends CalculusFieldElement<T>> void doTestEquivalentManeuver(final Field<T> field,
230 final double mass,
231 final AttitudeProvider maneuverLaw,
232 final ConstantThrustManeuver maneuver,
233 final AttitudeProvider accelerationLaw,
234 final ParametricAcceleration parametricAcceleration,
235 final double positionTolerance)
236 {
237
238 FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(field,
239 new SpacecraftState(initialOrbit,
240 maneuverLaw.getAttitude(initialOrbit,
241 initialOrbit.getDate(),
242 initialOrbit.getFrame()),
243 mass));
244
245 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(10).getTolerances(
246 initialState.getOrbit(),
247 initialState.getOrbit().getType());
248
249
250
251 AdaptiveStepsizeFieldIntegrator<T> integrator0 =
252 new DormandPrince853FieldIntegrator<>(field, 0.001, 100, tolerance[0], tolerance[1]);
253 integrator0.setInitialStepSize(60);
254 final FieldNumericalPropagator<T> propagator0 = new FieldNumericalPropagator<>(field, integrator0);
255 propagator0.setOrbitType(OrbitType.EQUINOCTIAL);
256 propagator0.setPositionAngleType(PositionAngleType.TRUE);
257 propagator0.setInitialState(initialState);
258 propagator0.setAttitudeProvider(maneuverLaw);
259 propagator0.addForceModel(maneuver);
260 final FieldEphemerisGenerator<T> generator0 = propagator0.getEphemerisGenerator();
261 propagator0.propagate(initialState.getDate(), initialState.getDate().shiftedBy(1000.0));
262 FieldBoundedPropagator<T> ephemeris0 = generator0.getGeneratedEphemeris();
263
264
265 AdaptiveStepsizeFieldIntegrator<T> integrator1 =
266 new DormandPrince853FieldIntegrator<>(field, 0.001, 100, tolerance[0], tolerance[1]);
267 integrator1.setInitialStepSize(60);
268 final FieldNumericalPropagator<T> propagator1 = new FieldNumericalPropagator<>(field, integrator1);
269 propagator1.setOrbitType(propagator0.getOrbitType());
270 propagator1.setPositionAngleType(propagator0.getPositionAngleType());
271 propagator1.setInitialState(initialState);
272 propagator1.setAttitudeProvider(accelerationLaw);
273 propagator1.addForceModel(parametricAcceleration);
274 final FieldEphemerisGenerator<T> generator1 = propagator1.getEphemerisGenerator();
275 propagator1.propagate(initialState.getDate(), initialState.getDate().shiftedBy(1000.0));
276 FieldBoundedPropagator<T> ephemeris1 = generator1.getGeneratedEphemeris();
277
278 for (double dt = 1; dt < 999; dt += 10) {
279 FieldAbsoluteDate<T> t = initialState.getDate().shiftedBy(dt);
280 FieldVector3D<T> p0 = ephemeris0.propagate(t).getPosition();
281 FieldVector3D<T> p1 = ephemeris1.propagate(t).getPosition();
282 Assertions.assertEquals(0, FieldVector3D.distance(p0, p1).getReal(), positionTolerance);
283 }
284
285 }
286
287 @Test
288 void testParameterDerivative() {
289
290
291 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
292 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
293 final SpacecraftState state =
294 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
295 FramesFactory.getGCRF(),
296 new AbsoluteDate(2005, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
297 Constants.EIGEN5C_EARTH_MU));
298
299 final AccelerationModel accelerationModel = new PolynomialAccelerationModel("ppa",
300 state.getDate().shiftedBy(-2.0),
301 3);
302 final ParametricAcceleration ppa = new ParametricAcceleration(Vector3D.PLUS_K, false, accelerationModel);
303
304 ppa.init(state, state.getDate().shiftedBy(3600.0));
305 ppa.getParametersDrivers().get(0).setValue(0.00001);
306 ppa.getParametersDrivers().get(1).setValue(0.00002);
307 ppa.getParametersDrivers().get(2).setValue(0.00003);
308 ppa.getParametersDrivers().get(3).setValue(0.00004);
309
310 checkParameterDerivative(state, ppa, "ppa[0]", 1.0e-3, 7.0e-12);
311 checkParameterDerivative(state, ppa, "ppa[1]", 1.0e-3, 9.0e-13);
312 checkParameterDerivative(state, ppa, "ppa[2]", 1.0e-3, 2.0e-12);
313 checkParameterDerivative(state, ppa, "ppa[3]", 1.0e-3, 3.0e-12);
314
315 }
316
317 @BeforeEach
318 public void setUp() {
319 try {
320 Utils.setDataRoot("regular-data");
321
322 final double a = 24396159;
323 final double e = 0.72831215;
324 final double i = FastMath.toRadians(7);
325 final double omega = FastMath.toRadians(180);
326 final double OMEGA = FastMath.toRadians(261);
327 final double lv = 0;
328
329 final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 01, 01),
330 new TimeComponents(23, 30, 00.000),
331 TimeScalesFactory.getUTC());
332 initialOrbit =
333 new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
334 FramesFactory.getEME2000(), initDate, Constants.EIGEN5C_EARTH_MU);
335 } catch (OrekitException oe) {
336 Assertions.fail(oe.getLocalizedMessage());
337 }
338
339 }
340
341 }