1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.attitudes;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.Field;
21 import org.hipparchus.complex.Complex;
22 import org.hipparchus.complex.ComplexField;
23 import org.hipparchus.fitting.PolynomialCurveFitter;
24 import org.hipparchus.fitting.WeightedObservedPoint;
25 import org.hipparchus.geometry.euclidean.threed.FieldRotation;
26 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
27 import org.hipparchus.geometry.euclidean.threed.Line;
28 import org.hipparchus.geometry.euclidean.threed.Rotation;
29 import org.hipparchus.geometry.euclidean.threed.Vector3D;
30 import org.hipparchus.util.Binary64Field;
31 import org.hipparchus.util.FastMath;
32 import org.hipparchus.util.MathArrays;
33 import org.junit.jupiter.api.AfterEach;
34 import org.junit.jupiter.api.Assertions;
35 import org.junit.jupiter.api.BeforeEach;
36 import org.junit.jupiter.api.Test;
37 import org.orekit.Utils;
38 import org.orekit.annotation.DefaultDataContext;
39 import org.orekit.bodies.GeodeticPoint;
40 import org.orekit.bodies.OneAxisEllipsoid;
41 import org.orekit.errors.OrekitException;
42 import org.orekit.frames.FramesFactory;
43 import org.orekit.frames.StaticTransform;
44 import org.orekit.frames.Transform;
45 import org.orekit.orbits.CircularOrbit;
46 import org.orekit.orbits.FieldCircularOrbit;
47 import org.orekit.orbits.FieldKeplerianOrbit;
48 import org.orekit.orbits.PositionAngleType;
49 import org.orekit.propagation.FieldSpacecraftState;
50 import org.orekit.propagation.SpacecraftState;
51 import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
52 import org.orekit.propagation.analytical.FieldEcksteinHechlerPropagator;
53 import org.orekit.time.AbsoluteDate;
54 import org.orekit.time.DateComponents;
55 import org.orekit.time.FieldAbsoluteDate;
56 import org.orekit.time.TimeComponents;
57 import org.orekit.time.TimeScalesFactory;
58 import org.orekit.utils.AngularCoordinates;
59 import org.orekit.utils.Constants;
60 import org.orekit.utils.FieldAngularCoordinates;
61 import org.orekit.utils.FieldPVCoordinates;
62 import org.orekit.utils.IERSConventions;
63 import org.orekit.utils.PVCoordinates;
64 import org.orekit.utils.TimeStampedFieldPVCoordinates;
65 import org.orekit.utils.TimeStampedPVCoordinates;
66
67 import java.util.ArrayList;
68 import java.util.List;
69
70
71 class BodyCenterPointingTest {
72
73
74 private AbsoluteDate date;
75
76
77 private CircularOrbit circ;
78
79
80 private OneAxisEllipsoid earth;
81
82
83 private Transform eme2000ToItrf;
84
85
86 private BodyCenterPointing earthCenterAttitudeLaw;
87
88 @Test
89 void testGetPosition() {
90
91
92 final Vector3D actualPosition = earthCenterAttitudeLaw.getTargetPosition(circ, date, circ.getFrame());
93
94 final Vector3D expectedPosition = earthCenterAttitudeLaw.getTargetPV(circ, date, circ.getFrame()).getPosition();
95 Assertions.assertEquals(expectedPosition, actualPosition);
96 }
97
98 @Test
99 void testGetPositionField() {
100
101 final FieldCircularOrbit<Complex> fieldCircularOrbit = new FieldCircularOrbit<>(ComplexField.getInstance(), circ);
102
103 final FieldVector3D<Complex> actualPosition = earthCenterAttitudeLaw.getTargetPosition(fieldCircularOrbit, fieldCircularOrbit.getDate(), circ.getFrame());
104
105 final FieldVector3D<Complex> expectedPosition = earthCenterAttitudeLaw.getTargetPV(fieldCircularOrbit, fieldCircularOrbit.getDate(), circ.getFrame()).getPosition();
106 Assertions.assertEquals(0., expectedPosition.subtract(actualPosition).getNorm().getReal(), 1e-10);
107 }
108
109
110
111 @Test
112 void testTarget() {
113
114
115 TimeStampedPVCoordinates target = earthCenterAttitudeLaw.getTargetPV(circ, date, circ.getFrame());
116
117
118 GeodeticPoint gp = earth.transform(target.getPosition(), circ.getFrame(), date);
119 Assertions.assertEquals(0.0, gp.getAltitude(), 1.0e-10);
120 Assertions.assertEquals(date, target.getDate());
121
122 }
123
124
125
126 @Test
127 void testBodyCenterInPointingDirection() {
128
129
130 PVCoordinates pvSatEME2000 = circ.getPVCoordinates();
131
132
133
134
135 Rotation rotSatEME2000 = earthCenterAttitudeLaw.getAttitude(circ, date, circ.getFrame()).getRotation();
136
137
138 Vector3D zSatEME2000 = rotSatEME2000.applyInverseTo(Vector3D.PLUS_K);
139
140
141 Vector3D zSatITRF2008C = eme2000ToItrf.transformVector(zSatEME2000);
142
143
144 PVCoordinates pvSatITRF2008C = eme2000ToItrf.transformPVCoordinates(pvSatEME2000);
145
146
147 Line pointingLine = new Line(pvSatITRF2008C.getPosition(),
148 pvSatITRF2008C.getPosition().add(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
149 zSatITRF2008C),
150 2.0e-8);
151
152
153 Assertions.assertTrue(pointingLine.contains(Vector3D.ZERO));
154
155 }
156
157 @Test
158 @DefaultDataContext
159 void testQDot() {
160
161 Utils.setDataRoot("regular-data");
162 final double ehMu = 3.9860047e14;
163 final double ae = 6.378137e6;
164 final double c20 = -1.08263e-3;
165 final double c30 = 2.54e-6;
166 final double c40 = 1.62e-6;
167 final double c50 = 2.3e-7;
168 final double c60 = -5.5e-7;
169 final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
170 final Vector3D position = new Vector3D(3220103., 69623., 6449822.);
171 final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
172 final CircularOrbit initialOrbit = new CircularOrbit(new PVCoordinates(position, velocity),
173 FramesFactory.getEME2000(), date, ehMu);
174
175 EcksteinHechlerPropagator propagator =
176 new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
177 propagator.setAttitudeProvider(earthCenterAttitudeLaw);
178
179 List<WeightedObservedPoint> w0 = new ArrayList<>();
180 List<WeightedObservedPoint> w1 = new ArrayList<>();
181 List<WeightedObservedPoint> w2 = new ArrayList<>();
182 List<WeightedObservedPoint> w3 = new ArrayList<>();
183 for (double dt = -1; dt < 1; dt += 0.01) {
184 Rotation rP = propagator.propagate(date.shiftedBy(dt)).getAttitude().getRotation();
185 w0.add(new WeightedObservedPoint(1, dt, rP.getQ0()));
186 w1.add(new WeightedObservedPoint(1, dt, rP.getQ1()));
187 w2.add(new WeightedObservedPoint(1, dt, rP.getQ2()));
188 w3.add(new WeightedObservedPoint(1, dt, rP.getQ3()));
189 }
190
191 double q0DotRef = PolynomialCurveFitter.create(2).fit(w0)[1];
192 double q1DotRef = PolynomialCurveFitter.create(2).fit(w1)[1];
193 double q2DotRef = PolynomialCurveFitter.create(2).fit(w2)[1];
194 double q3DotRef = PolynomialCurveFitter.create(2).fit(w3)[1];
195
196 Attitude a0 = propagator.propagate(date).getAttitude();
197 double q0 = a0.getRotation().getQ0();
198 double q1 = a0.getRotation().getQ1();
199 double q2 = a0.getRotation().getQ2();
200 double q3 = a0.getRotation().getQ3();
201 double oX = a0.getSpin().getX();
202 double oY = a0.getSpin().getY();
203 double oZ = a0.getSpin().getZ();
204
205
206 double q0Dot = 0.5 * MathArrays.linearCombination(-q1, oX, -q2, oY, -q3, oZ);
207 double q1Dot = 0.5 * MathArrays.linearCombination( q0, oX, -q3, oY, q2, oZ);
208 double q2Dot = 0.5 * MathArrays.linearCombination( q3, oX, q0, oY, -q1, oZ);
209 double q3Dot = 0.5 * MathArrays.linearCombination(-q2, oX, q1, oY, q0, oZ);
210
211 Assertions.assertEquals(q0DotRef, q0Dot, 5.0e-9);
212 Assertions.assertEquals(q1DotRef, q1Dot, 5.0e-9);
213 Assertions.assertEquals(q2DotRef, q2Dot, 5.0e-9);
214 Assertions.assertEquals(q3DotRef, q3Dot, 5.0e-9);
215
216 }
217
218 @Test
219 @DefaultDataContext
220 void testSpin() {
221
222 Utils.setDataRoot("regular-data");
223 final double ehMu = 3.9860047e14;
224 final double ae = 6.378137e6;
225 final double c20 = -1.08263e-3;
226 final double c30 = 2.54e-6;
227 final double c40 = 1.62e-6;
228 final double c50 = 2.3e-7;
229 final double c60 = -5.5e-7;
230 final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
231 final Vector3D position = new Vector3D(3220103., 69623., 6449822.);
232 final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
233 final CircularOrbit initialOrbit = new CircularOrbit(new PVCoordinates(position, velocity),
234 FramesFactory.getEME2000(), date, ehMu);
235
236 EcksteinHechlerPropagator propagator =
237 new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
238 propagator.setAttitudeProvider(earthCenterAttitudeLaw);
239
240 double h = 0.01;
241 SpacecraftState s0 = propagator.propagate(date);
242 SpacecraftState sMinus = propagator.propagate(date.shiftedBy(-h));
243 SpacecraftState sPlus = propagator.propagate(date.shiftedBy(h));
244
245
246 double errorAngleMinus = Rotation.distance(sMinus.shiftedBy(h).getAttitude().getRotation(),
247 s0.getAttitude().getRotation());
248 double evolutionAngleMinus = Rotation.distance(sMinus.getAttitude().getRotation(),
249 s0.getAttitude().getRotation());
250 Assertions.assertEquals(0.0, errorAngleMinus, 1.0e-6 * evolutionAngleMinus);
251 double errorAnglePlus = Rotation.distance(s0.getAttitude().getRotation(),
252 sPlus.shiftedBy(-h).getAttitude().getRotation());
253 double evolutionAnglePlus = Rotation.distance(s0.getAttitude().getRotation(),
254 sPlus.getAttitude().getRotation());
255 Assertions.assertEquals(0.0, errorAnglePlus, 1.0e-6 * evolutionAnglePlus);
256
257 Vector3D spin0 = s0.getAttitude().getSpin();
258 Vector3D reference = AngularCoordinates.estimateRate(sMinus.getAttitude().getRotation(),
259 sPlus.getAttitude().getRotation(),
260 2 * h);
261 Assertions.assertTrue(spin0.getNorm() > 1.0e-3);
262 Assertions.assertEquals(0.0, spin0.subtract(reference).getNorm(), 1.0e-13);
263
264 }
265
266 @Test
267 @DefaultDataContext
268 void testTargetField() {
269 doTestTarget(Binary64Field.getInstance());
270 }
271
272 @Test
273 @DefaultDataContext
274 public void doxBodyCenterInPointingDirectionTest() {
275 doTestBodyCenterInPointingDirection(Binary64Field.getInstance());
276 }
277
278 @Test
279 @DefaultDataContext
280 void testQDotField() {
281 doTestQDot(Binary64Field.getInstance());
282 }
283
284 @Test
285 @DefaultDataContext
286 void testSpinField() {
287 doTestSpin(Binary64Field.getInstance());
288 }
289
290 @DefaultDataContext
291 private <T extends CalculusFieldElement<T>>void doTestTarget(final Field<T> field) {
292
293 T mu = field.getZero().add(3.9860047e14);
294 T zero = field.getZero();
295
296 final T raan = zero.add(FastMath.toRadians(270.));
297 final T a =zero.add(7178000.0);
298 final T e =zero.add(7e-5);
299 final T i =zero.add(FastMath.toRadians(50.));
300 final T pa=zero.add(FastMath.toRadians(45.));
301 final T m =zero.add(FastMath.toRadians(5.3-270));
302
303
304 FieldAbsoluteDate<T> date= new FieldAbsoluteDate<>(field, new DateComponents(2008, 4, 7),
305 TimeComponents.H00,
306 TimeScalesFactory.getUTC());
307
308 FieldKeplerianOrbit<T> circ = new FieldKeplerianOrbit<>(a, e, i, pa, raan,
309 m, PositionAngleType.MEAN,
310 FramesFactory.getEME2000(), date, mu);
311
312 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
313 Constants.WGS84_EARTH_FLATTENING,
314 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
315
316 BodyCenterPointing earthCenterAttitudeLaw= new BodyCenterPointing(circ.getFrame(), earth);
317
318
319 TimeStampedFieldPVCoordinates<T> target = earthCenterAttitudeLaw.getTargetPV(circ, date, circ.getFrame());
320
321
322 GeodeticPoint gp = earth.transform(target.getPosition().toVector3D(), circ.getFrame(), date.toAbsoluteDate());
323
324 Assertions.assertEquals(0.0, gp.getAltitude(), 1.0e-8);
325 Assertions.assertEquals(date, target.getDate());
326
327 }
328
329 @DefaultDataContext
330 private <T extends CalculusFieldElement<T>> void doTestBodyCenterInPointingDirection(final Field<T> field) {
331
332 T zero = field.getZero();
333 T mu = zero.add(3.9860047e14);
334
335 final T raan = zero.add(FastMath.toRadians(270.));
336 final T a =zero.add(7178000.0);
337 final T e =zero.add(7E-5);
338 final T i =zero.add(FastMath.toRadians(50.));
339 final T pa=zero.add(FastMath.toRadians(45.));
340
341 final T m =zero.add(FastMath.toRadians(5.300-270.));
342
343
344 FieldAbsoluteDate<T> date= new FieldAbsoluteDate<>(field, new DateComponents(2008, 4, 7),
345 TimeComponents.H00,
346 TimeScalesFactory.getUTC());
347
348 FieldKeplerianOrbit<T> circ = new FieldKeplerianOrbit<>(a, e, i, pa, raan,
349 m, PositionAngleType.MEAN,
350 FramesFactory.getEME2000(), date, mu);
351
352
353 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
354 Constants.WGS84_EARTH_FLATTENING,
355 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
356
357 StaticTransform eme2000ToItrf = FramesFactory.getEME2000().getStaticTransformTo(earth.getBodyFrame(), date.toAbsoluteDate());
358
359
360 BodyCenterPointing earthCenterAttitudeLaw= new BodyCenterPointing(circ.getFrame(), earth);
361
362 FieldVector3D<T> positionSatEME2000 = circ.getPosition();
363
364
365
366
367 FieldRotation<T> rotSatEME2000 = earthCenterAttitudeLaw.getAttitude(circ, date, circ.getFrame()).getRotation();
368
369
370
371 FieldVector3D<T> zSatEME2000 = rotSatEME2000.applyInverseTo(Vector3D.PLUS_K);
372
373
374 FieldVector3D<T> zSatITRF2008C = eme2000ToItrf.transformVector(zSatEME2000);
375
376
377 FieldVector3D<T> positionSatITRF2008C = eme2000ToItrf.transformPosition(positionSatEME2000);
378
379
380
381 Line pointingLine = new Line(positionSatITRF2008C.toVector3D(),
382 positionSatITRF2008C.toVector3D().add(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
383 zSatITRF2008C.toVector3D()),
384 2.0e-8);
385
386
387 Assertions.assertTrue(pointingLine.contains(Vector3D.ZERO));
388
389 }
390
391 @DefaultDataContext
392 private <T extends CalculusFieldElement<T>> void doTestQDot(final Field<T> field) {
393
394 final double ae = 6.378137e6;
395 final double c20 = -1.08263e-3;
396 final double c30 = 2.54e-6;
397 final double c40 = 1.62e-6;
398 final double c50 = 2.3e-7;
399 final double c60 = -5.5e-7;
400
401
402 T zero = field.getZero();
403 final T a = zero.add(7178000.0);
404 final T e = zero.add(7e-5);
405 final T i = zero.add(FastMath.toRadians(50.));
406 final T pa = zero.add(FastMath.toRadians(45.));
407 final T raan = zero.add(FastMath.toRadians(270.));
408 final T m = zero.add(FastMath.toRadians(5.3-270));
409 final T ehMu = zero.add(3.9860047e14);
410
411
412 FieldAbsoluteDate<T> date_comp= new FieldAbsoluteDate<>(field, new DateComponents(2008, 4, 7),
413 TimeComponents.H00,
414 TimeScalesFactory.getUTC());
415
416 FieldKeplerianOrbit<T> circ = new FieldKeplerianOrbit<>(a, e, i, pa, raan,
417 m, PositionAngleType.MEAN,
418 FramesFactory.getEME2000(), date_comp, ehMu);
419
420 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
421 Constants.WGS84_EARTH_FLATTENING,
422 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
423
424
425 BodyCenterPointing earthCenterAttitudeLaw= new BodyCenterPointing(circ.getFrame(), earth);
426
427
428 final FieldAbsoluteDate<T> date = FieldAbsoluteDate.getJ2000Epoch(field).shiftedBy(584.);
429 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(6449822.));
430 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
431 final FieldCircularOrbit<T> initialOrbit = new FieldCircularOrbit<>(new FieldPVCoordinates<>(position, velocity),
432 FramesFactory.getEME2000(), date, ehMu);
433
434 FieldEcksteinHechlerPropagator<T> propagator =
435 new FieldEcksteinHechlerPropagator<>(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
436 propagator.setAttitudeProvider(earthCenterAttitudeLaw);
437
438 List<WeightedObservedPoint> w0 = new ArrayList<>();
439 List<WeightedObservedPoint> w1 = new ArrayList<>();
440 List<WeightedObservedPoint> w2 = new ArrayList<>();
441 List<WeightedObservedPoint> w3 = new ArrayList<>();
442 for (double dt = -1; dt < 1; dt += 0.01) {
443 FieldRotation<T> rP = propagator.propagate(date.shiftedBy(dt)).getAttitude().getRotation();
444 w0.add(new WeightedObservedPoint(1, dt, rP.getQ0().getReal()));
445 w1.add(new WeightedObservedPoint(1, dt, rP.getQ1().getReal()));
446 w2.add(new WeightedObservedPoint(1, dt, rP.getQ2().getReal()));
447 w3.add(new WeightedObservedPoint(1, dt, rP.getQ3().getReal()));
448 }
449
450 double q0DotRef = PolynomialCurveFitter.create(2).fit(w0)[1];
451 double q1DotRef = PolynomialCurveFitter.create(2).fit(w1)[1];
452 double q2DotRef = PolynomialCurveFitter.create(2).fit(w2)[1];
453 double q3DotRef = PolynomialCurveFitter.create(2).fit(w3)[1];
454
455 FieldAttitude<T> a0 = propagator.propagate(date).getAttitude();
456 T q0 = a0.getRotation().getQ0();
457 T q1 = a0.getRotation().getQ1();
458 T q2 = a0.getRotation().getQ2();
459 T q3 = a0.getRotation().getQ3();
460 T oX = a0.getSpin().getX();
461 T oY = a0.getSpin().getY();
462 T oZ = a0.getSpin().getZ();
463
464
465 double q0Dot = 0.5 * MathArrays.linearCombination(-q1.getReal(), oX.getReal(), -q2.getReal(), oY.getReal(), -q3.getReal(), oZ.getReal());
466 double q1Dot = 0.5 * MathArrays.linearCombination( q0.getReal(), oX.getReal(), -q3.getReal(), oY.getReal(), q2.getReal(), oZ.getReal());
467 double q2Dot = 0.5 * MathArrays.linearCombination( q3.getReal(), oX.getReal(), q0.getReal(), oY.getReal(), -q1.getReal(), oZ.getReal());
468 double q3Dot = 0.5 * MathArrays.linearCombination(-q2.getReal(), oX.getReal(), q1.getReal(), oY.getReal(), q0.getReal(), oZ.getReal());
469
470 Assertions.assertEquals(q0DotRef, q0Dot, 5.0e-9);
471 Assertions.assertEquals(q1DotRef, q1Dot, 5.0e-9);
472 Assertions.assertEquals(q2DotRef, q2Dot, 5.0e-9);
473 Assertions.assertEquals(q3DotRef, q3Dot, 5.0e-9);
474
475 }
476
477 @DefaultDataContext
478 private <T extends CalculusFieldElement<T>> void doTestSpin(final Field<T> field) {
479
480 final double ae = 6.378137e6;
481 final double c20 = -1.08263e-3;
482 final double c30 = 2.54e-6;
483 final double c40 = 1.62e-6;
484 final double c50 = 2.3e-7;
485 final double c60 = -5.5e-7;
486
487
488 final T zero = field.getZero();
489 final T a = zero.add(7178000.0);
490 final T e = zero.add(7e-5);
491 final T i = zero.add(FastMath.toRadians(50.));
492 final T pa = zero.add(FastMath.toRadians(45.));
493 final T raan = zero.add(FastMath.toRadians(270.));
494 final T m = zero.add(FastMath.toRadians(5.3-270));
495 final T ehMu = zero.add(3.9860047e14);
496
497
498 FieldAbsoluteDate<T> date_R = new FieldAbsoluteDate<>(field, new DateComponents(2008, 4, 7),
499 TimeComponents.H00,
500 TimeScalesFactory.getUTC());
501
502 FieldKeplerianOrbit<T> circ = new FieldKeplerianOrbit<>(a, e, i, pa, raan,
503 m, PositionAngleType.MEAN,
504 FramesFactory.getEME2000(), date_R, ehMu);
505
506 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
507 Constants.WGS84_EARTH_FLATTENING,
508 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
509
510
511 BodyCenterPointing earthCenterAttitudeLaw= new BodyCenterPointing(circ.getFrame(), earth);
512
513
514 final FieldAbsoluteDate<T> date = FieldAbsoluteDate.getJ2000Epoch(field).shiftedBy(584.);
515 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(6449822.));
516 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
517 final FieldCircularOrbit<T> initialOrbit = new FieldCircularOrbit<>(new FieldPVCoordinates<>(position, velocity),
518 FramesFactory.getEME2000(), date, ehMu);
519
520 FieldEcksteinHechlerPropagator<T> propagator =
521 new FieldEcksteinHechlerPropagator<>(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
522 propagator.setAttitudeProvider(earthCenterAttitudeLaw);
523
524 double h = 0.01;
525
526 FieldSpacecraftState<T> s0 = propagator.propagate(date);
527 FieldSpacecraftState<T> sMinus = propagator.propagate(date.shiftedBy(-h));
528 FieldSpacecraftState<T> sPlus = propagator.propagate(date.shiftedBy(h));
529
530
531 T errorAngleMinus = FieldRotation.distance(sMinus.shiftedBy(zero.add(h)).getAttitude().getRotation(),
532 s0.getAttitude().getRotation());
533 T evolutionAngleMinus = FieldRotation.distance(sMinus.getAttitude().getRotation(),
534 s0.getAttitude().getRotation());
535 Assertions.assertEquals(0.0, errorAngleMinus.getReal(), 1.0e-6 * evolutionAngleMinus.getReal());
536 T errorAnglePlus = FieldRotation.distance(s0.getAttitude().getRotation(),
537 sPlus.shiftedBy(zero.add(-h)).getAttitude().getRotation());
538 T evolutionAnglePlus = FieldRotation.distance(s0.getAttitude().getRotation(),
539 sPlus.getAttitude().getRotation());
540 Assertions.assertEquals(0.0, errorAnglePlus.getReal(), 1.0e-6 * evolutionAnglePlus.getReal());
541
542 FieldVector3D<T> spin0 = s0.getAttitude().getSpin();
543 FieldVector3D<T> reference = FieldAngularCoordinates.estimateRate(sMinus.getAttitude().getRotation(),
544 sPlus.getAttitude().getRotation(),
545 2 * h);
546 Assertions.assertTrue(spin0.getNorm().getReal() > 1.0e-3);
547 Assertions.assertEquals(0.0, spin0.subtract(reference).getNorm().getReal(), 1.4e-13);
548
549 }
550
551 @BeforeEach
552 @DefaultDataContext
553 public void setUp() {
554 try {
555
556 Utils.setDataRoot("regular-data");
557
558
559 date = new AbsoluteDate(new DateComponents(2008, 4, 7),
560 TimeComponents.H00,
561 TimeScalesFactory.getUTC());
562
563
564 final double mu = 3.9860047e14;
565 final double raan = 270.;
566 circ =
567 new CircularOrbit(7178000.0, 0.5e-4, -0.5e-4, FastMath.toRadians(50.), FastMath.toRadians(raan),
568 FastMath.toRadians(5.300 - raan), PositionAngleType.MEAN,
569 FramesFactory.getEME2000(), date, mu);
570
571
572
573 earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
574 Constants.WGS84_EARTH_FLATTENING,
575 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
576
577
578 eme2000ToItrf = FramesFactory.getEME2000().getTransformTo(earth.getBodyFrame(), date);
579
580
581 earthCenterAttitudeLaw = new BodyCenterPointing(circ.getFrame(), earth);
582
583 } catch (OrekitException oe) {
584 Assertions.fail(oe.getMessage());
585 }
586
587 }
588
589 @AfterEach
590 public void tearDown() {
591 date = null;
592 earth = null;
593 eme2000ToItrf = null;
594 earthCenterAttitudeLaw = null;
595 circ = null;
596 }
597
598 }
599