1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces.radiation;
18
19
20 import java.io.FileNotFoundException;
21 import java.text.ParseException;
22
23 import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
24 import org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator;
25 import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
26 import org.apache.commons.math3.util.FastMath;
27 import org.junit.Assert;
28 import org.junit.Before;
29 import org.junit.Test;
30 import org.orekit.Utils;
31 import org.orekit.bodies.CelestialBodyFactory;
32 import org.orekit.bodies.OneAxisEllipsoid;
33 import org.orekit.errors.OrekitException;
34 import org.orekit.errors.OrekitMessages;
35 import org.orekit.forces.AbstractForceModelTest;
36 import org.orekit.forces.BoxAndSolarArraySpacecraft;
37 import org.orekit.frames.FramesFactory;
38 import org.orekit.orbits.CartesianOrbit;
39 import org.orekit.orbits.EquinoctialOrbit;
40 import org.orekit.orbits.KeplerianOrbit;
41 import org.orekit.orbits.Orbit;
42 import org.orekit.orbits.OrbitType;
43 import org.orekit.orbits.PositionAngle;
44 import org.orekit.propagation.SpacecraftState;
45 import org.orekit.propagation.analytical.KeplerianPropagator;
46 import org.orekit.propagation.numerical.NumericalPropagator;
47 import org.orekit.propagation.sampling.OrekitFixedStepHandler;
48 import org.orekit.time.AbsoluteDate;
49 import org.orekit.time.DateComponents;
50 import org.orekit.time.TimeComponents;
51 import org.orekit.time.TimeScalesFactory;
52 import org.orekit.utils.Constants;
53 import org.orekit.utils.IERSConventions;
54 import org.orekit.utils.PVCoordinates;
55 import org.orekit.utils.PVCoordinatesProvider;
56
57
58 public class SolarRadiationPressureTest extends AbstractForceModelTest {
59
60 @Test
61 public void testLighting() throws OrekitException, ParseException {
62
63 AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 3, 21),
64 new TimeComponents(13, 59, 27.816),
65 TimeScalesFactory.getUTC());
66 Orbit orbit = new EquinoctialOrbit(42164000,10e-3,10e-3,
67 FastMath.tan(0.001745329)*FastMath.cos(2*FastMath.PI/3), FastMath.tan(0.001745329)*FastMath.sin(2*FastMath.PI/3),
68 0.1, PositionAngle.TRUE, FramesFactory.getEME2000(), date, mu);
69 PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
70 OneAxisEllipsoid earth =
71 new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765,
72 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
73 SolarRadiationPressure SRP =
74 new SolarRadiationPressure(sun, earth.getEquatorialRadius(),
75 (RadiationSensitive) new IsotropicRadiationCNES95Convention(50.0, 0.5, 0.5));
76
77 double period = 2*FastMath.PI*FastMath.sqrt(orbit.getA()*orbit.getA()*orbit.getA()/orbit.getMu());
78 Assert.assertEquals(86164, period,1);
79
80
81 KeplerianPropagator k = new KeplerianPropagator(orbit);
82
83
84 AbsoluteDate currentDate;
85 double changed = 1;
86 int count=0;
87
88 for(int t=1;t<3*period;t+=1000) {
89 currentDate = date.shiftedBy(t);
90 try {
91
92 double ratio = SRP.getLightingRatio(k.propagate(currentDate).getPVCoordinates().getPosition(),
93 FramesFactory.getEME2000(), currentDate);
94
95 if(FastMath.floor(ratio)!=changed) {
96 changed = FastMath.floor(ratio);
97 if(changed == 0) {
98 count++;
99 }
100 }
101 } catch (OrekitException e) {
102 e.printStackTrace();
103 }
104 }
105 Assert.assertTrue(3==count);
106 }
107
108 @Test
109 public void testParameterDerivativeIsotropicSingle() throws OrekitException {
110
111 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
112 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
113 final SpacecraftState state =
114 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
115 FramesFactory.getGCRF(),
116 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
117 Constants.EIGEN5C_EARTH_MU));
118
119 RadiationSensitive rs = new IsotropicRadiationSingleCoefficient(2.5, 0.7);
120 SolarRadiationPressure forceModel =
121 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
122 rs);
123
124 checkParameterDerivative(state, forceModel, RadiationSensitive.REFLECTION_COEFFICIENT, 1.0, 2.0e-15);
125
126 try {
127 rs.radiationPressureAcceleration(state.getDate(), state.getFrame(),
128 state.getPVCoordinates().getPosition(),
129 state.getAttitude().getRotation(),
130 state.getMass(), Vector3D.ZERO,
131 RadiationSensitive.ABSORPTION_COEFFICIENT);
132 Assert.fail("an exception should have been thrown");
133 } catch (OrekitException oe) {
134 Assert.assertEquals(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, oe.getSpecifier());
135 }
136 try {
137 Assert.assertEquals(0.0, rs.getAbsorptionCoefficient(), 1.0e-15);
138 rs.setAbsorptionCoefficient(0.3);
139 Assert.fail("an exception should have been thrown");
140 } catch (UnsupportedOperationException uso) {
141
142 }
143 }
144
145 @Test
146 public void testParameterDerivativeIsotropicClassical() throws OrekitException {
147
148 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
149 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
150 final SpacecraftState state =
151 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
152 FramesFactory.getGCRF(),
153 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
154 Constants.EIGEN5C_EARTH_MU));
155
156 RadiationSensitive rs = new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2);
157 SolarRadiationPressure forceModel =
158 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
159 rs);
160
161 checkParameterDerivative(state, forceModel, RadiationSensitive.ABSORPTION_COEFFICIENT, 1.0, 5.0e-16);
162 checkParameterDerivative(state, forceModel, RadiationSensitive.REFLECTION_COEFFICIENT, 1.0, 2.0e-15);
163
164 try {
165 rs.radiationPressureAcceleration(state.getDate(), state.getFrame(),
166 state.getPVCoordinates().getPosition(),
167 state.getAttitude().getRotation(),
168 state.getMass(), Vector3D.ZERO,
169 "UNKNOWN");
170 Assert.fail("an exception should have been thrown");
171 } catch (OrekitException oe) {
172 Assert.assertEquals(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, oe.getSpecifier());
173 }
174
175 }
176
177 @Test
178 public void testParameterDerivativeIsotropicCnes() throws OrekitException {
179
180 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
181 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
182 final SpacecraftState state =
183 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
184 FramesFactory.getGCRF(),
185 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
186 Constants.EIGEN5C_EARTH_MU));
187
188 RadiationSensitive rs = new IsotropicRadiationCNES95Convention(2.5, 0.7, 0.2);
189 SolarRadiationPressure forceModel =
190 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
191 rs);
192
193 checkParameterDerivative(state, forceModel, RadiationSensitive.ABSORPTION_COEFFICIENT, 1.0, 5.0e-16);
194 checkParameterDerivative(state, forceModel, RadiationSensitive.REFLECTION_COEFFICIENT, 1.0, 2.0e-15);
195
196 try {
197 rs.radiationPressureAcceleration(state.getDate(), state.getFrame(),
198 state.getPVCoordinates().getPosition(),
199 state.getAttitude().getRotation(),
200 state.getMass(), Vector3D.ZERO,
201 "UNKNOWN");
202 Assert.fail("an exception should have been thrown");
203 } catch (OrekitException oe) {
204 Assert.assertEquals(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, oe.getSpecifier());
205 }
206
207 }
208
209 @Test
210 public void testStateJacobianIsotropicSingle()
211 throws OrekitException {
212
213
214 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
215 new TimeComponents(13, 59, 27.816),
216 TimeScalesFactory.getUTC());
217 double i = FastMath.toRadians(98.7);
218 double omega = FastMath.toRadians(93.0);
219 double OMEGA = FastMath.toRadians(15.0 * 22.5);
220 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
221 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
222 Constants.EIGEN5C_EARTH_MU);
223 OrbitType integrationType = OrbitType.CARTESIAN;
224 double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
225
226 NumericalPropagator propagator =
227 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
228 tolerances[0], tolerances[1]));
229 propagator.setOrbitType(integrationType);
230 SolarRadiationPressure forceModel =
231 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
232 new IsotropicRadiationSingleCoefficient(2.5, 0.7));
233 propagator.addForceModel(forceModel);
234 SpacecraftState state0 = new SpacecraftState(orbit);
235
236 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
237 1e3, tolerances[0], 2.0e-6);
238
239 }
240
241 @Test
242 public void testStateJacobianIsotropicClassical()
243 throws OrekitException {
244
245
246 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
247 new TimeComponents(13, 59, 27.816),
248 TimeScalesFactory.getUTC());
249 double i = FastMath.toRadians(98.7);
250 double omega = FastMath.toRadians(93.0);
251 double OMEGA = FastMath.toRadians(15.0 * 22.5);
252 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
253 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
254 Constants.EIGEN5C_EARTH_MU);
255 OrbitType integrationType = OrbitType.CARTESIAN;
256 double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
257
258 NumericalPropagator propagator =
259 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
260 tolerances[0], tolerances[1]));
261 propagator.setOrbitType(integrationType);
262 SolarRadiationPressure forceModel =
263 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
264 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
265 propagator.addForceModel(forceModel);
266 SpacecraftState state0 = new SpacecraftState(orbit);
267
268 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
269 1e3, tolerances[0], 2.0e-6);
270
271 }
272
273 @Test
274 public void testStateJacobianIsotropicCnes()
275 throws OrekitException {
276
277
278 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
279 new TimeComponents(13, 59, 27.816),
280 TimeScalesFactory.getUTC());
281 double i = FastMath.toRadians(98.7);
282 double omega = FastMath.toRadians(93.0);
283 double OMEGA = FastMath.toRadians(15.0 * 22.5);
284 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
285 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
286 Constants.EIGEN5C_EARTH_MU);
287 OrbitType integrationType = OrbitType.CARTESIAN;
288 double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
289
290 NumericalPropagator propagator =
291 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
292 tolerances[0], tolerances[1]));
293 propagator.setOrbitType(integrationType);
294 SolarRadiationPressure forceModel =
295 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
296 new IsotropicRadiationCNES95Convention(2.5, 0.7, 0.2));
297 propagator.addForceModel(forceModel);
298 SpacecraftState state0 = new SpacecraftState(orbit);
299
300 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
301 1e3, tolerances[0], 2.0e-6);
302
303 }
304
305 @Test
306 public void testParameterDerivativeBox() throws OrekitException {
307
308 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
309 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
310 final SpacecraftState state =
311 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
312 FramesFactory.getGCRF(),
313 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
314 Constants.EIGEN5C_EARTH_MU));
315
316 SolarRadiationPressure forceModel =
317 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
318 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
319 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
320
321 checkParameterDerivative(state, forceModel, RadiationSensitive.ABSORPTION_COEFFICIENT, 1.0, 4.0e-16);
322 checkParameterDerivative(state, forceModel, RadiationSensitive.REFLECTION_COEFFICIENT, 1.0, 3.0e-15);
323
324 }
325
326 @Test
327 public void testStateJacobianBox()
328 throws OrekitException {
329
330
331 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
332 new TimeComponents(13, 59, 27.816),
333 TimeScalesFactory.getUTC());
334 double i = FastMath.toRadians(98.7);
335 double omega = FastMath.toRadians(93.0);
336 double OMEGA = FastMath.toRadians(15.0 * 22.5);
337 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
338 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
339 Constants.EIGEN5C_EARTH_MU);
340 OrbitType integrationType = OrbitType.CARTESIAN;
341 double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
342
343 NumericalPropagator propagator =
344 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
345 tolerances[0], tolerances[1]));
346 propagator.setOrbitType(integrationType);
347 SolarRadiationPressure forceModel =
348 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
349 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
350 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
351 propagator.addForceModel(forceModel);
352 SpacecraftState state0 = new SpacecraftState(orbit);
353
354 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
355 1e3, tolerances[0], 2.0e-5);
356
357 }
358
359 @Test
360 public void testRoughOrbitalModifs() throws ParseException, OrekitException, FileNotFoundException {
361
362
363 AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 7, 1),
364 new TimeComponents(13, 59, 27.816),
365 TimeScalesFactory.getUTC());
366 Orbit orbit = new EquinoctialOrbit(42164000,10e-3,10e-3,
367 FastMath.tan(0.001745329)*FastMath.cos(2*FastMath.PI/3),
368 FastMath.tan(0.001745329)*FastMath.sin(2*FastMath.PI/3),
369 0.1, PositionAngle.TRUE, FramesFactory.getEME2000(), date, mu);
370 final double period = orbit.getKeplerianPeriod();
371 Assert.assertEquals(86164, period, 1);
372 PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
373
374
375 OneAxisEllipsoid earth =
376 new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765,
377 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
378 SolarRadiationPressure SRP =
379 new SolarRadiationPressure(sun, earth.getEquatorialRadius(),
380 new IsotropicRadiationCNES95Convention(500.0, 0.7, 0.7));
381
382
383 double[] absTolerance = {
384 0.1, 1.0e-9, 1.0e-9, 1.0e-5, 1.0e-5, 1.0e-5, 0.001
385 };
386 double[] relTolerance = {
387 1.0e-4, 1.0e-4, 1.0e-4, 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-7
388 };
389 AdaptiveStepsizeIntegrator integrator =
390 new DormandPrince853Integrator(900.0, 60000, absTolerance, relTolerance);
391 integrator.setInitialStepSize(3600);
392 final NumericalPropagator calc = new NumericalPropagator(integrator);
393 calc.addForceModel(SRP);
394
395
396 calc.setMasterMode(FastMath.floor(period), new SolarStepHandler());
397 AbsoluteDate finalDate = date.shiftedBy(10 * period);
398 calc.setInitialState(new SpacecraftState(orbit, 1500.0));
399 calc.propagate(finalDate);
400 Assert.assertTrue(calc.getCalls() < 7100);
401 }
402
403 public static void checkRadius(double radius , double min , double max) {
404 Assert.assertTrue(radius >= min);
405 Assert.assertTrue(radius <= max);
406 }
407
408 private double mu = 3.98600E14;
409
410 private static class SolarStepHandler implements OrekitFixedStepHandler {
411
412 private SolarStepHandler() {
413 }
414
415 public void init(SpacecraftState s0, AbsoluteDate t) {
416 }
417
418 public void handleStep(SpacecraftState currentState, boolean isLast) {
419 final double dex = currentState.getEquinoctialEx() - 0.01071166;
420 final double dey = currentState.getEquinoctialEy() - 0.00654848;
421 final double alpha = FastMath.toDegrees(FastMath.atan2(dey, dex));
422 Assert.assertTrue(alpha > 100.0);
423 Assert.assertTrue(alpha < 112.0);
424 checkRadius(FastMath.sqrt(dex * dex + dey * dey), 0.003524, 0.003541);
425 }
426
427 }
428
429 @Before
430 public void setUp() {
431 Utils.setDataRoot("regular-data");
432 }
433
434 }