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