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