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