1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces.radiation;
18
19 import org.hipparchus.Field;
20 import org.hipparchus.analysis.differentiation.DSFactory;
21 import org.hipparchus.analysis.differentiation.DerivativeStructure;
22 import org.hipparchus.analysis.differentiation.Gradient;
23 import org.hipparchus.geometry.euclidean.threed.Vector3D;
24 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
25 import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaFieldIntegrator;
26 import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator;
27 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
28 import org.hipparchus.util.FastMath;
29 import org.junit.jupiter.api.Assertions;
30 import org.junit.jupiter.api.BeforeEach;
31 import org.junit.jupiter.api.Test;
32 import org.orekit.Utils;
33 import org.orekit.bodies.CelestialBodyFactory;
34 import org.orekit.forces.AbstractForceModelTest;
35 import org.orekit.frames.Frame;
36 import org.orekit.frames.FramesFactory;
37 import org.orekit.orbits.CartesianOrbit;
38 import org.orekit.orbits.FieldKeplerianOrbit;
39 import org.orekit.orbits.KeplerianOrbit;
40 import org.orekit.orbits.Orbit;
41 import org.orekit.orbits.OrbitType;
42 import org.orekit.orbits.PositionAngleType;
43 import org.orekit.propagation.FieldSpacecraftState;
44 import org.orekit.propagation.SpacecraftState;
45 import org.orekit.propagation.ToleranceProvider;
46 import org.orekit.propagation.analytical.tle.TLE;
47 import org.orekit.propagation.analytical.tle.TLEPropagator;
48 import org.orekit.propagation.numerical.FieldNumericalPropagator;
49 import org.orekit.propagation.numerical.NumericalPropagator;
50 import org.orekit.propagation.sampling.OrekitFixedStepHandler;
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.ExtendedPositionProvider;
58 import org.orekit.utils.PVCoordinates;
59
60
61
62
63
64
65
66 class KnockeRediffusedForceModelTest extends AbstractForceModelTest{
67
68 @BeforeEach
69 public void setUp() {
70 Utils.setDataRoot("regular-data");
71 }
72
73
74 @Test
75 void testJacobianVsFiniteDifferences() {
76
77
78 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 3, 1),
79 new TimeComponents(13, 59, 27.816),
80 TimeScalesFactory.getUTC());
81 double i = FastMath.toRadians(98.7);
82 double omega = FastMath.toRadians(93.0);
83 double OMEGA = FastMath.toRadians(15.0 * 22.5);
84 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
85 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
86 Constants.EIGEN5C_EARTH_MU);
87
88
89 final ExtendedPositionProvider sun = CelestialBodyFactory.getSun();
90
91
92 final RadiationSensitive radiationSensitive = new IsotropicRadiationSingleCoefficient(1, 1.5);
93
94
95 final KnockeRediffusedForceModel forceModel = new KnockeRediffusedForceModel(sun,
96 radiationSensitive,
97 Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
98 FastMath.toRadians(30));
99
100 SpacecraftState state = new SpacecraftState(orbit,
101 Utils.defaultLaw().getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
102 checkStateJacobianVsFiniteDifferences(state, forceModel, Utils.defaultLaw(), 1.0, 5.5e-9, false);
103
104 }
105
106 @Test
107 void testParameterIsotropicSingle() {
108
109 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
110 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
111 final SpacecraftState state =
112 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
113 FramesFactory.getGCRF(),
114 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
115 Constants.EIGEN5C_EARTH_MU));
116
117
118 final ExtendedPositionProvider sun = CelestialBodyFactory.getSun();
119
120
121 final RadiationSensitive radiationSensitive = new IsotropicRadiationSingleCoefficient(1, 1.5);
122
123
124 final KnockeRediffusedForceModel forceModel = new KnockeRediffusedForceModel(sun,
125 radiationSensitive,
126 Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
127 FastMath.toRadians(30));
128
129 checkParameterDerivative(state, forceModel, RadiationSensitive.REFLECTION_COEFFICIENT, 0.25, 1.8e-16);
130
131 }
132
133 @Test
134 void testGlobalStateJacobianIsotropicSingle()
135 {
136
137
138 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 3, 1),
139 new TimeComponents(13, 59, 27.816),
140 TimeScalesFactory.getUTC());
141 double i = FastMath.toRadians(98.7);
142 double omega = FastMath.toRadians(93.0);
143 double OMEGA = FastMath.toRadians(15.0 * 22.5);
144 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
145 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
146 Constants.EIGEN5C_EARTH_MU);
147 OrbitType integrationType = OrbitType.CARTESIAN;
148 double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(0.01).getTolerances(orbit, integrationType);
149
150 NumericalPropagator propagator =
151 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
152 tolerances[0], tolerances[1]));
153 propagator.setOrbitType(integrationType);
154
155
156 final ExtendedPositionProvider sun = CelestialBodyFactory.getSun();
157
158
159 final RadiationSensitive radiationSensitive = new IsotropicRadiationSingleCoefficient(1, 1.5);
160
161
162 final KnockeRediffusedForceModel forceModel = new KnockeRediffusedForceModel(sun,
163 radiationSensitive,
164 Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
165 FastMath.toRadians(30));
166 propagator.addForceModel(forceModel);
167 SpacecraftState state0 = new SpacecraftState(orbit);
168
169 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
170 1e3, tolerances[0], 3.2e-8);
171
172 }
173
174 @Test
175 void testRealField() {
176
177
178
179 DSFactory factory = new DSFactory(6, 4);
180 DerivativeStructure a_0 = factory.variable(0, 7e6);
181 DerivativeStructure e_0 = factory.variable(1, 0.01);
182 DerivativeStructure i_0 = factory.variable(2, 1.2);
183 DerivativeStructure R_0 = factory.variable(3, 0.7);
184 DerivativeStructure O_0 = factory.variable(4, 0.5);
185 DerivativeStructure n_0 = factory.variable(5, 0.1);
186
187 Field<DerivativeStructure> field = a_0.getField();
188 DerivativeStructure zero = field.getZero();
189
190
191 FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
192
193
194 Frame EME = FramesFactory.getEME2000();
195
196
197 FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
198 PositionAngleType.MEAN,
199 EME,
200 J2000,
201 zero.add(Constants.EIGEN5C_EARTH_MU));
202
203
204 FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
205 SpacecraftState iSR = initialState.toSpacecraftState();
206
207
208 ClassicalRungeKuttaFieldIntegrator<DerivativeStructure> integrator =
209 new ClassicalRungeKuttaFieldIntegrator<>(field, zero.add(6));
210 ClassicalRungeKuttaIntegrator RIntegrator =
211 new ClassicalRungeKuttaIntegrator(6);
212 OrbitType type = OrbitType.KEPLERIAN;
213
214
215 FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
216 FNP.setOrbitType(type);
217 FNP.setInitialState(initialState);
218
219 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
220 NP.setOrbitType(type);
221 NP.setInitialState(iSR);
222
223
224 final ExtendedPositionProvider sun = CelestialBodyFactory.getSun();
225
226
227 final RadiationSensitive radiationSensitive = new IsotropicRadiationSingleCoefficient(1, 1.5);
228
229
230 final KnockeRediffusedForceModel forceModel = new KnockeRediffusedForceModel(sun,
231 radiationSensitive,
232 Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
233 FastMath.toRadians(30));
234
235 FNP.addForceModel(forceModel);
236 NP.addForceModel(forceModel);
237
238
239 checkRealFieldPropagation(FKO, PositionAngleType.MEAN, 300., NP, FNP,
240 1.0e-15, 1.3e-8, 6.7e-11, 1.4e-10,
241 1, false);
242 }
243
244
245 @Test
246 void testRealFieldGradient() {
247
248
249
250 final Gradient a = Gradient.variable(6, 0, 7e6);
251 final Gradient e = Gradient.variable(6, 1, 0.01);
252 final Gradient i = Gradient.variable(6, 2, 1.2);
253 final Gradient pa = Gradient.variable(6, 3, 0.7);
254 final Gradient raan = Gradient.variable(6, 4, 0.5);
255 final Gradient meanAnomaly = Gradient.variable(6, 5, 0.1);
256
257 final Field<Gradient> field = a.getField();
258 final Gradient zero = field.getZero();
259
260
261 final FieldAbsoluteDate<Gradient> J2000 = new FieldAbsoluteDate<>(field);
262
263
264 final Frame EME = FramesFactory.getEME2000();
265
266
267 final FieldKeplerianOrbit<Gradient> FKO = new FieldKeplerianOrbit<>(a, e, i, pa, raan, meanAnomaly,
268 PositionAngleType.MEAN,
269 EME,
270 J2000,
271 zero.add(Constants.EIGEN5C_EARTH_MU));
272
273
274 final FieldSpacecraftState<Gradient> initialState = new FieldSpacecraftState<>(FKO);
275 final SpacecraftState iSR = initialState.toSpacecraftState();
276
277
278 final ClassicalRungeKuttaFieldIntegrator<Gradient> integrator =
279 new ClassicalRungeKuttaFieldIntegrator<>(field, zero.add(6));
280 final ClassicalRungeKuttaIntegrator RIntegrator =
281 new ClassicalRungeKuttaIntegrator(6);
282 final OrbitType type = OrbitType.KEPLERIAN;
283
284
285 final FieldNumericalPropagator<Gradient> FNP = new FieldNumericalPropagator<>(field, integrator);
286 FNP.setOrbitType(type);
287 FNP.setInitialState(initialState);
288
289 final NumericalPropagator NP = new NumericalPropagator(RIntegrator);
290 NP.setOrbitType(type);
291 NP.setInitialState(iSR);
292
293
294 final ExtendedPositionProvider sun = CelestialBodyFactory.getSun();
295
296
297 final RadiationSensitive radiationSensitive = new IsotropicRadiationSingleCoefficient(1, 1.5);
298
299
300 final KnockeRediffusedForceModel forceModel = new KnockeRediffusedForceModel(sun,
301 radiationSensitive,
302 Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
303 FastMath.toRadians(30));
304
305 FNP.addForceModel(forceModel);
306 NP.addForceModel(forceModel);
307
308
309 checkRealFieldPropagationGradient(FKO, PositionAngleType.MEAN, 300., NP, FNP,
310 1.0e-15, 1.3e-2, 9.6e-5, 1.4e-4,
311 1, false);
312 }
313
314
315
316
317
318
319 @Test
320 void testRoughAcceleration() {
321
322
323 final double mass = 406.9;
324 final double crossSection = 7E-4 * mass;
325 final double K = 1 + 0.13;
326
327 final TLE tle = new TLE ("1 08820U 76039 A 77047.52561960 .00000002 +00000-0 +00000-0 0 9994",
328 "2 08820 109.8332 127.3884 0044194 201.3006 158.6132 06.38663945018402");
329
330
331 final KeplerianOrbit keplerianTLE = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(TLEPropagator.selectExtrapolator(tle).
332 propagate(tle.getDate()).getOrbit());
333 final double a = keplerianTLE.getA();
334 final double e = keplerianTLE.getE();
335 final double i = keplerianTLE.getI();
336 final double pa = keplerianTLE.getPerigeeArgument();
337 final double raan = keplerianTLE.getRightAscensionOfAscendingNode();
338 final double nu = keplerianTLE.getTrueAnomaly();
339
340
341 final AbsoluteDate date0 = new AbsoluteDate(new DateComponents(1970, 1, 18),
342 new TimeComponents(0, 0, 0.0),
343 TimeScalesFactory.getUTC());
344
345
346 final Frame frame = FramesFactory.getTEME();
347
348 final KeplerianOrbit keplerian = new KeplerianOrbit(a, e, i, pa, raan, nu, PositionAngleType.TRUE, frame, date0, Constants.IERS2010_EARTH_MU);
349 final SpacecraftState initState = new SpacecraftState(keplerian).withMass(mass);
350
351
352
353
354 final ExtendedPositionProvider sun = CelestialBodyFactory.getSun();
355
356
357 final double equatorialRadius = Constants.EGM96_EARTH_EQUATORIAL_RADIUS;
358
359
360
361 final double angularResolution = FastMath.toRadians(15);
362 final RadiationSensitive radiationSensitive = new IsotropicRadiationSingleCoefficient(crossSection, K);
363 final KnockeRediffusedForceModel knockeModel = new KnockeRediffusedForceModel(sun, radiationSensitive,
364 equatorialRadius,
365 angularResolution);
366
367
368 final double duration = keplerian.getKeplerianPeriod();
369
370
371 final double minStep = duration * 0.01;
372 final double maxStep = duration * 0.5;
373 final double handlerStep = duration / 20;
374 final double positionTolerance = 1e-3;
375 final double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(positionTolerance).getTolerances(keplerian, OrbitType.KEPLERIAN);
376 AdaptiveStepsizeIntegrator integrator =
377 new DormandPrince853Integrator(minStep, maxStep, tolerances[0], tolerances[1]);
378
379 final NumericalPropagator propagator = new NumericalPropagator(integrator);
380
381 propagator.setInitialState(initState);
382 propagator.addForceModel(knockeModel);
383 propagator.setStepHandler(handlerStep, new KnockeStepHandler(knockeModel));
384
385 final SpacecraftState finalState = propagator.propagate(date0.shiftedBy(duration));
386
387 Assertions.assertEquals(date0.shiftedBy(duration), finalState.getDate());
388 }
389
390
391 private static class KnockeStepHandler implements OrekitFixedStepHandler {
392
393
394 private final KnockeRediffusedForceModel knockeModel;
395
396
397
398 KnockeStepHandler(final KnockeRediffusedForceModel knockeModel) {
399
400 this.knockeModel = knockeModel;
401 }
402
403 @Override
404 public void handleStep(SpacecraftState currentState) {
405
406
407
408 final Vector3D knockeAcceleration = knockeModel.acceleration(currentState, knockeModel.getParameters(currentState.getDate()));
409
410
411 final Vector3D radialUnit = currentState.getOrbit().getPosition().normalize();
412
413
414 final Vector3D velocity = currentState.getOrbit().getPVCoordinates().getVelocity();
415 final Vector3D alongTrackUnit = velocity.subtract(radialUnit.scalarMultiply(velocity.dotProduct(radialUnit))).normalize();
416
417
418 final Vector3D crossTrackUnit = radialUnit.crossProduct(alongTrackUnit);
419
420
421 final double radialAcceleration = knockeAcceleration.dotProduct(radialUnit);
422 final double alongTrackAcceleration = knockeAcceleration.dotProduct(alongTrackUnit);
423 final double crossTrackAcceleration = knockeAcceleration.dotProduct(crossTrackUnit);
424
425
426 Assertions.assertEquals(2.5e-10, radialAcceleration, 1.5e-10);
427 Assertions.assertEquals(0.0, alongTrackAcceleration, 5e-11);
428 Assertions.assertEquals(0.0, crossTrackAcceleration, 5e-12);
429 }
430
431 }
432 }