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.geometry.euclidean.threed.Rotation;
24  import org.hipparchus.geometry.euclidean.threed.Vector3D;
25  import org.hipparchus.util.FastMath;
26  import org.junit.Assert;
27  import org.junit.Test;
28  import org.orekit.Utils;
29  import org.orekit.bodies.OneAxisEllipsoid;
30  import org.orekit.frames.FramesFactory;
31  import org.orekit.orbits.CircularOrbit;
32  import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
33  import org.orekit.time.AbsoluteDate;
34  import org.orekit.utils.AngularCoordinates;
35  import org.orekit.utils.Constants;
36  import org.orekit.utils.IERSConventions;
37  import org.orekit.utils.PVCoordinates;
38  
39  public class AttitudeTest {
40  
41      @Test
42      public void testZeroRate() {
43          Attitude attitude = new Attitude(AbsoluteDate.J2000_EPOCH, FramesFactory.getEME2000(),
44                                           new Rotation(0.48, 0.64, 0.36, 0.48, false),
45                                           Vector3D.ZERO, Vector3D.ZERO);
46          Assert.assertEquals(Vector3D.ZERO, attitude.getSpin());
47          double dt = 10.0;
48          Attitude shifted = attitude.shiftedBy(dt);
49          Assert.assertEquals(Vector3D.ZERO, shifted.getRotationAcceleration());
50          Assert.assertEquals(Vector3D.ZERO, shifted.getSpin());
51          Assert.assertEquals(0.0, Rotation.distance(attitude.getRotation(), shifted.getRotation()), 1.0e-15);
52      }
53  
54      @Test
55      public void testShift() {
56          double rate = 2 * FastMath.PI / (12 * 60);
57          Attitude attitude = new Attitude(AbsoluteDate.J2000_EPOCH, FramesFactory.getEME2000(),
58                                           Rotation.IDENTITY,
59                                           new Vector3D(rate, Vector3D.PLUS_K), Vector3D.ZERO);
60          Assert.assertEquals(rate, attitude.getSpin().getNorm(), 1.0e-10);
61          double dt = 10.0;
62          double alpha = rate * dt;
63          Attitude shifted = attitude.shiftedBy(dt);
64          Assert.assertEquals(rate, shifted.getSpin().getNorm(), 1.0e-10);
65          Assert.assertEquals(alpha, Rotation.distance(attitude.getRotation(), shifted.getRotation()), 1.0e-10);
66  
67          Vector3D xSat = shifted.getRotation().applyInverseTo(Vector3D.PLUS_I);
68          Assert.assertEquals(0.0, xSat.subtract(new Vector3D(FastMath.cos(alpha), FastMath.sin(alpha), 0)).getNorm(), 1.0e-10);
69          Vector3D ySat = shifted.getRotation().applyInverseTo(Vector3D.PLUS_J);
70          Assert.assertEquals(0.0, ySat.subtract(new Vector3D(-FastMath.sin(alpha), FastMath.cos(alpha), 0)).getNorm(), 1.0e-10);
71          Vector3D zSat = shifted.getRotation().applyInverseTo(Vector3D.PLUS_K);
72          Assert.assertEquals(0.0, zSat.subtract(Vector3D.PLUS_K).getNorm(), 1.0e-10);
73  
74      }
75  
76      @Test
77      public void testSpin() {
78          double rate = 2 * FastMath.PI / (12 * 60);
79          Attitude attitude = new Attitude(AbsoluteDate.J2000_EPOCH, FramesFactory.getEME2000(),
80                                           new Rotation(0.48, 0.64, 0.36, 0.48, false),
81                                           new Vector3D(rate, Vector3D.PLUS_K), Vector3D.ZERO);
82          Assert.assertEquals(rate, attitude.getSpin().getNorm(), 1.0e-10);
83          double dt = 10.0;
84          Attitude shifted = attitude.shiftedBy(dt);
85          Assert.assertEquals(rate, shifted.getSpin().getNorm(), 1.0e-10);
86          Assert.assertEquals(rate * dt, Rotation.distance(attitude.getRotation(), shifted.getRotation()), 1.0e-10);
87  
88          Vector3D shiftedX  = shifted.getRotation().applyInverseTo(Vector3D.PLUS_I);
89          Vector3D shiftedY  = shifted.getRotation().applyInverseTo(Vector3D.PLUS_J);
90          Vector3D shiftedZ  = shifted.getRotation().applyInverseTo(Vector3D.PLUS_K);
91          Vector3D originalX = attitude.getRotation().applyInverseTo(Vector3D.PLUS_I);
92          Vector3D originalY = attitude.getRotation().applyInverseTo(Vector3D.PLUS_J);
93          Vector3D originalZ = attitude.getRotation().applyInverseTo(Vector3D.PLUS_K);
94          Assert.assertEquals( FastMath.cos(rate * dt), Vector3D.dotProduct(shiftedX, originalX), 1.0e-10);
95          Assert.assertEquals( FastMath.sin(rate * dt), Vector3D.dotProduct(shiftedX, originalY), 1.0e-10);
96          Assert.assertEquals( 0.0,                 Vector3D.dotProduct(shiftedX, originalZ), 1.0e-10);
97          Assert.assertEquals(-FastMath.sin(rate * dt), Vector3D.dotProduct(shiftedY, originalX), 1.0e-10);
98          Assert.assertEquals( FastMath.cos(rate * dt), Vector3D.dotProduct(shiftedY, originalY), 1.0e-10);
99          Assert.assertEquals( 0.0,                 Vector3D.dotProduct(shiftedY, originalZ), 1.0e-10);
100         Assert.assertEquals( 0.0,                 Vector3D.dotProduct(shiftedZ, originalX), 1.0e-10);
101         Assert.assertEquals( 0.0,                 Vector3D.dotProduct(shiftedZ, originalY), 1.0e-10);
102         Assert.assertEquals( 1.0,                 Vector3D.dotProduct(shiftedZ, originalZ), 1.0e-10);
103 
104         Vector3D forward = AngularCoordinates.estimateRate(attitude.getRotation(), shifted.getRotation(), dt);
105         Assert.assertEquals(0.0, forward.subtract(attitude.getSpin()).getNorm(), 1.0e-10);
106 
107         Vector3D reversed = AngularCoordinates.estimateRate(shifted.getRotation(), attitude.getRotation(), dt);
108         Assert.assertEquals(0.0, reversed.add(attitude.getSpin()).getNorm(), 1.0e-10);
109 
110     }
111 
112     @Test
113     public void testInterpolation() {
114 
115         Utils.setDataRoot("regular-data");
116         final double ehMu  = 3.9860047e14;
117         final double ae  = 6.378137e6;
118         final double c20 = -1.08263e-3;
119         final double c30 = 2.54e-6;
120         final double c40 = 1.62e-6;
121         final double c50 = 2.3e-7;
122         final double c60 = -5.5e-7;
123 
124         final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
125         final Vector3D position = new Vector3D(3220103., 69623., 6449822.);
126         final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
127         final CircularOrbit initialOrbit = new CircularOrbit(new PVCoordinates(position, velocity),
128                                                              FramesFactory.getEME2000(), date, ehMu);
129 
130         EcksteinHechlerPropagator propagator =
131                 new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
132         OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
133                                                       Constants.WGS84_EARTH_FLATTENING,
134                                                       FramesFactory.getITRF(IERSConventions.IERS_2010, true));
135         propagator.setAttitudeProvider(new BodyCenterPointing(initialOrbit.getFrame(), earth));
136         final Attitude initialAttitude = propagator.propagate(initialOrbit.getDate()).getAttitude();
137 
138         // set up a 5 points sample
139         List<Attitude> sample = new ArrayList<Attitude>();
140         for (double dt = 0; dt < 251.0; dt += 60.0) {
141             sample.add(propagator.propagate(date.shiftedBy(dt)).getAttitude());
142         }
143 
144         // well inside the sample, interpolation should be better than quadratic shift
145         double maxShiftAngleError = 0;
146         double maxInterpolationAngleError = 0;
147         double maxShiftRateError = 0;
148         double maxInterpolationRateError = 0;
149         for (double dt = 0; dt < 240.0; dt += 1.0) {
150             AbsoluteDate t                 = initialOrbit.getDate().shiftedBy(dt);
151             Attitude propagated            = propagator.propagate(t).getAttitude();
152             double shiftAngleError         = Rotation.distance(propagated.getRotation(),
153                                                                initialAttitude.shiftedBy(dt).getRotation());
154             double interpolationAngleError = Rotation.distance(propagated.getRotation(),
155                                                                initialAttitude.interpolate(t, sample).getRotation());
156             double shiftRateError          = Vector3D.distance(propagated.getSpin(),
157                                                                initialAttitude.shiftedBy(dt).getSpin());
158             double interpolationRateError  = Vector3D.distance(propagated.getSpin(),
159                                                                initialAttitude.interpolate(t, sample).getSpin());
160             maxShiftAngleError             = FastMath.max(maxShiftAngleError, shiftAngleError);
161             maxInterpolationAngleError     = FastMath.max(maxInterpolationAngleError, interpolationAngleError);
162             maxShiftRateError              = FastMath.max(maxShiftRateError, shiftRateError);
163             maxInterpolationRateError      = FastMath.max(maxInterpolationRateError, interpolationRateError);
164         }
165         Assert.assertTrue(maxShiftAngleError         > 4.0e-6);
166         Assert.assertTrue(maxInterpolationAngleError < 1.5e-13);
167         Assert.assertTrue(maxShiftRateError          > 6.0e-8);
168         Assert.assertTrue(maxInterpolationRateError  < 2.5e-14);
169 
170         // past sample end, interpolation error should increase, but still be far better than quadratic shift
171         maxShiftAngleError = 0;
172         maxInterpolationAngleError = 0;
173         maxShiftRateError = 0;
174         maxInterpolationRateError = 0;
175         for (double dt = 250.0; dt < 300.0; dt += 1.0) {
176             AbsoluteDate t                 = initialOrbit.getDate().shiftedBy(dt);
177             Attitude propagated            = propagator.propagate(t).getAttitude();
178             double shiftAngleError         = Rotation.distance(propagated.getRotation(),
179                                                                initialAttitude.shiftedBy(dt).getRotation());
180             double interpolationAngleError = Rotation.distance(propagated.getRotation(),
181                                                                initialAttitude.interpolate(t, sample).getRotation());
182             double shiftRateError          = Vector3D.distance(propagated.getSpin(),
183                                                                initialAttitude.shiftedBy(dt).getSpin());
184             double interpolationRateError  = Vector3D.distance(propagated.getSpin(),
185                                                                initialAttitude.interpolate(t, sample).getSpin());
186             maxShiftAngleError             = FastMath.max(maxShiftAngleError, shiftAngleError);
187             maxInterpolationAngleError     = FastMath.max(maxInterpolationAngleError, interpolationAngleError);
188             maxShiftRateError              = FastMath.max(maxShiftRateError, shiftRateError);
189             maxInterpolationRateError      = FastMath.max(maxInterpolationRateError, interpolationRateError);
190         }
191         Assert.assertTrue(maxShiftAngleError         > 9.0e-6);
192         Assert.assertTrue(maxInterpolationAngleError < 6.0e-11);
193         Assert.assertTrue(maxShiftRateError          > 9.0e-8);
194         Assert.assertTrue(maxInterpolationRateError  < 4.0e-12);
195 
196     }
197 
198 }
199