1   /* Copyright 2002-2020 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.RealFieldElement;
25  import org.hipparchus.geometry.euclidean.threed.Rotation;
26  import org.hipparchus.geometry.euclidean.threed.Vector3D;
27  import org.hipparchus.util.Decimal64Field;
28  import org.hipparchus.util.FastMath;
29  import org.junit.After;
30  import org.junit.Assert;
31  import org.junit.Before;
32  import org.junit.Test;
33  import org.orekit.Utils;
34  import org.orekit.bodies.OneAxisEllipsoid;
35  import org.orekit.errors.OrekitException;
36  import org.orekit.frames.Frame;
37  import org.orekit.frames.FramesFactory;
38  import org.orekit.orbits.CircularOrbit;
39  import org.orekit.orbits.FieldOrbit;
40  import org.orekit.orbits.Orbit;
41  import org.orekit.orbits.PositionAngle;
42  import org.orekit.propagation.FieldSpacecraftState;
43  import org.orekit.propagation.Propagator;
44  import org.orekit.propagation.SpacecraftState;
45  import org.orekit.propagation.analytical.KeplerianPropagator;
46  import org.orekit.propagation.sampling.OrekitFixedStepHandler;
47  import org.orekit.time.AbsoluteDate;
48  import org.orekit.time.DateComponents;
49  import org.orekit.time.FieldAbsoluteDate;
50  import org.orekit.time.TimeComponents;
51  import org.orekit.time.TimeScalesFactory;
52  import org.orekit.utils.AngularDerivativesFilter;
53  import org.orekit.utils.IERSConventions;
54  import org.orekit.utils.TimeStampedAngularCoordinates;
55  
56  public class TabulatedProviderTest {
57  
58      // Computation date
59      private AbsoluteDate date;
60  
61      // Reference frame = ITRF
62      private Frame itrf;
63  
64      // Satellite position
65      CircularOrbit circOrbit;
66  
67      // Earth shape
68      OneAxisEllipsoid earthShape;
69  
70      @Test
71      public void testWithoutRate() {
72          double             samplingRate      = 10.0;
73          double             checkingRate      = 1.0;
74          int                n                 = 8;
75          AttitudeProvider   referenceProvider = new NadirPointing(circOrbit.getFrame(), earthShape);
76          List<TimeStampedAngularCoordinates> sample = createSample(samplingRate, referenceProvider);
77          final double       margin            = samplingRate * n / 2;
78          final AbsoluteDate start             = sample.get(0).getDate().shiftedBy(margin);
79          final AbsoluteDate end               = sample.get(sample.size() - 1).getDate().shiftedBy(-margin);
80          TabulatedProvider  provider          = new TabulatedProvider(circOrbit.getFrame(), sample, n,
81                                                                       AngularDerivativesFilter.USE_R);
82          Assert.assertEquals(0.0, checkError(start, end, checkingRate, referenceProvider, provider), 2.2e-14);
83      }
84  
85      @Test
86      public void testWithRate() {
87          double             samplingRate      = 10.0;
88          double             checkingRate      = 1.0;
89          int                n                 = 8;
90          AttitudeProvider   referenceProvider = new NadirPointing(circOrbit.getFrame(), earthShape);
91          List<TimeStampedAngularCoordinates> sample = createSample(samplingRate, referenceProvider);
92          final double       margin            = samplingRate * n / 2;
93          final AbsoluteDate start             = sample.get(0).getDate().shiftedBy(margin);
94          final AbsoluteDate end               = sample.get(sample.size() - 1).getDate().shiftedBy(-margin);
95          TabulatedProvider  provider          = new TabulatedProvider(circOrbit.getFrame(), sample, n,
96                                                                       AngularDerivativesFilter.USE_RR);
97          Assert.assertEquals(0.0, checkError(start, end, checkingRate, referenceProvider, provider), 1.3e-11);
98      }
99  
100     @Test
101     public void testWithAcceleration() {
102         double             samplingRate      = 10.0;
103         double             checkingRate      = 1.0;
104         int                n                 = 8;
105         AttitudeProvider   referenceProvider = new NadirPointing(circOrbit.getFrame(), earthShape);
106         List<TimeStampedAngularCoordinates> sample = createSample(samplingRate, referenceProvider);
107         final double       margin            = samplingRate * n / 2;
108         final AbsoluteDate start             = sample.get(0).getDate().shiftedBy(margin);
109         final AbsoluteDate end               = sample.get(sample.size() - 1).getDate().shiftedBy(-margin);
110         TabulatedProvider  provider          = new TabulatedProvider(circOrbit.getFrame(), sample, n,
111                                                                      AngularDerivativesFilter.USE_RRA);
112         Assert.assertEquals(0.0, checkError(start, end, checkingRate, referenceProvider, provider), 4.3e-9);
113         checkField(Decimal64Field.getInstance(), provider, circOrbit, circOrbit.getDate(), circOrbit.getFrame());
114     }
115 
116     private List<TimeStampedAngularCoordinates> createSample(double samplingRate, AttitudeProvider referenceProvider) {
117 
118         // reference propagator, using a yaw compensation law
119         final KeplerianPropagator referencePropagator = new KeplerianPropagator(circOrbit);
120         referencePropagator.setAttitudeProvider(referenceProvider);
121 
122         // create sample
123         final List<TimeStampedAngularCoordinates> sample = new ArrayList<TimeStampedAngularCoordinates>();
124         referencePropagator.setMasterMode(samplingRate, new OrekitFixedStepHandler() {
125 
126             public void handleStep(SpacecraftState currentState, boolean isLast) {
127                 sample.add(currentState.getAttitude().getOrientation());
128             }
129 
130         });
131         referencePropagator.propagate(circOrbit.getDate().shiftedBy(2 * circOrbit.getKeplerianPeriod()));
132 
133         return sample;
134 
135     }
136 
137     private double checkError(final AbsoluteDate start, AbsoluteDate end, double checkingRate,
138                               final AttitudeProvider referenceProvider, TabulatedProvider provider) {
139 
140         // prepare an interpolating provider, using only internal steps
141         // (i.e. ignoring interpolation near boundaries)
142         Propagator interpolatingPropagator = new KeplerianPropagator(circOrbit.shiftedBy(start.durationFrom(circOrbit.getDate())));
143         interpolatingPropagator.setAttitudeProvider(provider);
144 
145         // compute interpolation error on the internal steps .
146         final double[] error = new double[1];
147         interpolatingPropagator.setMasterMode(checkingRate, new OrekitFixedStepHandler() {
148 
149             public void init(SpacecraftState s0, AbsoluteDate t, double step) {
150                 error[0] = 0.0;
151             }
152 
153             public void handleStep(SpacecraftState currentState, boolean isLast) {
154                 Attitude interpolated = currentState.getAttitude();
155                 Attitude reference    = referenceProvider.getAttitude(currentState.getOrbit(),
156                                                                       currentState.getDate(),
157                                                                       currentState.getFrame());
158                 double localError = Rotation.distance(interpolated.getRotation(), reference.getRotation());
159                 error[0] = FastMath.max(error[0], localError);
160             }
161 
162         });
163 
164         interpolatingPropagator.propagate(end);
165 
166         return error[0];
167 
168     }
169 
170     private <T extends RealFieldElement<T>> void checkField(final Field<T> field, final AttitudeProvider provider,
171                                                             final Orbit orbit, final AbsoluteDate date,
172                                                             final Frame frame) {
173         Attitude attitudeD = provider.getAttitude(orbit, date, frame);
174         final FieldOrbit<T> orbitF = new FieldSpacecraftState<>(field, new SpacecraftState(orbit)).getOrbit();
175         final FieldAbsoluteDate<T> dateF = new FieldAbsoluteDate<>(field, date);
176         FieldAttitude<T> attitudeF = provider.getAttitude(orbitF, dateF, frame);
177         Assert.assertEquals(0.0, Rotation.distance(attitudeD.getRotation(), attitudeF.getRotation().toRotation()), 1.0e-15);
178         Assert.assertEquals(0.0, Vector3D.distance(attitudeD.getSpin(), attitudeF.getSpin().toVector3D()), 1.0e-15);
179         Assert.assertEquals(0.0, Vector3D.distance(attitudeD.getRotationAcceleration(), attitudeF.getRotationAcceleration().toVector3D()), 1.0e-15);
180     }
181 
182     @Before
183     public void setUp() {
184         try {
185             Utils.setDataRoot("regular-data");
186 
187             // Computation date
188             date = new AbsoluteDate(new DateComponents(2008, 04, 07),
189                                     TimeComponents.H00,
190                                     TimeScalesFactory.getUTC());
191 
192             // Body mu
193             final double mu = 3.9860047e14;
194 
195             // Reference frame = ITRF
196             itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
197 
198             //  Satellite position
199             circOrbit =
200                 new CircularOrbit(7178000.0, 0.5e-4, -0.5e-4, FastMath.toRadians(50.), FastMath.toRadians(270.),
201                                        FastMath.toRadians(5.300), PositionAngle.MEAN,
202                                        FramesFactory.getEME2000(), date, mu);
203 
204             // Elliptic earth shape
205             earthShape =
206                 new OneAxisEllipsoid(6378136.460, 1 / 298.257222101, itrf);
207 
208         } catch (OrekitException oe) {
209             Assert.fail(oe.getMessage());
210         }
211 
212     }
213 
214     @After
215     public void tearDown() {
216         date = null;
217         itrf = null;
218         circOrbit = null;
219         earthShape = null;
220     }
221 
222 }
223