1   /* Copyright 2002-2022 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.attitudes;
18  
19  
20  import java.util.ArrayList;
21  import java.util.List;
22  
23  import org.hipparchus.Field;
24  import org.hipparchus.CalculusFieldElement;
25  import org.hipparchus.fitting.PolynomialCurveFitter;
26  import org.hipparchus.fitting.WeightedObservedPoint;
27  import org.hipparchus.geometry.euclidean.threed.FieldRotation;
28  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
29  import org.hipparchus.geometry.euclidean.threed.Line;
30  import org.hipparchus.geometry.euclidean.threed.Rotation;
31  import org.hipparchus.geometry.euclidean.threed.Vector3D;
32  import org.hipparchus.util.Decimal64Field;
33  import org.hipparchus.util.FastMath;
34  import org.hipparchus.util.MathArrays;
35  import org.junit.After;
36  import org.junit.Assert;
37  import org.junit.Before;
38  import org.junit.Test;
39  import org.orekit.Utils;
40  import org.orekit.bodies.GeodeticPoint;
41  import org.orekit.bodies.OneAxisEllipsoid;
42  import org.orekit.errors.OrekitException;
43  import org.orekit.frames.FramesFactory;
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.PositionAngle;
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  
68  public class BodyCenterPointingTest {
69  
70      // Computation date
71      private AbsoluteDate date;
72  
73      // Orbit
74      private CircularOrbit circ;
75  
76      // WGS84 Earth model
77      private OneAxisEllipsoid earth;
78  
79      // Transform from EME2000 to ITRF2008
80      private Transform eme2000ToItrf;
81  
82      // Earth center pointing attitude provider
83      private BodyCenterPointing earthCenterAttitudeLaw;
84  
85      /** Test if target is on Earth surface
86       */
87      @Test
88      public void testTarget() {
89  
90          // Call get target method
91          TimeStampedPVCoordinates target = earthCenterAttitudeLaw.getTargetPV(circ, date, circ.getFrame());
92  
93          // Check that target is on Earth surface
94          GeodeticPoint gp = earth.transform(target.getPosition(), circ.getFrame(), date);
95          Assert.assertEquals(0.0, gp.getAltitude(), 1.0e-10);
96          Assert.assertEquals(date, target.getDate());
97  
98      }
99  
100     /** Test if body center belongs to the direction pointed by the satellite
101      */
102     @Test
103     public void testBodyCenterInPointingDirection() {
104 
105         // Transform satellite position to position/velocity parameters in EME2000 frame
106         PVCoordinates pvSatEME2000 = circ.getPVCoordinates();
107 
108         //  Pointing direction
109         // ********************
110         // Get satellite attitude rotation, i.e rotation from EME2000 frame to satellite frame
111         Rotation rotSatEME2000 = earthCenterAttitudeLaw.getAttitude(circ, date, circ.getFrame()).getRotation();
112 
113         // Transform Z axis from satellite frame to EME2000
114         Vector3D zSatEME2000 = rotSatEME2000.applyInverseTo(Vector3D.PLUS_K);
115 
116         // Transform Z axis from EME2000 to ITRF2008
117         Vector3D zSatITRF2008C = eme2000ToItrf.transformVector(zSatEME2000);
118 
119         // Transform satellite position/velocity from EME2000 to ITRF2008
120         PVCoordinates pvSatITRF2008C = eme2000ToItrf.transformPVCoordinates(pvSatEME2000);
121 
122        // Line containing satellite point and following pointing direction
123         Line pointingLine = new Line(pvSatITRF2008C.getPosition(),
124                                      pvSatITRF2008C.getPosition().add(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
125                                                                       zSatITRF2008C),
126                                      2.0e-8);
127 
128         // Check that the line contains Earth center
129         Assert.assertTrue(pointingLine.contains(Vector3D.ZERO));
130 
131     }
132 
133     @Test
134     public void testQDot() {
135 
136         Utils.setDataRoot("regular-data");
137         final double ehMu  = 3.9860047e14;
138         final double ae  = 6.378137e6;
139         final double c20 = -1.08263e-3;
140         final double c30 = 2.54e-6;
141         final double c40 = 1.62e-6;
142         final double c50 = 2.3e-7;
143         final double c60 = -5.5e-7;
144         final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
145         final Vector3D position = new Vector3D(3220103., 69623., 6449822.);
146         final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
147         final CircularOrbit initialOrbit = new CircularOrbit(new PVCoordinates(position, velocity),
148                                                              FramesFactory.getEME2000(), date, ehMu);
149 
150         EcksteinHechlerPropagator propagator =
151                 new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
152         propagator.setAttitudeProvider(earthCenterAttitudeLaw);
153 
154         List<WeightedObservedPoint> w0 = new ArrayList<WeightedObservedPoint>();
155         List<WeightedObservedPoint> w1 = new ArrayList<WeightedObservedPoint>();
156         List<WeightedObservedPoint> w2 = new ArrayList<WeightedObservedPoint>();
157         List<WeightedObservedPoint> w3 = new ArrayList<WeightedObservedPoint>();
158         for (double dt = -1; dt < 1; dt += 0.01) {
159             Rotation rP = propagator.propagate(date.shiftedBy(dt)).getAttitude().getRotation();
160             w0.add(new WeightedObservedPoint(1, dt, rP.getQ0()));
161             w1.add(new WeightedObservedPoint(1, dt, rP.getQ1()));
162             w2.add(new WeightedObservedPoint(1, dt, rP.getQ2()));
163             w3.add(new WeightedObservedPoint(1, dt, rP.getQ3()));
164         }
165 
166         double q0DotRef = PolynomialCurveFitter.create(2).fit(w0)[1];
167         double q1DotRef = PolynomialCurveFitter.create(2).fit(w1)[1];
168         double q2DotRef = PolynomialCurveFitter.create(2).fit(w2)[1];
169         double q3DotRef = PolynomialCurveFitter.create(2).fit(w3)[1];
170 
171         Attitude a0 = propagator.propagate(date).getAttitude();
172         double   q0 = a0.getRotation().getQ0();
173         double   q1 = a0.getRotation().getQ1();
174         double   q2 = a0.getRotation().getQ2();
175         double   q3 = a0.getRotation().getQ3();
176         double   oX = a0.getSpin().getX();
177         double   oY = a0.getSpin().getY();
178         double   oZ = a0.getSpin().getZ();
179 
180         // first time-derivatives of the quaternion
181         double q0Dot = 0.5 * MathArrays.linearCombination(-q1, oX, -q2, oY, -q3, oZ);
182         double q1Dot = 0.5 * MathArrays.linearCombination( q0, oX, -q3, oY,  q2, oZ);
183         double q2Dot = 0.5 * MathArrays.linearCombination( q3, oX,  q0, oY, -q1, oZ);
184         double q3Dot = 0.5 * MathArrays.linearCombination(-q2, oX,  q1, oY,  q0, oZ);
185 
186         Assert.assertEquals(q0DotRef, q0Dot, 5.0e-9);
187         Assert.assertEquals(q1DotRef, q1Dot, 5.0e-9);
188         Assert.assertEquals(q2DotRef, q2Dot, 5.0e-9);
189         Assert.assertEquals(q3DotRef, q3Dot, 5.0e-9);
190 
191     }
192 
193     @Test
194     public void testSpin() {
195 
196         Utils.setDataRoot("regular-data");
197         final double ehMu  = 3.9860047e14;
198         final double ae  = 6.378137e6;
199         final double c20 = -1.08263e-3;
200         final double c30 = 2.54e-6;
201         final double c40 = 1.62e-6;
202         final double c50 = 2.3e-7;
203         final double c60 = -5.5e-7;
204         final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
205         final Vector3D position = new Vector3D(3220103., 69623., 6449822.);
206         final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
207         final CircularOrbit initialOrbit = new CircularOrbit(new PVCoordinates(position, velocity),
208                                                              FramesFactory.getEME2000(), date, ehMu);
209 
210         EcksteinHechlerPropagator propagator =
211                 new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
212         propagator.setAttitudeProvider(earthCenterAttitudeLaw);
213 
214         double h = 0.01;
215         SpacecraftState s0     = propagator.propagate(date);
216         SpacecraftState sMinus = propagator.propagate(date.shiftedBy(-h));
217         SpacecraftState sPlus  = propagator.propagate(date.shiftedBy(h));
218 
219         // check spin is consistent with attitude evolution
220         double errorAngleMinus     = Rotation.distance(sMinus.shiftedBy(h).getAttitude().getRotation(),
221                                                        s0.getAttitude().getRotation());
222         double evolutionAngleMinus = Rotation.distance(sMinus.getAttitude().getRotation(),
223                                                        s0.getAttitude().getRotation());
224         Assert.assertEquals(0.0, errorAngleMinus, 1.0e-6 * evolutionAngleMinus);
225         double errorAnglePlus      = Rotation.distance(s0.getAttitude().getRotation(),
226                                                        sPlus.shiftedBy(-h).getAttitude().getRotation());
227         double evolutionAnglePlus  = Rotation.distance(s0.getAttitude().getRotation(),
228                                                        sPlus.getAttitude().getRotation());
229         Assert.assertEquals(0.0, errorAnglePlus, 1.0e-6 * evolutionAnglePlus);
230 
231         Vector3D spin0 = s0.getAttitude().getSpin();
232         Vector3D reference = AngularCoordinates.estimateRate(sMinus.getAttitude().getRotation(),
233                                                              sPlus.getAttitude().getRotation(),
234                                                              2 * h);
235         Assert.assertTrue(spin0.getNorm() > 1.0e-3);
236         Assert.assertEquals(0.0, spin0.subtract(reference).getNorm(), 1.0e-13);
237 
238     }
239 
240     @Test
241     public void testTargetField() {
242         doTestTarget(Decimal64Field.getInstance());
243     }
244     @Test
245     public void doxBodyCenterInPointingDirectionTest() {
246         doTestBodyCenterInPointingDirection(Decimal64Field.getInstance());
247     }
248 
249     @Test
250     public void testQDotField() {
251         doTestQDot(Decimal64Field.getInstance());
252     }
253 
254     @Test
255     public void testSpinField() {
256         doTestSpin(Decimal64Field.getInstance());
257     }
258 
259     private <T extends CalculusFieldElement<T>>void doTestTarget(final Field<T> field) {
260 
261         T mu = field.getZero().add(3.9860047e14);
262         T zero = field.getZero();
263         // Satellite position as circular parameters
264         final T raan = zero.add(FastMath.toRadians(270.));
265         final T a =zero.add(7178000.0);
266         final T e =zero.add(7e-5);
267         final T i =zero.add(FastMath.toRadians(50.));
268         final T pa=zero.add(FastMath.toRadians(45.));
269         final T m =zero.add(FastMath.toRadians(5.3-270));
270 
271      // Computation date
272         FieldAbsoluteDate<T> date= new FieldAbsoluteDate<>(field, new DateComponents(2008, 04, 07),
273                                                            TimeComponents.H00,
274                                                            TimeScalesFactory.getUTC());
275         // Orbit
276         FieldKeplerianOrbit<T> circ = new FieldKeplerianOrbit<>(a, e, i, pa, raan,
277                                                                 m, PositionAngle.MEAN,
278                                                                 FramesFactory.getEME2000(), date, mu);
279         // WGS84 Earth model
280         OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
281                                                       Constants.WGS84_EARTH_FLATTENING,
282                                                       FramesFactory.getITRF(IERSConventions.IERS_2010, true));
283         // Earth center pointing attitude provider
284         BodyCenterPointing earthCenterAttitudeLaw= new BodyCenterPointing(circ.getFrame(), earth);
285 
286         // Call get target method
287         TimeStampedFieldPVCoordinates<T> target = earthCenterAttitudeLaw.getTargetPV(circ, date, circ.getFrame());
288 
289         // Check that target is on Earth surface
290         GeodeticPoint gp = earth.transform(target.getPosition().toVector3D(), circ.getFrame(), date.toAbsoluteDate());
291 
292         Assert.assertEquals(0.0, gp.getAltitude(), 1.0e-8); //less precision because i suppose we are working with Keplerian instead of circular
293         Assert.assertEquals(date, target.getDate());
294 
295     }
296 
297     private <T extends CalculusFieldElement<T>> void doTestBodyCenterInPointingDirection(final Field<T> field)  {
298 
299         T zero = field.getZero();
300         T mu = zero.add(3.9860047e14);
301         // Satellite position as circular parameters
302         final T raan = zero.add(FastMath.toRadians(270.));
303         final T a =zero.add(7178000.0);
304         final T e =zero.add(7E-5);
305         final T i =zero.add(FastMath.toRadians(50.));
306         final T pa=zero.add(FastMath.toRadians(45.));
307 
308         final T m =zero.add(FastMath.toRadians(5.300-270.));
309 
310      // Computation date
311         FieldAbsoluteDate<T> date= new FieldAbsoluteDate<>(field, new DateComponents(2008, 04, 07),
312                                                            TimeComponents.H00,
313                                                            TimeScalesFactory.getUTC());
314         // Orbit
315         FieldKeplerianOrbit<T> circ = new FieldKeplerianOrbit<>(a, e, i, pa, raan,
316                                                                 m, PositionAngle.MEAN,
317                                                                 FramesFactory.getEME2000(), date, mu);
318 
319         // WGS84 Earth model
320         OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
321                                                       Constants.WGS84_EARTH_FLATTENING,
322                                                       FramesFactory.getITRF(IERSConventions.IERS_2010, true));
323         // Transform from EME2000 to ITRF2008
324         Transform eme2000ToItrf = FramesFactory.getEME2000().getTransformTo(earth.getBodyFrame(), date.toAbsoluteDate());
325 
326         // Earth center pointing attitude provider
327         BodyCenterPointing earthCenterAttitudeLaw= new BodyCenterPointing(circ.getFrame(), earth);
328         // Transform satellite position to position/velocity parameters in EME2000 frame
329         FieldPVCoordinates<T> pvSatEME2000 = circ.getPVCoordinates();
330 
331         //  Pointing direction
332         // ********************
333         // Get satellite attitude rotation, i.e rotation from EME2000 frame to satellite frame
334         FieldRotation<T> rotSatEME2000 = earthCenterAttitudeLaw.getAttitude(circ, date, circ.getFrame()).getRotation();
335 
336         //checked this values with the bodycenterpointing test values
337         // Transform Z axis from satellite frame to EME2000
338         FieldVector3D<T> zSatEME2000 = rotSatEME2000.applyInverseTo(Vector3D.PLUS_K);
339 
340         // Transform Z axis from EME2000 to ITRF2008
341         FieldVector3D<T> zSatITRF2008C = eme2000ToItrf.transformVector(zSatEME2000);
342 
343         // Transform satellite position/velocity from EME2000 to ITRF2008
344         FieldPVCoordinates<T> pvSatITRF2008C = eme2000ToItrf.transformPVCoordinates(pvSatEME2000);
345 
346 
347        // Line containing satellite point and following pointing direction
348         Line pointingLine = new Line(pvSatITRF2008C.getPosition().toVector3D(),
349                                      pvSatITRF2008C.getPosition().toVector3D().add(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
350                                                                       zSatITRF2008C.toVector3D()),
351                                      2.0e-8);
352 
353         // Check that the line contains Earth center
354         Assert.assertTrue(pointingLine.contains(Vector3D.ZERO));
355 
356     }
357 
358     private <T extends CalculusFieldElement<T>> void doTestQDot(final Field<T> field) {
359 
360         final double ae  = 6.378137e6;
361         final double c20 = -1.08263e-3;
362         final double c30 = 2.54e-6;
363         final double c40 = 1.62e-6;
364         final double c50 = 2.3e-7;
365         final double c60 = -5.5e-7;
366 
367         // Satellite position as circular parameters
368         T zero = field.getZero();
369         final T a     = zero.add(7178000.0);
370         final T e     = zero.add(7e-5);
371         final T i     = zero.add(FastMath.toRadians(50.));
372         final T pa    = zero.add(FastMath.toRadians(45.));
373         final T raan  = zero.add(FastMath.toRadians(270.));
374         final T m     = zero.add(FastMath.toRadians(5.3-270));
375         final T ehMu  = zero.add(3.9860047e14);
376 
377      // Computation date
378         FieldAbsoluteDate<T> date_comp= new FieldAbsoluteDate<>(field, new DateComponents(2008, 04, 07),
379                                                                 TimeComponents.H00,
380                                                                 TimeScalesFactory.getUTC());
381         // Orbit
382         FieldKeplerianOrbit<T> circ = new FieldKeplerianOrbit<>(a, e, i, pa, raan,
383                                                                 m, PositionAngle.MEAN,
384                                                                 FramesFactory.getEME2000(), date_comp, ehMu);
385         // WGS84 Earth model
386         OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
387                                                       Constants.WGS84_EARTH_FLATTENING,
388                                                       FramesFactory.getITRF(IERSConventions.IERS_2010, true));
389 
390         // Earth center pointing attitude provider
391         BodyCenterPointing earthCenterAttitudeLaw= new BodyCenterPointing(circ.getFrame(), earth);
392 
393 
394         final FieldAbsoluteDate<T> date = FieldAbsoluteDate.getJ2000Epoch(field).shiftedBy(584.);
395         final FieldVector3D<T> position = new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(6449822.));
396         final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
397         final FieldCircularOrbit<T> initialOrbit = new FieldCircularOrbit<>(new FieldPVCoordinates<>(position, velocity),
398                                                              FramesFactory.getEME2000(), date, ehMu);
399 
400         FieldEcksteinHechlerPropagator<T> propagator =
401                 new FieldEcksteinHechlerPropagator<>(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
402         propagator.setAttitudeProvider(earthCenterAttitudeLaw);
403 
404         List<WeightedObservedPoint> w0 = new ArrayList<WeightedObservedPoint>();
405         List<WeightedObservedPoint> w1 = new ArrayList<WeightedObservedPoint>();
406         List<WeightedObservedPoint> w2 = new ArrayList<WeightedObservedPoint>();
407         List<WeightedObservedPoint> w3 = new ArrayList<WeightedObservedPoint>();
408         for (double dt = -1; dt < 1; dt += 0.01) {
409             FieldRotation<T> rP = propagator.propagate(date.shiftedBy(dt)).getAttitude().getRotation();
410             w0.add(new WeightedObservedPoint(1, dt, rP.getQ0().getReal()));
411             w1.add(new WeightedObservedPoint(1, dt, rP.getQ1().getReal()));
412             w2.add(new WeightedObservedPoint(1, dt, rP.getQ2().getReal()));
413             w3.add(new WeightedObservedPoint(1, dt, rP.getQ3().getReal()));
414         }
415 
416         double q0DotRef = PolynomialCurveFitter.create(2).fit(w0)[1];
417         double q1DotRef = PolynomialCurveFitter.create(2).fit(w1)[1];
418         double q2DotRef = PolynomialCurveFitter.create(2).fit(w2)[1];
419         double q3DotRef = PolynomialCurveFitter.create(2).fit(w3)[1];
420 
421         FieldAttitude<T> a0 = propagator.propagate(date).getAttitude();
422         T   q0 = a0.getRotation().getQ0();
423         T   q1 = a0.getRotation().getQ1();
424         T   q2 = a0.getRotation().getQ2();
425         T   q3 = a0.getRotation().getQ3();
426         T   oX = a0.getSpin().getX();
427         T   oY = a0.getSpin().getY();
428         T   oZ = a0.getSpin().getZ();
429 
430         // first time-derivatives of the quaternion
431         double q0Dot = 0.5 * MathArrays.linearCombination(-q1.getReal(), oX.getReal(), -q2.getReal(), oY.getReal(), -q3.getReal(), oZ.getReal());
432         double q1Dot = 0.5 * MathArrays.linearCombination( q0.getReal(), oX.getReal(), -q3.getReal(), oY.getReal(),  q2.getReal(), oZ.getReal());
433         double q2Dot = 0.5 * MathArrays.linearCombination( q3.getReal(), oX.getReal(),  q0.getReal(), oY.getReal(), -q1.getReal(), oZ.getReal());
434         double q3Dot = 0.5 * MathArrays.linearCombination(-q2.getReal(), oX.getReal(),  q1.getReal(), oY.getReal(),  q0.getReal(), oZ.getReal());
435 
436         Assert.assertEquals(q0DotRef, q0Dot, 5.0e-9);
437         Assert.assertEquals(q1DotRef, q1Dot, 5.0e-9);
438         Assert.assertEquals(q2DotRef, q2Dot, 5.0e-9);
439         Assert.assertEquals(q3DotRef, q3Dot, 5.0e-9);
440 
441     }
442 
443     private <T extends CalculusFieldElement<T>> void doTestSpin(final Field<T> field) {
444 
445         final double ae  = 6.378137e6;
446         final double c20 = -1.08263e-3;
447         final double c30 = 2.54e-6;
448         final double c40 = 1.62e-6;
449         final double c50 = 2.3e-7;
450         final double c60 = -5.5e-7;
451 
452         // Satellite position as circular parameters
453         final T zero  = field.getZero();
454         final T a     = zero.add(7178000.0);
455         final T e     = zero.add(7e-5);
456         final T i     = zero.add(FastMath.toRadians(50.));
457         final T pa    = zero.add(FastMath.toRadians(45.));
458         final T raan  = zero.add(FastMath.toRadians(270.));
459         final T m     = zero.add(FastMath.toRadians(5.3-270));
460         final T ehMu  = zero.add(3.9860047e14);
461 
462         // Computation date
463         FieldAbsoluteDate<T> date_R = new FieldAbsoluteDate<>(field, new DateComponents(2008, 04, 07),
464                                                               TimeComponents.H00,
465                                                               TimeScalesFactory.getUTC());
466         // Orbit
467         FieldKeplerianOrbit<T> circ = new FieldKeplerianOrbit<>(a, e, i, pa, raan,
468                                                                 m, PositionAngle.MEAN,
469                                                                 FramesFactory.getEME2000(), date_R, ehMu);
470         // WGS84 Earth model
471         OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
472                                                       Constants.WGS84_EARTH_FLATTENING,
473                                                       FramesFactory.getITRF(IERSConventions.IERS_2010, true));
474 
475        // Earth center pointing attitude provider
476         BodyCenterPointing earthCenterAttitudeLaw= new BodyCenterPointing(circ.getFrame(), earth);
477 
478 
479         final FieldAbsoluteDate<T> date = FieldAbsoluteDate.getJ2000Epoch(field).shiftedBy(584.);
480         final FieldVector3D<T> position = new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(6449822.));
481         final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
482         final FieldCircularOrbit<T> initialOrbit = new FieldCircularOrbit<>(new FieldPVCoordinates<>(position, velocity),
483                                                              FramesFactory.getEME2000(), date, ehMu);
484 
485         FieldEcksteinHechlerPropagator<T> propagator =
486                 new FieldEcksteinHechlerPropagator<>(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
487         propagator.setAttitudeProvider(earthCenterAttitudeLaw);
488 
489         double h = 0.01;
490 
491         FieldSpacecraftState<T> s0     = propagator.propagate(date);
492         FieldSpacecraftState<T> sMinus = propagator.propagate(date.shiftedBy(-h));
493         FieldSpacecraftState<T> sPlus  = propagator.propagate(date.shiftedBy(h));
494 
495         // check spin is consistent with attitude evolution
496         T errorAngleMinus     = FieldRotation.distance(sMinus.shiftedBy(zero.add(h)).getAttitude().getRotation(),
497                                                        s0.getAttitude().getRotation());
498         T evolutionAngleMinus = FieldRotation.distance(sMinus.getAttitude().getRotation(),
499                                                        s0.getAttitude().getRotation());
500         Assert.assertEquals(0.0, errorAngleMinus.getReal(), 1.0e-6 * evolutionAngleMinus.getReal());
501         T errorAnglePlus      = FieldRotation.distance(s0.getAttitude().getRotation(),
502                                                        sPlus.shiftedBy(zero.add(-h)).getAttitude().getRotation());
503         T evolutionAnglePlus  = FieldRotation.distance(s0.getAttitude().getRotation(),
504                                                        sPlus.getAttitude().getRotation());
505         Assert.assertEquals(0.0, errorAnglePlus.getReal(), 1.0e-6 * evolutionAnglePlus.getReal());
506 
507         FieldVector3D<T> spin0 = s0.getAttitude().getSpin();
508         FieldVector3D<T> reference = FieldAngularCoordinates.estimateRate(sMinus.getAttitude().getRotation(),
509                                                              sPlus.getAttitude().getRotation(),
510                                                              2 * h);
511         Assert.assertTrue(spin0.getNorm().getReal() > 1.0e-3);
512         Assert.assertEquals(0.0, spin0.subtract(reference).getNorm().getReal(), 1.0e-13);
513 
514     }
515 
516     @Before
517     public void setUp() {
518         try {
519 
520             Utils.setDataRoot("regular-data");
521 
522             // Computation date
523             date = new AbsoluteDate(new DateComponents(2008, 04, 07),
524                                     TimeComponents.H00,
525                                     TimeScalesFactory.getUTC());
526 
527             // Satellite position as circular parameters
528             final double mu = 3.9860047e14;
529             final double raan = 270.;
530             circ =
531                 new CircularOrbit(7178000.0, 0.5e-4, -0.5e-4, FastMath.toRadians(50.), FastMath.toRadians(raan),
532                                        FastMath.toRadians(5.300 - raan), PositionAngle.MEAN,
533                                        FramesFactory.getEME2000(), date, mu);
534 
535 
536             // WGS84 Earth model
537             earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
538                                          Constants.WGS84_EARTH_FLATTENING,
539                                          FramesFactory.getITRF(IERSConventions.IERS_2010, true));
540 
541             // Transform from EME2000 to ITRF2008
542             eme2000ToItrf = FramesFactory.getEME2000().getTransformTo(earth.getBodyFrame(), date);
543 
544             // Create earth center pointing attitude provider
545             earthCenterAttitudeLaw = new BodyCenterPointing(circ.getFrame(), earth);
546 
547         } catch (OrekitException oe) {
548             Assert.fail(oe.getMessage());
549         }
550 
551     }
552 
553     @After
554     public void tearDown() {
555         date = null;
556         earth = null;
557         eme2000ToItrf = null;
558         earthCenterAttitudeLaw = null;
559         circ = null;
560     }
561 
562 }
563