1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces;
18
19
20 import java.util.List;
21
22 import org.hamcrest.MatcherAssert;
23 import org.hipparchus.CalculusFieldElement;
24 import org.hipparchus.Field;
25 import org.hipparchus.geometry.euclidean.threed.FieldRotation;
26 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
27 import org.hipparchus.geometry.euclidean.threed.Rotation;
28 import org.hipparchus.geometry.euclidean.threed.Vector3D;
29 import org.hipparchus.util.Decimal64;
30 import org.hipparchus.util.Decimal64Field;
31 import org.hipparchus.util.FastMath;
32 import org.hipparchus.util.MathArrays;
33 import org.hipparchus.util.Precision;
34 import org.junit.Assert;
35 import org.junit.Before;
36 import org.junit.Test;
37 import org.orekit.OrekitMatchers;
38 import org.orekit.Utils;
39 import org.orekit.attitudes.LofOffset;
40 import org.orekit.bodies.CelestialBody;
41 import org.orekit.bodies.CelestialBodyFactory;
42 import org.orekit.errors.OrekitException;
43 import org.orekit.forces.drag.DragSensitive;
44 import org.orekit.forces.radiation.RadiationSensitive;
45 import org.orekit.frames.Frame;
46 import org.orekit.frames.FramesFactory;
47 import org.orekit.frames.LOFType;
48 import org.orekit.orbits.CircularOrbit;
49 import org.orekit.orbits.Orbit;
50 import org.orekit.orbits.PositionAngle;
51 import org.orekit.propagation.Propagator;
52 import org.orekit.propagation.SpacecraftState;
53 import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
54 import org.orekit.time.AbsoluteDate;
55 import org.orekit.time.DateComponents;
56 import org.orekit.time.FieldAbsoluteDate;
57 import org.orekit.time.TimeComponents;
58 import org.orekit.time.TimeScalesFactory;
59 import org.orekit.utils.Constants;
60 import org.orekit.utils.ParameterDriver;
61 import org.orekit.utils.TimeStampedPVCoordinates;
62
63 public class BoxAndSolarArraySpacecraftTest {
64
65 @Test
66 public void testParametersDrivers() {
67
68 CelestialBody sun = CelestialBodyFactory.getSun();
69 BoxAndSolarArraySpacecraft.Facet[] facets = new BoxAndSolarArraySpacecraft.Facet[] {
70 new BoxAndSolarArraySpacecraft.Facet(Vector3D.MINUS_I, 3.0),
71 new BoxAndSolarArraySpacecraft.Facet(Vector3D.PLUS_I, 3.0),
72 new BoxAndSolarArraySpacecraft.Facet(Vector3D.MINUS_J, 3.0),
73 new BoxAndSolarArraySpacecraft.Facet(Vector3D.PLUS_J, 3.0),
74 new BoxAndSolarArraySpacecraft.Facet(Vector3D.MINUS_K, 3.0),
75 new BoxAndSolarArraySpacecraft.Facet(Vector3D.PLUS_K, 3.0)
76 };
77
78 BoxAndSolarArraySpacecraft s1 =
79 new BoxAndSolarArraySpacecraft(1.5, 3.5, 2.5, sun, 20.0, Vector3D.PLUS_J, 2.0, 0.8, 0.1);
80 Assert.assertEquals(1, s1.getDragParametersDrivers().size());
81 Assert.assertEquals(DragSensitive.DRAG_COEFFICIENT, s1.getDragParametersDrivers().get(0).getName());
82 Assert.assertEquals(2.0, s1.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
83 Assert.assertEquals(2, s1.getRadiationParametersDrivers().size());
84 Assert.assertEquals(RadiationSensitive.ABSORPTION_COEFFICIENT, s1.getRadiationParametersDrivers().get(0).getName());
85 Assert.assertEquals(0.8, s1.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
86 Assert.assertEquals(RadiationSensitive.REFLECTION_COEFFICIENT, s1.getRadiationParametersDrivers().get(1).getName());
87 Assert.assertEquals(0.1, s1.getRadiationParametersDrivers().get(1).getValue(), 1.0e-15);
88
89 BoxAndSolarArraySpacecraft s2 =
90 new BoxAndSolarArraySpacecraft(1.5, 3.5, 2.5, sun, 20.0, Vector3D.PLUS_J, 2.0, 0.4, 0.8, 0.1);
91 Assert.assertEquals(2, s2.getDragParametersDrivers().size());
92 Assert.assertEquals(DragSensitive.DRAG_COEFFICIENT, s2.getDragParametersDrivers().get(0).getName());
93 Assert.assertEquals(2.0, s2.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
94 Assert.assertEquals(DragSensitive.LIFT_RATIO, s2.getDragParametersDrivers().get(1).getName());
95 Assert.assertEquals(0.4, s2.getDragParametersDrivers().get(1).getValue(), 1.0e-15);
96 Assert.assertEquals(2, s2.getRadiationParametersDrivers().size());
97 Assert.assertEquals(RadiationSensitive.ABSORPTION_COEFFICIENT, s2.getRadiationParametersDrivers().get(0).getName());
98 Assert.assertEquals(0.8, s2.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
99 Assert.assertEquals(RadiationSensitive.REFLECTION_COEFFICIENT, s2.getRadiationParametersDrivers().get(1).getName());
100 Assert.assertEquals(0.1, s2.getRadiationParametersDrivers().get(1).getValue(), 1.0e-15);
101
102 BoxAndSolarArraySpacecraft s3 =
103 new BoxAndSolarArraySpacecraft(facets, sun, 20.0, Vector3D.PLUS_J, 2.0, 0.8, 0.1);
104 Assert.assertEquals(1, s3.getDragParametersDrivers().size());
105 Assert.assertEquals(DragSensitive.DRAG_COEFFICIENT, s3.getDragParametersDrivers().get(0).getName());
106 Assert.assertEquals(2.0, s3.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
107 Assert.assertEquals(2, s3.getRadiationParametersDrivers().size());
108 Assert.assertEquals(RadiationSensitive.ABSORPTION_COEFFICIENT, s3.getRadiationParametersDrivers().get(0).getName());
109 Assert.assertEquals(0.8, s3.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
110 Assert.assertEquals(RadiationSensitive.REFLECTION_COEFFICIENT, s3.getRadiationParametersDrivers().get(1).getName());
111 Assert.assertEquals(0.1, s3.getRadiationParametersDrivers().get(1).getValue(), 1.0e-15);
112
113 BoxAndSolarArraySpacecraft s4 =
114 new BoxAndSolarArraySpacecraft(facets, sun, 20.0, Vector3D.PLUS_J, 2.0, 0.4, 0.8, 0.1);
115 Assert.assertEquals(2, s4.getDragParametersDrivers().size());
116 Assert.assertEquals(DragSensitive.DRAG_COEFFICIENT, s4.getDragParametersDrivers().get(0).getName());
117 Assert.assertEquals(2.0, s4.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
118 Assert.assertEquals(DragSensitive.LIFT_RATIO, s4.getDragParametersDrivers().get(1).getName());
119 Assert.assertEquals(0.4, s4.getDragParametersDrivers().get(1).getValue(), 1.0e-15);
120 Assert.assertEquals(2, s4.getRadiationParametersDrivers().size());
121 Assert.assertEquals(RadiationSensitive.ABSORPTION_COEFFICIENT, s4.getRadiationParametersDrivers().get(0).getName());
122 Assert.assertEquals(0.8, s4.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
123 Assert.assertEquals(RadiationSensitive.REFLECTION_COEFFICIENT, s4.getRadiationParametersDrivers().get(1).getName());
124 Assert.assertEquals(0.1, s4.getRadiationParametersDrivers().get(1).getValue(), 1.0e-15);
125
126 BoxAndSolarArraySpacecraft s5 =
127 new BoxAndSolarArraySpacecraft(1.5, 3.5, 2.5, sun, 20.0, Vector3D.PLUS_J,
128 AbsoluteDate.J2000_EPOCH, Vector3D.PLUS_I, 7.292e-5,
129 2.0, 0.8, 0.1);
130 Assert.assertEquals(1, s5.getDragParametersDrivers().size());
131 Assert.assertEquals(DragSensitive.DRAG_COEFFICIENT, s5.getDragParametersDrivers().get(0).getName());
132 Assert.assertEquals(2.0, s5.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
133 Assert.assertEquals(2, s5.getRadiationParametersDrivers().size());
134 Assert.assertEquals(RadiationSensitive.ABSORPTION_COEFFICIENT, s5.getRadiationParametersDrivers().get(0).getName());
135 Assert.assertEquals(0.8, s5.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
136 Assert.assertEquals(RadiationSensitive.REFLECTION_COEFFICIENT, s5.getRadiationParametersDrivers().get(1).getName());
137 Assert.assertEquals(0.1, s5.getRadiationParametersDrivers().get(1).getValue(), 1.0e-15);
138
139 BoxAndSolarArraySpacecraft s6 =
140 new BoxAndSolarArraySpacecraft(1.5, 3.5, 2.5, sun, 20.0, Vector3D.PLUS_J,
141 AbsoluteDate.J2000_EPOCH, Vector3D.PLUS_I, 7.292e-5,
142 2.0, 0.4, 0.8, 0.1);
143 Assert.assertEquals(2, s6.getDragParametersDrivers().size());
144 Assert.assertEquals(DragSensitive.DRAG_COEFFICIENT, s6.getDragParametersDrivers().get(0).getName());
145 Assert.assertEquals(2.0, s6.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
146 Assert.assertEquals(DragSensitive.LIFT_RATIO, s6.getDragParametersDrivers().get(1).getName());
147 Assert.assertEquals(0.4, s6.getDragParametersDrivers().get(1).getValue(), 1.0e-15);
148 Assert.assertEquals(2, s6.getRadiationParametersDrivers().size());
149 Assert.assertEquals(RadiationSensitive.ABSORPTION_COEFFICIENT, s6.getRadiationParametersDrivers().get(0).getName());
150 Assert.assertEquals(0.8, s6.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
151 Assert.assertEquals(RadiationSensitive.REFLECTION_COEFFICIENT, s6.getRadiationParametersDrivers().get(1).getName());
152 Assert.assertEquals(0.1, s6.getRadiationParametersDrivers().get(1).getValue(), 1.0e-15);
153
154 BoxAndSolarArraySpacecraft s7 =
155 new BoxAndSolarArraySpacecraft(facets, sun, 20.0, Vector3D.PLUS_J,
156 AbsoluteDate.J2000_EPOCH, Vector3D.PLUS_I, 7.292e-5,
157 2.0, 0.8, 0.1);
158 Assert.assertEquals(1, s7.getDragParametersDrivers().size());
159 Assert.assertEquals(DragSensitive.DRAG_COEFFICIENT, s7.getDragParametersDrivers().get(0).getName());
160 Assert.assertEquals(2.0, s7.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
161 Assert.assertEquals(2, s7.getRadiationParametersDrivers().size());
162 Assert.assertEquals(RadiationSensitive.ABSORPTION_COEFFICIENT, s7.getRadiationParametersDrivers().get(0).getName());
163 Assert.assertEquals(0.8, s7.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
164 Assert.assertEquals(RadiationSensitive.REFLECTION_COEFFICIENT, s7.getRadiationParametersDrivers().get(1).getName());
165 Assert.assertEquals(0.1, s7.getRadiationParametersDrivers().get(1).getValue(), 1.0e-15);
166
167 BoxAndSolarArraySpacecraft s8 =
168 new BoxAndSolarArraySpacecraft(facets, sun, 20.0, Vector3D.PLUS_J,
169 AbsoluteDate.J2000_EPOCH, Vector3D.PLUS_I, 7.292e-5,
170 2.0, 0.4, 0.8, 0.1);
171 Assert.assertEquals(2, s8.getDragParametersDrivers().size());
172 Assert.assertEquals(DragSensitive.DRAG_COEFFICIENT, s8.getDragParametersDrivers().get(0).getName());
173 Assert.assertEquals(2.0, s8.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
174 Assert.assertEquals(DragSensitive.LIFT_RATIO, s8.getDragParametersDrivers().get(1).getName());
175 Assert.assertEquals(0.4, s8.getDragParametersDrivers().get(1).getValue(), 1.0e-15);
176 Assert.assertEquals(2, s8.getRadiationParametersDrivers().size());
177 Assert.assertEquals(RadiationSensitive.ABSORPTION_COEFFICIENT, s8.getRadiationParametersDrivers().get(0).getName());
178 Assert.assertEquals(0.8, s8.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
179 Assert.assertEquals(RadiationSensitive.REFLECTION_COEFFICIENT, s8.getRadiationParametersDrivers().get(1).getName());
180 Assert.assertEquals(0.1, s8.getRadiationParametersDrivers().get(1).getValue(), 1.0e-15);
181
182 }
183
184 @Test
185 public void testBestPointing() {
186
187 AbsoluteDate initialDate = propagator.getInitialState().getDate();
188 CelestialBody sun = CelestialBodyFactory.getSun();
189 BoxAndSolarArraySpacecraft s =
190 new BoxAndSolarArraySpacecraft(1.5, 3.5, 2.5, sun, 20.0, Vector3D.PLUS_J, 0.0, 0.0, 0.0);
191 for (double dt = 0; dt < 4000; dt += 60) {
192
193 SpacecraftState state = propagator.propagate(initialDate.shiftedBy(dt));
194
195 Vector3D sunInert = sun.getPVCoordinates(initialDate, state.getFrame()).getPosition();
196 Vector3D momentum = state.getPVCoordinates().getMomentum();
197 double sunElevation = FastMath.PI / 2 - Vector3D.angle(sunInert, momentum);
198 Assert.assertEquals(15.1, FastMath.toDegrees(sunElevation), 0.1);
199
200 Vector3D n = s.getNormal(state.getDate(), state.getFrame(),
201 state.getPVCoordinates().getPosition(),
202 state.getAttitude().getRotation());
203 Assert.assertEquals(0.0, n.getY(), 1.0e-10);
204
205
206 Vector3D sunSat = state.getAttitude().getRotation().applyTo(sunInert);
207 double misAlignment = Vector3D.angle(sunSat, n);
208 Assert.assertEquals(sunElevation, misAlignment, 1.0e-3);
209
210 }
211 }
212
213 @Test
214 public void testCorrectFixedRate() {
215
216 AbsoluteDate initialDate = propagator.getInitialState().getDate();
217 CelestialBody sun = CelestialBodyFactory.getSun();
218 BoxAndSolarArraySpacecraft s =
219 new BoxAndSolarArraySpacecraft(1.5, 3.5, 2.5, sun, 20.0, Vector3D.PLUS_J,
220 initialDate,
221 new Vector3D(0.46565509814462996, 0.0, 0.884966287251619),
222 propagator.getInitialState().getKeplerianMeanMotion(),
223 0.0, 0.0, 0.0);
224
225 for (double dt = 0; dt < 4000; dt += 60) {
226
227 SpacecraftState state = propagator.propagate(initialDate.shiftedBy(dt));
228
229 Vector3D sunInert = sun.getPVCoordinates(initialDate, state.getFrame()).getPosition();
230 Vector3D momentum = state.getPVCoordinates().getMomentum();
231 double sunElevation = FastMath.PI / 2 - Vector3D.angle(sunInert, momentum);
232 Assert.assertEquals(15.1, FastMath.toDegrees(sunElevation), 0.1);
233
234 Vector3D n = s.getNormal(state.getDate(), state.getFrame(),
235 state.getPVCoordinates().getPosition(),
236 state.getAttitude().getRotation());
237 Assert.assertEquals(0.0, n.getY(), 1.0e-10);
238
239
240 Vector3D sunSat = state.getAttitude().getRotation().applyTo(sunInert);
241 double misAlignment = Vector3D.angle(sunSat, n);
242 Assert.assertEquals(sunElevation, misAlignment, 1.0e-3);
243
244 }
245 }
246
247 @Test
248 public void testTooSlowFixedRate() {
249
250 AbsoluteDate initialDate = propagator.getInitialState().getDate();
251 CelestialBody sun = CelestialBodyFactory.getSun();
252 BoxAndSolarArraySpacecraft s =
253 new BoxAndSolarArraySpacecraft(1.5, 3.5, 2.5, sun, 20.0, Vector3D.PLUS_J,
254 initialDate,
255 new Vector3D(0.46565509814462996, 0.0, 0.884966287251619),
256 0.1 * propagator.getInitialState().getKeplerianMeanMotion(),
257 0.0, 0.0, 0.0);
258
259 double maxDelta = 0;
260 for (double dt = 0; dt < 4000; dt += 60) {
261
262 SpacecraftState state = propagator.propagate(initialDate.shiftedBy(dt));
263
264 Vector3D sunInert = sun.getPVCoordinates(initialDate, state.getFrame()).getPosition();
265 Vector3D momentum = state.getPVCoordinates().getMomentum();
266 double sunElevation = FastMath.PI / 2 - Vector3D.angle(sunInert, momentum);
267 Assert.assertEquals(15.1, FastMath.toDegrees(sunElevation), 0.1);
268
269 Vector3D n = s.getNormal(state.getDate(), state.getFrame(),
270 state.getPVCoordinates().getPosition(),
271 state.getAttitude().getRotation());
272 Assert.assertEquals(0.0, n.getY(), 1.0e-10);
273
274
275 Vector3D sunSat = state.getAttitude().getRotation().applyTo(sunInert);
276 double misAlignment = Vector3D.angle(sunSat, n);
277 maxDelta = FastMath.max(maxDelta, FastMath.abs(sunElevation - misAlignment));
278
279 }
280 Assert.assertTrue(FastMath.toDegrees(maxDelta) > 120.0);
281
282 }
283
284 @Test
285 public void testNoLiftWithoutReflection() {
286
287 AbsoluteDate initialDate = propagator.getInitialState().getDate();
288 CelestialBody sun = CelestialBodyFactory.getSun();
289 BoxAndSolarArraySpacecraft s =
290 new BoxAndSolarArraySpacecraft(1.5, 3.5, 2.5, sun, 20.0, Vector3D.PLUS_J, 1.0, 0.0, 1.0, 0.0);
291
292 Vector3D earthRot = new Vector3D(0.0, 0.0, 7.292115e-4);
293 for (double dt = 0; dt < 4000; dt += 60) {
294
295 AbsoluteDate date = initialDate.shiftedBy(dt);
296 SpacecraftState state = propagator.propagate(date);
297
298
299 Vector3D p = state.getPVCoordinates().getPosition();
300 Vector3D v = state.getPVCoordinates().getVelocity();
301 Vector3D vAtm = Vector3D.crossProduct(earthRot, p);
302 Vector3D relativeVelocity = vAtm.subtract(v);
303
304 Vector3D drag = s.dragAcceleration(state.getDate(), state.getFrame(),
305 state.getPVCoordinates().getPosition(),
306 state.getAttitude().getRotation(),
307 state.getMass(), 0.001, relativeVelocity,
308 getDragParameters(s));
309 Assert.assertEquals(0.0, Vector3D.angle(relativeVelocity, drag), 1.0e-15);
310
311 Vector3D sunDirection = sun.getPVCoordinates(date, state.getFrame()).getPosition().normalize();
312 Vector3D flux = new Vector3D(-4.56e-6, sunDirection);
313 Vector3D radiation = s.radiationPressureAcceleration(state.getDate(), state.getFrame(),
314 state.getPVCoordinates().getPosition(),
315 state.getAttitude().getRotation(),
316 state.getMass(), flux,
317 getRadiationParameters(s));
318 Assert.assertEquals(0.0, Vector3D.angle(flux, radiation), 1.0e-9);
319
320 }
321
322 }
323
324 @Test
325 public void testOnlyLiftWithoutReflection() {
326
327 AbsoluteDate initialDate = propagator.getInitialState().getDate();
328 CelestialBody sun = CelestialBodyFactory.getSun();
329 BoxAndSolarArraySpacecraft s =
330 new BoxAndSolarArraySpacecraft(1.5, 3.5, 2.5, sun, 20.0, Vector3D.PLUS_J, 1.0, 1.0, 1.0, 0.0);
331
332 Vector3D earthRot = new Vector3D(0.0, 0.0, 7.292115e-4);
333 for (double dt = 0; dt < 4000; dt += 60) {
334
335 AbsoluteDate date = initialDate.shiftedBy(dt);
336 SpacecraftState state = propagator.propagate(date);
337
338
339 Vector3D p = state.getPVCoordinates().getPosition();
340 Vector3D v = state.getPVCoordinates().getVelocity();
341 Vector3D vAtm = Vector3D.crossProduct(earthRot, p);
342 Vector3D relativeVelocity = vAtm.subtract(v);
343
344 Vector3D drag = s.dragAcceleration(state.getDate(), state.getFrame(),
345 state.getPVCoordinates().getPosition(),
346 state.getAttitude().getRotation(),
347 state.getMass(), 0.001, relativeVelocity,
348 getDragParameters(s));
349 Assert.assertTrue(Vector3D.angle(relativeVelocity, drag) > 0.167);
350 Assert.assertTrue(Vector3D.angle(relativeVelocity, drag) < 0.736);
351
352 Vector3D sunDirection = sun.getPVCoordinates(date, state.getFrame()).getPosition().normalize();
353 Vector3D flux = new Vector3D(-4.56e-6, sunDirection);
354 Vector3D radiation = s.radiationPressureAcceleration(state.getDate(), state.getFrame(),
355 state.getPVCoordinates().getPosition(),
356 state.getAttitude().getRotation(),
357 state.getMass(), flux,
358 getRadiationParameters(s));
359 Assert.assertEquals(0.0, Vector3D.angle(flux, radiation), 1.0e-9);
360
361 }
362
363 }
364
365 @Test
366 public void testLiftVsNoLift()
367 throws NoSuchFieldException, SecurityException,
368 IllegalArgumentException, IllegalAccessException {
369
370 CelestialBody sun = CelestialBodyFactory.getSun();
371
372
373
374
375
376 BoxAndSolarArraySpacecraft.Facet[] facets = new BoxAndSolarArraySpacecraft.Facet[] {
377 new BoxAndSolarArraySpacecraft.Facet(Vector3D.MINUS_I, 3.0),
378 new BoxAndSolarArraySpacecraft.Facet(Vector3D.PLUS_I, 3.0),
379 new BoxAndSolarArraySpacecraft.Facet(Vector3D.MINUS_J, 3.0),
380 new BoxAndSolarArraySpacecraft.Facet(Vector3D.PLUS_J, 3.0),
381 new BoxAndSolarArraySpacecraft.Facet(Vector3D.MINUS_K, 3.0),
382 new BoxAndSolarArraySpacecraft.Facet(Vector3D.PLUS_K, 3.0)
383 };
384 BoxAndSolarArraySpacecraft cube =
385 new BoxAndSolarArraySpacecraft(facets, sun, 0.0, Vector3D.PLUS_J, 1.0, 1.0, 1.0, 0.0);
386
387 AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
388 Frame frame = FramesFactory.getEME2000();
389 Vector3D position = new Vector3D(1234567.8, 9876543.21, 121212.3434);
390 double mass = 1000.0;
391 double density = 0.001;
392 Rotation rotation = Rotation.IDENTITY;
393
394
395 Vector3D headOnVelocity = new Vector3D(2000, 0.0, 0.0);
396 Vector3D newHeadOnDrag = cube.dragAcceleration(date, frame, position, rotation, mass, density, headOnVelocity,
397 getDragParameters(cube));
398 Vector3D oldHeadOnDrag = oldDragAcceleration(cube, date, frame, position, rotation, mass, density, headOnVelocity);
399 MatcherAssert.assertThat(newHeadOnDrag, OrekitMatchers.vectorCloseTo(oldHeadOnDrag.scalarMultiply(2), 1));
400
401
402
403
404
405
406
407 Vector3D diagonalVelocity = new Vector3D(2000, 2000, 2000);
408 Vector3D newDiagDrag= cube.dragAcceleration(date, frame, position, rotation, mass, density, diagonalVelocity,
409 getDragParameters(cube));
410 Vector3D oldDiagDrag = oldDragAcceleration(cube, date, frame, position, rotation, mass, density, diagonalVelocity);
411 double oldMissingCoeff = 2.0 / FastMath.sqrt(3.0);
412 Vector3D fixedOldDrag = new Vector3D(oldMissingCoeff, oldDiagDrag);
413 MatcherAssert.assertThat(newDiagDrag, OrekitMatchers.vectorCloseTo(fixedOldDrag, 1));
414
415 }
416
417
418
419
420 private Vector3D oldDragAcceleration(final BoxAndSolarArraySpacecraft bsa,
421 final AbsoluteDate date, final Frame frame, final Vector3D position,
422 final Rotation rotation, final double mass,
423 final double density, final Vector3D relativeVelocity)
424 throws IllegalArgumentException, IllegalAccessException,
425 NoSuchFieldException, SecurityException {
426
427 java.lang.reflect.Field facetsField = BoxAndSolarArraySpacecraft.class.getDeclaredField("facets");
428 facetsField.setAccessible(true);
429 @SuppressWarnings("unchecked")
430 final List<BoxAndSolarArraySpacecraft.Facet> facets = (List<BoxAndSolarArraySpacecraft.Facet>) facetsField.get(bsa);
431
432 java.lang.reflect.Field saAreaField = BoxAndSolarArraySpacecraft.class.getDeclaredField("solarArrayArea");
433 saAreaField.setAccessible(true);
434 final double solarArrayArea = (Double) saAreaField.get(bsa);
435
436 final double dragCoeff = bsa.getDragParametersDrivers().get(0).getValue();
437
438
439 final Vector3D v = rotation.applyTo(relativeVelocity);
440
441
442 final Vector3D solarArrayFacet = new Vector3D(solarArrayArea, bsa.getNormal(date, frame, position, rotation));
443 double sv = FastMath.abs(Vector3D.dotProduct(solarArrayFacet, v));
444
445
446 for (final BoxAndSolarArraySpacecraft.Facet facet : facets) {
447 final double dot = Vector3D.dotProduct(facet.getNormal(), v);
448 if (dot < 0) {
449
450 sv -= facet.getArea() * dot;
451 }
452 }
453
454 return new Vector3D(sv * density * dragCoeff / (2.0 * mass), relativeVelocity);
455
456 }
457
458 @Test
459 public void testPlaneSpecularReflection() {
460
461 AbsoluteDate initialDate = propagator.getInitialState().getDate();
462 CelestialBody sun = CelestialBodyFactory.getSun();
463 BoxAndSolarArraySpacecraft s =
464 new BoxAndSolarArraySpacecraft(0, 0, 0, sun, 20.0, Vector3D.PLUS_J, 0.0, 0.0, 1.0);
465
466 for (double dt = 0; dt < 4000; dt += 60) {
467
468 AbsoluteDate date = initialDate.shiftedBy(dt);
469 SpacecraftState state = propagator.propagate(date);
470
471 Vector3D sunDirection = sun.getPVCoordinates(date, state.getFrame()).getPosition().normalize();
472 Vector3D flux = new Vector3D(-4.56e-6, sunDirection);
473 Vector3D acceleration = s.radiationPressureAcceleration(state.getDate(), state.getFrame(),
474 state.getPVCoordinates().getPosition(),
475 state.getAttitude().getRotation(),
476 state.getMass(), flux,
477 getRadiationParameters(s));
478 Vector3D normal = state.getAttitude().getRotation().applyInverseTo(s.getNormal(state.getDate(), state.getFrame(),
479 state.getPVCoordinates().getPosition(),
480 state.getAttitude().getRotation()));
481
482
483 Assert.assertEquals(15.1, FastMath.toDegrees(Vector3D.angle(sunDirection, normal)), 0.11);
484
485
486 Assert.assertEquals(180.0, FastMath.toDegrees(Vector3D.angle(acceleration, normal)), 1.0e-3);
487
488 }
489
490 }
491
492 @Test
493 public void testPlaneAbsorption() {
494
495 AbsoluteDate initialDate = propagator.getInitialState().getDate();
496 CelestialBody sun = CelestialBodyFactory.getSun();
497 BoxAndSolarArraySpacecraft s =
498 new BoxAndSolarArraySpacecraft(0, 0, 0, sun, 20.0, Vector3D.PLUS_J, 0.0, 1.0, 0.0);
499
500 for (double dt = 0; dt < 4000; dt += 60) {
501
502 AbsoluteDate date = initialDate.shiftedBy(dt);
503 SpacecraftState state = propagator.propagate(date);
504
505 Vector3D sunDirection = sun.getPVCoordinates(date, state.getFrame()).getPosition().normalize();
506 Vector3D flux = new Vector3D(-4.56e-6, sunDirection);
507 Vector3D acceleration =
508 s.radiationPressureAcceleration(state.getDate(), state.getFrame(),
509 state.getPVCoordinates().getPosition(),
510 state.getAttitude().getRotation(),
511 state.getMass(), flux,
512 getRadiationParameters(s));
513 Vector3D normal = state.getAttitude().getRotation().applyInverseTo(s.getNormal(state.getDate(), state.getFrame(),
514 state.getPVCoordinates().getPosition(),
515 state.getAttitude().getRotation()));
516
517
518 Assert.assertEquals(15.1, FastMath.toDegrees(Vector3D.angle(sunDirection, normal)), 0.11);
519
520
521 Assert.assertEquals(180.0, FastMath.toDegrees(Vector3D.angle(acceleration, sunDirection)), 1.0e-3);
522
523 }
524
525 }
526
527
528 @Test
529 public void testNullIllumination() {
530 SpacecraftState state = propagator.getInitialState();
531 CelestialBody sun = CelestialBodyFactory.getSun();
532 BoxAndSolarArraySpacecraft s =
533 new BoxAndSolarArraySpacecraft(0, 0, 0, sun, 20.0, Vector3D.PLUS_J, 0.0, 1.0, 0.0);
534
535
536 Field<Decimal64> field = Decimal64Field.getInstance();
537 Decimal64[] srpParam = getRadiationParameters(s, field);
538
539 FieldAbsoluteDate<Decimal64> date = new FieldAbsoluteDate<>(field, state.getDate());
540 FieldVector3D<Decimal64> position = new FieldVector3D<Decimal64>(field.getOne(), state.getPVCoordinates().getPosition());
541 FieldRotation<Decimal64> rotation = new FieldRotation<>(field, state.getAttitude().getRotation());
542 Decimal64 mass = new Decimal64(state.getMass());
543 FieldVector3D<Decimal64> flux = new FieldVector3D<Decimal64>(field.getOne(),
544 new Vector3D(Precision.SAFE_MIN / 2, Vector3D.PLUS_I));
545
546
547 FieldVector3D<Decimal64> a = s.radiationPressureAcceleration(date, state.getFrame(),
548 position, rotation, mass,
549 flux, srpParam);
550 Assert.assertEquals(0.0, a.getNorm().getReal(), Double.MIN_VALUE);
551 }
552
553
554 @Test
555 public void testBackwardIllumination() {
556 SpacecraftState state = propagator.getInitialState();
557 CelestialBody sun = CelestialBodyFactory.getSun();
558 BoxAndSolarArraySpacecraft s =
559 new BoxAndSolarArraySpacecraft(0, 0, 0, sun, 20.0, Vector3D.PLUS_J, 0.0, 1.0, 0.0);
560
561
562 Field<Decimal64> field = Decimal64Field.getInstance();
563 Decimal64[] srpParam = getRadiationParameters(s, field);
564
565 FieldAbsoluteDate<Decimal64> date = new FieldAbsoluteDate<>(field, state.getDate());
566 FieldVector3D<Decimal64> position = new FieldVector3D<Decimal64>(field.getOne(), state.getPVCoordinates().getPosition());
567 FieldRotation<Decimal64> rotation = new FieldRotation<>(field, state.getAttitude().getRotation());
568 Decimal64 mass = new Decimal64(state.getMass());
569
570
571 FieldVector3D<Decimal64> flux = s.getNormal(date, state.getFrame(), position, rotation);
572
573
574 FieldVector3D<Decimal64> aPlus = s.radiationPressureAcceleration(date, state.getFrame(),
575 position, rotation, mass,
576 flux, srpParam);
577
578 FieldVector3D<Decimal64> aMinus = s.radiationPressureAcceleration(date, state.getFrame(),
579 position, rotation, mass,
580 flux.negate(), srpParam);
581
582 Assert.assertEquals(0.0, aPlus.add(aMinus).getNorm().getReal(), Double.MIN_VALUE);
583 }
584
585 @Test
586 public void testNormalOptimalRotationDouble() {
587 AbsoluteDate initialDate = propagator.getInitialState().getDate();
588 CelestialBody sun = CelestialBodyFactory.getSun();
589 BoxAndSolarArraySpacecraft s =
590 new BoxAndSolarArraySpacecraft(0, 0, 0, sun, 20.0, Vector3D.PLUS_J, 0.0, 1.0, 0.0);
591 for (double dt = 0; dt < 4000; dt += 60) {
592 AbsoluteDate date = initialDate.shiftedBy(dt);
593 SpacecraftState state = propagator.propagate(date);
594 Vector3D normal = s.getNormal(state.getDate(), state.getFrame(),
595 state.getPVCoordinates().getPosition(),
596 state.getAttitude().getRotation());
597 Assert.assertEquals(0, Vector3D.dotProduct(normal, Vector3D.PLUS_J), 1.0e-16);
598 }
599 }
600
601 @Test
602 public void testNormalOptimalRotationField() {
603 AbsoluteDate initialDate = propagator.getInitialState().getDate();
604 CelestialBody sun = CelestialBodyFactory.getSun();
605 BoxAndSolarArraySpacecraft s =
606 new BoxAndSolarArraySpacecraft(0, 0, 0, sun, 20.0, Vector3D.PLUS_J, 0.0, 1.0, 0.0);
607 Field<Decimal64> field = Decimal64Field.getInstance();
608 for (double dt = 0; dt < 4000; dt += 60) {
609 AbsoluteDate date = initialDate.shiftedBy(dt);
610 SpacecraftState state = propagator.propagate(date);
611 FieldVector3D<Decimal64> normal = s.getNormal(new FieldAbsoluteDate<>(field, state.getDate()),
612 state.getFrame(),
613 new FieldVector3D<>(field, state.getPVCoordinates().getPosition()),
614 new FieldRotation<>(field, state.getAttitude().getRotation()));
615 Assert.assertEquals(0, FieldVector3D.dotProduct(normal, Vector3D.PLUS_J).getReal(), 1.0e-16);
616 }
617 }
618
619 @Test
620 public void testNormalFixedRateDouble() {
621 AbsoluteDate initialDate = propagator.getInitialState().getDate();
622 CelestialBody sun = CelestialBodyFactory.getSun();
623 BoxAndSolarArraySpacecraft s =
624 new BoxAndSolarArraySpacecraft(0, 0, 0, sun, 20.0, Vector3D.PLUS_J,
625 initialDate, Vector3D.PLUS_K, 1.0e-3,
626 0.0, 1.0, 0.0);
627 for (double dt = 0; dt < 4000; dt += 60) {
628 AbsoluteDate date = initialDate.shiftedBy(dt);
629 SpacecraftState state = propagator.propagate(date);
630 Vector3D normal = s.getNormal(state.getDate(), state.getFrame(),
631 state.getPVCoordinates().getPosition(),
632 state.getAttitude().getRotation());
633 Assert.assertEquals(0, Vector3D.dotProduct(normal, Vector3D.PLUS_J), 1.0e-16);
634 }
635 }
636
637 @Test
638 public void testNormalFixedRateField() {
639 AbsoluteDate initialDate = propagator.getInitialState().getDate();
640 CelestialBody sun = CelestialBodyFactory.getSun();
641 BoxAndSolarArraySpacecraft s =
642 new BoxAndSolarArraySpacecraft(0, 0, 0, sun, 20.0, Vector3D.PLUS_J,
643 initialDate, Vector3D.PLUS_K, 1.0e-3,
644 0.0, 1.0, 0.0);
645 Field<Decimal64> field = Decimal64Field.getInstance();
646 for (double dt = 0; dt < 4000; dt += 60) {
647 AbsoluteDate date = initialDate.shiftedBy(dt);
648 SpacecraftState state = propagator.propagate(date);
649 FieldVector3D<Decimal64> normal = s.getNormal(new FieldAbsoluteDate<>(field, state.getDate()),
650 state.getFrame(),
651 new FieldVector3D<>(field, state.getPVCoordinates().getPosition()),
652 new FieldRotation<>(field, state.getAttitude().getRotation()));
653 Assert.assertEquals(0, FieldVector3D.dotProduct(normal, Vector3D.PLUS_J).getReal(), 1.0e-16);
654 }
655 }
656
657 @Test
658 public void testNormalSunAlignedDouble() {
659 BoxAndSolarArraySpacecraft s =
660 new BoxAndSolarArraySpacecraft(0, 0, 0,
661 (date, frame) -> new TimeStampedPVCoordinates(date, new Vector3D(0, 1e6, 0), Vector3D.ZERO),
662 20.0, Vector3D.PLUS_J, 0.0, 1.0, 0.0);
663 Vector3D normal = s.getNormal(AbsoluteDate.J2000_EPOCH, FramesFactory.getEME2000(),
664 Vector3D.ZERO, Rotation.IDENTITY);
665 Assert.assertEquals(0, Vector3D.dotProduct(normal, Vector3D.PLUS_J), 1.0e-16);
666 }
667
668 @Test
669 public void testNormalSunAlignedField() {
670 BoxAndSolarArraySpacecraft s =
671 new BoxAndSolarArraySpacecraft(0, 0, 0,
672 (date, frame) -> new TimeStampedPVCoordinates(date, new Vector3D(0, 1e6, 0), Vector3D.ZERO),
673 20.0, Vector3D.PLUS_J, 0.0, 1.0, 0.0);
674 Field<Decimal64> field = Decimal64Field.getInstance();
675 FieldVector3D<Decimal64> normal = s.getNormal(FieldAbsoluteDate.getJ2000Epoch(field),
676 FramesFactory.getEME2000(),
677 FieldVector3D.getZero(field),
678 FieldRotation.getIdentity(field));
679 Assert.assertEquals(0, FieldVector3D.dotProduct(normal, Vector3D.PLUS_J).getReal(), 1.0e-16);
680 }
681
682
683
684
685 @Test
686 public void testFieldAcceleration() {
687
688 AbsoluteDate initialDate = propagator.getInitialState().getDate();
689 CelestialBody sun = CelestialBodyFactory.getSun();
690
691
692 Vector3D earthRot = new Vector3D(0., 0., Constants.GRIM5C1_EARTH_ANGULAR_VELOCITY);
693 double density = 1.e-3;
694 double refFlux = 4.56e-6;
695
696
697
698
699 BoxAndSolarArraySpacecraft s =
700 new BoxAndSolarArraySpacecraft(1., 2., 3., sun, 20.0, Vector3D.PLUS_J, 2.0, 0.3, 0.5, 0.4);
701
702 for (double dt = 0; dt < 4000; dt += 60) {
703 AbsoluteDate date = initialDate.shiftedBy(dt);
704 SpacecraftState state = propagator.propagate(date);
705
706
707 Vector3D position = state.getPVCoordinates().getPosition();
708 Vector3D velocity = state.getPVCoordinates().getVelocity();
709 Vector3D vAtm = Vector3D.crossProduct(earthRot, position);
710 Vector3D relativeVelocity = vAtm.subtract(velocity);
711
712 Frame frame = state.getFrame();
713 Rotation rotation = state.getAttitude().getRotation();
714 double mass = state.getMass();
715 Vector3D flux = position.subtract(sun.getPVCoordinates(date, frame).getPosition()).normalize().scalarMultiply(refFlux);
716
717
718 Vector3D aDrag = s.dragAcceleration(date, frame, position, rotation, mass,
719 density, relativeVelocity,
720 getDragParameters(s));
721 Vector3D aSrp = s.radiationPressureAcceleration(date, frame, position, rotation, mass,
722 flux, getRadiationParameters(s));
723
724
725 Field<Decimal64> field = Decimal64Field.getInstance();
726
727 FieldAbsoluteDate<Decimal64> dateF = new FieldAbsoluteDate<>(field, date);
728 FieldVector3D<Decimal64> positionF = new FieldVector3D<Decimal64>(field.getOne(), position);
729 FieldRotation<Decimal64> rotationF = new FieldRotation<>(field, rotation);
730 Decimal64 massF = new Decimal64(mass);
731 FieldVector3D<Decimal64> fluxF = new FieldVector3D<Decimal64>(field.getOne(), flux);
732 Decimal64 densityF = new Decimal64(density);
733 FieldVector3D<Decimal64> relativeVelocityF = new FieldVector3D<Decimal64>(field.getOne(), relativeVelocity);
734
735
736
737 FieldVector3D<Decimal64> aDragF = s.dragAcceleration(dateF, frame,
738 positionF, rotationF, massF, densityF,
739 relativeVelocityF, getDragParameters(s, field));
740 FieldVector3D<Decimal64> aSrpF = s.radiationPressureAcceleration(dateF, frame,
741 positionF, rotationF, massF,
742 fluxF, getRadiationParameters(s, field));
743
744 Assert.assertEquals(0.0, Vector3D.distance(aDrag, aDragF.toVector3D()), Precision.EPSILON);
745 Assert.assertEquals(0.0, Vector3D.distance(aSrp, aSrpF.toVector3D()), Precision.EPSILON);
746 }
747 }
748
749
750 private double[] getDragParameters(final BoxAndSolarArraySpacecraft basa) {
751 final List<ParameterDriver> drivers = basa.getDragParametersDrivers();
752 final double[] parameters = new double[drivers.size()];
753 for (int i = 0; i < drivers.size(); ++i) {
754 parameters[i] = drivers.get(i).getValue();
755 }
756 return parameters;
757 }
758
759
760 private double[] getRadiationParameters(final BoxAndSolarArraySpacecraft basa) {
761 final List<ParameterDriver> drivers = basa.getRadiationParametersDrivers();
762 final double[] parameters = new double[drivers.size()];
763 for (int i = 0; i < drivers.size(); ++i) {
764 parameters[i] = drivers.get(i).getValue();
765 }
766 return parameters;
767 }
768
769
770 private <T extends CalculusFieldElement<T>> T[] getDragParameters(final BoxAndSolarArraySpacecraft basa,
771 final Field<T> field) {
772 final List<ParameterDriver> drivers = basa.getDragParametersDrivers();
773 final T[] parameters = MathArrays.buildArray(field, drivers.size());
774 for (int i = 0; i < drivers.size(); ++i) {
775 parameters[i] = field.getZero().add(drivers.get(i).getValue());
776 }
777 return parameters;
778 }
779
780
781 private <T extends CalculusFieldElement<T>> T[] getRadiationParameters(final BoxAndSolarArraySpacecraft basa,
782 final Field<T> field) {
783 final List<ParameterDriver> drivers = basa.getRadiationParametersDrivers();
784 final T[] parameters = MathArrays.buildArray(field, drivers.size());
785 for (int i = 0; i < drivers.size(); ++i) {
786 parameters[i] = field.getZero().add(drivers.get(i).getValue());
787 }
788 return parameters;
789 }
790
791 @Before
792 public void setUp() {
793 try {
794 Utils.setDataRoot("regular-data");
795 mu = 3.9860047e14;
796 double ae = 6.378137e6;
797 double c20 = -1.08263e-3;
798 double c30 = 2.54e-6;
799 double c40 = 1.62e-6;
800 double c50 = 2.3e-7;
801 double c60 = -5.5e-7;
802
803 AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 7, 1),
804 new TimeComponents(13, 59, 27.816),
805 TimeScalesFactory.getUTC());
806
807
808
809 Orbit circ =
810 new CircularOrbit(7178000.0, 0.5e-4, -0.5e-4, FastMath.toRadians(50.), FastMath.toRadians(280),
811 FastMath.toRadians(10.0), PositionAngle.MEAN,
812 FramesFactory.getEME2000(), date, mu);
813 propagator =
814 new EcksteinHechlerPropagator(circ,
815 new LofOffset(circ.getFrame(), LOFType.LVLH_CCSDS),
816 ae, mu, c20, c30, c40, c50, c60);
817 } catch (OrekitException oe) {
818 Assert.fail(oe.getLocalizedMessage());
819 }
820 }
821
822 private double mu;
823 private Propagator propagator;
824
825 }