1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces.maneuvers.propulsion;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.Field;
21 import org.hipparchus.analysis.polynomials.PolynomialFunction;
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.mockito.Mockito;
34 import org.orekit.Utils;
35 import org.orekit.attitudes.AttitudeProvider;
36 import org.orekit.attitudes.FrameAlignedProvider;
37 import org.orekit.forces.maneuvers.Control3DVectorCostType;
38 import org.orekit.forces.maneuvers.Maneuver;
39 import org.orekit.forces.maneuvers.trigger.DateBasedManeuverTriggers;
40 import org.orekit.frames.FramesFactory;
41 import org.orekit.orbits.FieldKeplerianOrbit;
42 import org.orekit.orbits.FieldOrbit;
43 import org.orekit.orbits.KeplerianOrbit;
44 import org.orekit.orbits.Orbit;
45 import org.orekit.orbits.PositionAngleType;
46 import org.orekit.propagation.FieldSpacecraftState;
47 import org.orekit.propagation.SpacecraftState;
48 import org.orekit.propagation.ToleranceProvider;
49 import org.orekit.propagation.numerical.FieldNumericalPropagator;
50 import org.orekit.propagation.numerical.NumericalPropagator;
51 import org.orekit.time.AbsoluteDate;
52 import org.orekit.time.DateComponents;
53 import org.orekit.time.FieldAbsoluteDate;
54 import org.orekit.time.TimeComponents;
55 import org.orekit.time.TimeScalesFactory;
56 import org.orekit.utils.Constants;
57 import org.orekit.utils.ParameterDriver;
58 import org.orekit.utils.TimeSpanMap;
59
60 import java.util.List;
61
62 class ProfileThrustPropulsionModelTest {
63
64 @Test
65 void testRoughBehaviour() {
66 doRoughBehaviour( 1.0, 2008.017, 28968115.974);
67 doRoughBehaviour( 10.0, 2009.229, 28950587.132);
68 doRoughBehaviour(100.0, 2021.350, 28777772.266);
69 }
70
71 @Test
72 void testRoughBehaviourField() {
73 doRoughBehaviourField(Binary64Field.getInstance(), 1.0, 2008.017, 28968115.974);
74 doRoughBehaviourField(Binary64Field.getInstance(), 10.0, 2009.229, 28950587.132);
75 doRoughBehaviourField(Binary64Field.getInstance(), 100.0, 2021.350, 28777772.266);
76 }
77
78 private void doRoughBehaviour(final double rampDuration, final double expectedM, final double expectedA) {
79 final double isp = 318;
80 final double mass = 2500;
81 final double a = 24396159;
82 final double e = 0.72831215;
83 final double i = FastMath.toRadians(7);
84 final double omega = FastMath.toRadians(180);
85 final double OMEGA = FastMath.toRadians(261);
86 final double lv = 0;
87
88 final double duration = 3653.99;
89 final double f = 420;
90 final double delta = FastMath.toRadians(-7.4978);
91 final double alpha = FastMath.toRadians(351);
92 final AttitudeProvider law = new FrameAlignedProvider(new Rotation(new Vector3D(alpha, delta),
93 Vector3D.PLUS_I));
94
95 final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 1, 1),
96 new TimeComponents(23, 30, 00.000),
97 TimeScalesFactory.getUTC());
98 final PositionAngleType positionAngleType = PositionAngleType.TRUE;
99 final Orbit initOrbit =
100 new KeplerianOrbit(a, e, i, omega, OMEGA, lv, positionAngleType,
101 FramesFactory.getEME2000(), initDate, Constants.EIGEN5C_EARTH_MU);
102 final SpacecraftState initialState =
103 new SpacecraftState(initOrbit, law.getAttitude(initOrbit, initOrbit.getDate(), initOrbit.getFrame()), mass);
104
105 final AbsoluteDate fireDate = new AbsoluteDate(new DateComponents(2004, 1, 2),
106 new TimeComponents(4, 15, 34.080),
107 TimeScalesFactory.getUTC());
108
109 final TimeSpanMap<PolynomialThrustSegment> profile = new TimeSpanMap<>(null);
110
111 final AbsoluteDate t0 = fireDate;
112 profile.addValidAfter(new PolynomialThrustSegment(t0,
113 new PolynomialFunction(0.0, f / rampDuration),
114 new PolynomialFunction(0.0),
115 new PolynomialFunction(0.0)),
116 t0, false);
117
118 final AbsoluteDate t1 = fireDate.shiftedBy(rampDuration);
119 profile.addValidAfter(new PolynomialThrustSegment(t1,
120 new PolynomialFunction(f),
121 new PolynomialFunction(0.0),
122 new PolynomialFunction(0.0)),
123 t1, false);
124
125 final AbsoluteDate t2 = fireDate.shiftedBy(duration - rampDuration);
126 profile.addValidAfter(new PolynomialThrustSegment(t2,
127 new PolynomialFunction(f, -f / rampDuration),
128 new PolynomialFunction(0.0),
129 new PolynomialFunction(0.0)),
130 t2, false);
131
132 final AbsoluteDate t3 = fireDate.shiftedBy(duration);
133 profile.addValidAfter(new PolynomialThrustSegment(t3,
134 new PolynomialFunction(0.0),
135 new PolynomialFunction(0.0),
136 new PolynomialFunction(0.0)),
137 t3, false);
138 final PropulsionModel propulsionModel = ProfileThrustPropulsionModel.of(profile, isp, Control3DVectorCostType.TWO_NORM, "ABM");
139
140 Assertions.assertEquals("ABM", propulsionModel.getName());
141 Assertions.assertEquals(0.0, thrust(initialState, t0.shiftedBy(-0.001), propulsionModel), 1.0e-8);
142 Assertions.assertEquals(0.5 * f, thrust(initialState, t0.shiftedBy(+0.5 * rampDuration), propulsionModel), 1.0e-8);
143 Assertions.assertEquals(f, thrust(initialState, t1.shiftedBy(+0.001), propulsionModel), 1.0e-8);
144 Assertions.assertEquals(0.5 * f, thrust(initialState, t3.shiftedBy(-0.5 * rampDuration), propulsionModel), 1.0e-8);
145 Assertions.assertEquals(0.0, thrust(initialState, t3.shiftedBy(+0.001), propulsionModel), 1.0e-8);
146
147 final Maneuver maneuver = new Maneuver(null, new DateBasedManeuverTriggers(fireDate, duration), propulsionModel);
148
149 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(1e-6).getTolerances(initOrbit, initOrbit.getType());
150 AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(1.0e-6, rampDuration, tolerance[0], tolerance[1]);
151 integrator.setInitialStepSize(0.1);
152 final NumericalPropagator propagator = new NumericalPropagator(integrator);
153 propagator.setOrbitType(initOrbit.getType());
154 propagator.setInitialState(initialState);
155 propagator.setAttitudeProvider(law);
156 propagator.addForceModel(maneuver);
157 final SpacecraftState finalorb = propagator.propagate(fireDate.shiftedBy(3800));
158
159 Assertions.assertEquals(expectedM, finalorb.getMass(), 0.001);
160 Assertions.assertEquals(expectedA, finalorb.getOrbit().getA(), 0.002);
161
162 }
163
164 private <T extends CalculusFieldElement<T>> void doRoughBehaviourField(final Field<T> field, final double rampDuration,
165 final double expectedM, final double expectedA) {
166 final double isp = 318;
167 final T zero = field.getZero();
168 final T mass = zero.newInstance(2500);
169 final T a = zero.newInstance(24396159);
170 final T e = zero.newInstance(0.72831215);
171 final T i = FastMath.toRadians(zero.newInstance(7));
172 final T omega = FastMath.toRadians(zero.newInstance(180));
173 final T OMEGA = FastMath.toRadians(zero.newInstance(261));
174 final T lv = zero;
175
176 final double duration = 3653.99;
177 final double f = 420;
178 final double delta = FastMath.toRadians(-7.4978);
179 final double alpha = FastMath.toRadians(351);
180 final AttitudeProvider law = new FrameAlignedProvider(new Rotation(new Vector3D(alpha, delta),
181 Vector3D.PLUS_I));
182
183 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field,
184 new DateComponents(2004, 1, 1),
185 new TimeComponents(23, 30, 00.000),
186 TimeScalesFactory.getUTC());
187 final PositionAngleType positionAngleType = PositionAngleType.TRUE;
188 final FieldOrbit<T> initOrbit =
189 new FieldKeplerianOrbit<>(a, e, i, omega, OMEGA, lv, positionAngleType,
190 FramesFactory.getEME2000(), initDate,
191 zero.newInstance(Constants.EIGEN5C_EARTH_MU));
192 final FieldSpacecraftState<T> initialState =
193 new FieldSpacecraftState<>(initOrbit, law.getAttitude(initOrbit, initOrbit.getDate(), initOrbit.getFrame()), mass);
194
195 final FieldAbsoluteDate<T> fireDate = new FieldAbsoluteDate<>(field,
196 new DateComponents(2004, 1, 2),
197 new TimeComponents(4, 15, 34.080),
198 TimeScalesFactory.getUTC());
199
200 final TimeSpanMap<PolynomialThrustSegment> profile = new TimeSpanMap<>(null);
201
202 final AbsoluteDate t0 = fireDate.toAbsoluteDate();
203 profile.addValidAfter(new PolynomialThrustSegment(t0,
204 new PolynomialFunction(0.0, f / rampDuration),
205 new PolynomialFunction(0.0),
206 new PolynomialFunction(0.0)),
207 t0, false);
208
209 final AbsoluteDate t1 = fireDate.shiftedBy(rampDuration).toAbsoluteDate();
210 profile.addValidAfter(new PolynomialThrustSegment(t1,
211 new PolynomialFunction(f),
212 new PolynomialFunction(0.0),
213 new PolynomialFunction(0.0)),
214 t1, false);
215
216 final AbsoluteDate t2 = fireDate.shiftedBy(duration - rampDuration).toAbsoluteDate();
217 profile.addValidAfter(new PolynomialThrustSegment(t2,
218 new PolynomialFunction(f, -f / rampDuration),
219 new PolynomialFunction(0.0),
220 new PolynomialFunction(0.0)),
221 t2, false);
222
223 final AbsoluteDate t3 = fireDate.shiftedBy(duration).toAbsoluteDate();
224 profile.addValidAfter(new PolynomialThrustSegment(t3,
225 new PolynomialFunction(0.0),
226 new PolynomialFunction(0.0),
227 new PolynomialFunction(0.0)),
228 t3, false);
229 final PropulsionModel propulsionModel = ProfileThrustPropulsionModel.of(profile, isp, Control3DVectorCostType.TWO_NORM, "ABM");
230
231 Assertions.assertEquals("ABM", propulsionModel.getName());
232 Assertions.assertEquals(0.0, thrust(initialState, fireDate.shiftedBy(-0.001), propulsionModel).getReal(), 1.0e-8);
233 Assertions.assertEquals(0.5 * f, thrust(initialState, fireDate.shiftedBy(0.5 * rampDuration), propulsionModel).getReal(), 1.0e-8);
234 Assertions.assertEquals(f, thrust(initialState, fireDate.shiftedBy( rampDuration + 0.001), propulsionModel).getReal(), 1.0e-8);
235 Assertions.assertEquals(0.5 * f, thrust(initialState, fireDate.shiftedBy(duration - 0.5 * rampDuration), propulsionModel).getReal(), 1.0e-8);
236 Assertions.assertEquals(0.0, thrust(initialState, fireDate.shiftedBy(duration + 0.001), propulsionModel).getReal(), 1.0e-8);
237
238 final Maneuver maneuver = new Maneuver(null, new DateBasedManeuverTriggers(fireDate.toAbsoluteDate(), duration), propulsionModel);
239
240 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(1e-6).getTolerances(initOrbit, initOrbit.getType());
241 AdaptiveStepsizeFieldIntegrator<T> integrator =
242 new DormandPrince853FieldIntegrator<>(field, 1.0e-6, rampDuration, tolerance[0], tolerance[1]);
243 integrator.setInitialStepSize(0.1);
244 final FieldNumericalPropagator<T> propagator = new FieldNumericalPropagator<>(field, integrator);
245 propagator.setOrbitType(initOrbit.getType());
246 propagator.setPositionAngleType(positionAngleType);
247 propagator.setInitialState(initialState);
248 propagator.setAttitudeProvider(law);
249 propagator.addForceModel(maneuver);
250 final FieldSpacecraftState<T> finalorb = propagator.propagate(fireDate.shiftedBy(3800));
251
252 Assertions.assertEquals(expectedM, finalorb.getMass().getReal(), 0.001);
253 Assertions.assertEquals(expectedA, finalorb.getOrbit().getA().getReal(), 0.002);
254
255 }
256
257 private double thrust(final SpacecraftState initialState, final AbsoluteDate targetDate, final PropulsionModel propulsionModel) {
258 final SpacecraftState state = initialState.shiftedBy(targetDate.durationFrom(initialState.getDate()));
259 return state.getMass() * propulsionModel.getAcceleration(state, state.getAttitude(), null).getNorm();
260 }
261
262 private <T extends CalculusFieldElement<T>> T thrust(final FieldSpacecraftState<T> initialState,
263 final FieldAbsoluteDate<T> targetDate,
264 final PropulsionModel propulsionModel) {
265 final FieldSpacecraftState<T> state = initialState.shiftedBy(targetDate.durationFrom(initialState.getDate()));
266 return state.getMass().multiply(propulsionModel.getAcceleration(state, state.getAttitude(), null).getNorm());
267 }
268
269 @Test
270 void testGetThrustVector() {
271
272 final ProfileThrustPropulsionModel propulsionModel = new ProfileThrustPropulsionModel(new TimeSpanMap<>(null),
273 1., "");
274 final AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
275 final SpacecraftState mockedState = Mockito.mock(SpacecraftState.class);
276 Mockito.when(mockedState.getDate()).thenReturn(date);
277
278 final Vector3D actualThrust = propulsionModel.getThrustVector(mockedState);
279
280 Assertions.assertEquals(propulsionModel.getThrustVector(mockedState, new double[0]), actualThrust);
281 }
282
283 @Test
284 void testGetFlowRate() {
285
286 final ProfileThrustPropulsionModel propulsionModel = new ProfileThrustPropulsionModel(new TimeSpanMap<>(null),
287 1., "");
288 final AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
289 final SpacecraftState mockedState = Mockito.mock(SpacecraftState.class);
290 Mockito.when(mockedState.getDate()).thenReturn(date);
291
292 final double actualRate = propulsionModel.getFlowRate(mockedState);
293
294 Assertions.assertEquals(propulsionModel.getFlowRate(mockedState, new double[0]), actualRate);
295 }
296
297 @Test
298 void testGetControl3DVectorCostTypeDefault() {
299
300 final ProfileThrustPropulsionModel propulsionModel = new ProfileThrustPropulsionModel(new TimeSpanMap<>(null),
301 1., "");
302
303 final Control3DVectorCostType actualCostType = propulsionModel.getControl3DVectorCostType();
304
305 Assertions.assertEquals(Control3DVectorCostType.TWO_NORM, actualCostType);
306 }
307
308 @Test
309 void testGetControl3DVectorCostType() {
310
311 final Control3DVectorCostType expectedCostType = Control3DVectorCostType.ONE_NORM;
312 final ProfileThrustPropulsionModel propulsionModel = new ProfileThrustPropulsionModel(new TimeSpanMap<>(null),
313 1., expectedCostType, "");
314
315 final Control3DVectorCostType actualCostType = propulsionModel.getControl3DVectorCostType();
316
317 Assertions.assertEquals(expectedCostType, actualCostType);
318 }
319
320 @Test
321 void testGetParametersDrivers() {
322
323 final ParameterDriver driver = Mockito.mock(ParameterDriver.class);
324 final String expectedName = "a";
325 Mockito.when(driver.getName()).thenReturn(expectedName);
326 final ProfileThrustPropulsionModel propulsionModel = new ProfileThrustPropulsionModel(new TimeSpanMap<>(null),
327 1., Control3DVectorCostType.NONE, "");
328
329 final List<ParameterDriver> propulsionDrivers = propulsionModel.getParametersDrivers();
330
331 Assertions.assertEquals(0, propulsionDrivers.size());
332 }
333
334 @BeforeEach
335 public void setUp() {
336 Utils.setDataRoot("regular-data");
337 }
338
339 }