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.propagation.analytical.tle;
18  
19  import java.io.FileNotFoundException;
20  import java.io.IOException;
21  import java.io.UnsupportedEncodingException;
22  import java.text.ParseException;
23  
24  import org.hipparchus.linear.RealMatrix;
25  import org.hipparchus.util.FastMath;
26  import org.junit.Assert;
27  import org.junit.Before;
28  import org.junit.Test;
29  import org.orekit.Utils;
30  import org.orekit.attitudes.Attitude;
31  import org.orekit.errors.OrekitException;
32  import org.orekit.errors.OrekitMessages;
33  import org.orekit.frames.Frame;
34  import org.orekit.orbits.CartesianOrbit;
35  import org.orekit.orbits.Orbit;
36  import org.orekit.orbits.OrbitType;
37  import org.orekit.orbits.PositionAngle;
38  import org.orekit.propagation.SpacecraftState;
39  import org.orekit.propagation.numerical.NumericalPropagator;
40  import org.orekit.time.AbsoluteDate;
41  import org.orekit.utils.ParameterDriver;
42  import org.orekit.utils.ParameterDriversList;
43  
44  @Deprecated
45  public class TLEPartialDerivativesTest {
46      
47      
48      // build two TLEs in order to test SGP4 and SDP4 algorithms
49      private TLE tleGPS;
50      
51      private TLE tleSPOT;
52  
53      @Test
54      public void testStateJacobianSDP4() throws FileNotFoundException, UnsupportedEncodingException {
55          doTestStateJacobian(2.53e-9, tleGPS);
56      }
57  
58      @Test
59      public void testStateJacobianSGP4() throws FileNotFoundException, UnsupportedEncodingException {
60          doTestStateJacobian(7.65e-10, tleSPOT);
61      }
62  
63      @Test
64      public void testBStarDerivatives() throws ParseException, IOException {
65          doTestParametersDerivatives(TLE.B_STAR,
66                                      5.16e-3,
67                                      tleSPOT);
68      }
69     
70      @Test(expected=OrekitException.class)
71      public void testNotInitialized() {
72          String line1 = "1 37753U 11036A   12090.13205652 -.00000006  00000-0  00000+0 0  2272";
73          String line2 = "2 37753  55.0032 176.5796 0004733  13.2285 346.8266  2.00565440  5153";       
74          TLE tle = new TLE(line1, line2);
75  
76          TLEPropagator propagator = TLEPropagator.selectExtrapolator(tle);
77          new TLEPartialDerivativesEquations("partials", propagator).getMapper();
78       }
79      
80      @Test(expected=OrekitException.class)
81      public void testTooSmallDimension() {
82          String line1 = "1 37753U 11036A   12090.13205652 -.00000006  00000-0  00000+0 0  2272";
83          String line2 = "2 37753  55.0032 176.5796 0004733  13.2285 346.8266  2.00565440  5153";       
84          TLE tle = new TLE(line1, line2);
85  
86          TLEPropagator propagator = TLEPropagator.selectExtrapolator(tle);
87          TLEPartialDerivativesEquations partials = new TLEPartialDerivativesEquations("partials", propagator);
88          partials.setInitialJacobians(propagator.getInitialState(),
89                                       new double[5][6], new double[6][2]);
90       }
91  
92      @Test(expected=OrekitException.class)
93      public void testTooLargeDimension() {
94          String line1 = "1 37753U 11036A   12090.13205652 -.00000006  00000-0  00000+0 0  2272";
95          String line2 = "2 37753  55.0032 176.5796 0004733  13.2285 346.8266  2.00565440  5153";       
96          TLE tle = new TLE(line1, line2);
97  
98          TLEPropagator propagator = TLEPropagator.selectExtrapolator(tle);
99          TLEPartialDerivativesEquations partials = new TLEPartialDerivativesEquations("partials", propagator);
100         partials.setInitialJacobians(propagator.getInitialState(),
101                                      new double[8][6], new double[6][2]);
102      }
103     
104     @Test(expected=OrekitException.class)
105     public void testMismatchedDimensions() {
106         String line1 = "1 37753U 11036A   12090.13205652 -.00000006  00000-0  00000+0 0  2272";
107         String line2 = "2 37753  55.0032 176.5796 0004733  13.2285 346.8266  2.00565440  5153";       
108         TLE tle = new TLE(line1, line2);
109 
110         TLEPropagator propagator = TLEPropagator.selectExtrapolator(tle);
111         TLEPartialDerivativesEquations partials = new TLEPartialDerivativesEquations("partials", propagator);
112         partials.setInitialJacobians(propagator.getInitialState(),
113                                      new double[6][6], new double[7][2]);
114      }
115     
116     @Test
117     public void testWrongParametersDimension() {
118         String line1 = "1 37753U 11036A   12090.13205652 -.00000006  00000-0  00000+0 0  2272";
119         String line2 = "2 37753  55.0032 176.5796 0004733  13.2285 346.8266  2.00565440  5153";       
120         TLE tle = new TLE(line1, line2);
121 
122         TLEPropagator propagator = TLEPropagator.selectExtrapolator(tle);
123         TLEPartialDerivativesEquations partials = new TLEPartialDerivativesEquations("partials", propagator);
124         try {
125             partials.setInitialJacobians(propagator.getInitialState(),
126                                          new double[6][6], new double[6][3]);
127             Assert.fail("an exception should have been thrown");
128         } catch (OrekitException oe) {
129             Assert.assertEquals(OrekitMessages.INITIAL_MATRIX_AND_PARAMETERS_NUMBER_MISMATCH,
130                                 oe.getSpecifier());
131         }
132     }
133     
134     private void doTestStateJacobian(double tolerance, TLE tle)
135         throws FileNotFoundException, UnsupportedEncodingException {
136 
137         // compute state Jacobian using PartialDerivatives
138         TLEPropagator propagator = TLEPropagator.selectExtrapolator(tle);
139         final SpacecraftState initialState = propagator.getInitialState();
140         final AbsoluteDate target = initialState.getDate().shiftedBy(initialState.getKeplerianPeriod());
141         Orbit orbit = propagator.getInitialState().getOrbit();
142         TLEPartialDerivativesEquations partials = new TLEPartialDerivativesEquations("partials", propagator);
143         final SpacecraftState initState = partials.setInitialJacobians(propagator.getInitialState());
144         final double[] stateVector = new double[6];
145         OrbitType.CARTESIAN.mapOrbitToArray(initState.getOrbit(), PositionAngle.MEAN, stateVector, null);
146         final TLEJacobiansMapper mapper = partials.getMapper();
147         propagator.resetInitialState(initState);
148         final SpacecraftState endState = propagator.propagate(target);
149         mapper.setReferenceState(endState);
150         RealMatrix dYdY0 = mapper.getStateTransitionMatrix(endState);
151 
152         // compute reference state Jacobian using finite differences
153         double[][] dYdY0Ref = new double[6][6];
154         TLEPropagator propagator2;
155         double[] steps = NumericalPropagator.tolerances(10, orbit, OrbitType.CARTESIAN)[0];
156         for (int i = 0; i < 6; ++i) {
157             propagator2 = TLEPropagator.selectExtrapolator(TLE.stateToTLE(shiftState(initialState, OrbitType.CARTESIAN, -4 * steps[i], i), tle));
158             SpacecraftState sM4h = propagator2.propagate(target);
159             propagator2 = TLEPropagator.selectExtrapolator(TLE.stateToTLE(shiftState(initialState, OrbitType.CARTESIAN, -3 * steps[i], i), tle));
160             SpacecraftState sM3h = propagator2.propagate(target);
161             propagator2 = TLEPropagator.selectExtrapolator(TLE.stateToTLE(shiftState(initialState, OrbitType.CARTESIAN, -2 * steps[i], i), tle));
162             SpacecraftState sM2h = propagator2.propagate(target);
163             propagator2 = TLEPropagator.selectExtrapolator(TLE.stateToTLE(shiftState(initialState, OrbitType.CARTESIAN, -1 * steps[i], i), tle));
164             SpacecraftState sM1h = propagator2.propagate(target);
165             propagator2 = TLEPropagator.selectExtrapolator(TLE.stateToTLE(shiftState(initialState, OrbitType.CARTESIAN, +1 * steps[i], i), tle));
166             SpacecraftState sP1h = propagator2.propagate(target);
167             propagator2 = TLEPropagator.selectExtrapolator(TLE.stateToTLE(shiftState(initialState, OrbitType.CARTESIAN, +2 * steps[i], i), tle));
168             SpacecraftState sP2h = propagator2.propagate(target);
169             propagator2 = TLEPropagator.selectExtrapolator(TLE.stateToTLE(shiftState(initialState, OrbitType.CARTESIAN, +3 * steps[i], i), tle));
170             SpacecraftState sP3h = propagator2.propagate(target);
171             propagator2 = TLEPropagator.selectExtrapolator(TLE.stateToTLE(shiftState(initialState, OrbitType.CARTESIAN, +4 * steps[i], i), tle));
172             SpacecraftState sP4h = propagator2.propagate(target);
173             fillJacobianColumn(dYdY0Ref, i, OrbitType.CARTESIAN, steps[i],
174                                sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
175         }
176 
177         for (int i = 0; i < 6; ++i) {
178             for (int j = 0; j < 6; ++j) {
179                 if (stateVector[i] != 0) {
180                     double error = FastMath.abs((dYdY0.getEntry(i, j) - dYdY0Ref[i][j]) / stateVector[i]) * steps[j];
181                     Assert.assertEquals(0, error, tolerance);
182                 }
183             }
184         }
185     }
186 
187     private void doTestParametersDerivatives(String parameterName, double tolerance, TLE tle) {
188 
189         // compute state Jacobian using PartialDerivatives
190         ParameterDriversList bound = new ParameterDriversList();
191 
192         for (final ParameterDriver driver : tle.getParametersDrivers()) {
193             if (driver.getName().equals(parameterName)) {
194                 driver.setSelected(true);
195                 bound.add(driver);
196             } else {
197                 driver.setSelected(false);
198             }
199         }
200         TLEPropagator propagator = TLEPropagator.selectExtrapolator(tle);
201         final SpacecraftState initialState = propagator.getInitialState();
202         final AbsoluteDate target = initialState.getDate().shiftedBy(initialState.getKeplerianPeriod());
203         TLEPartialDerivativesEquations partials = new TLEPartialDerivativesEquations("partials", propagator);
204         final SpacecraftState endState = partials.setInitialJacobians(propagator.propagate(target));
205         final double[] stateVector = new double[6];
206         OrbitType.CARTESIAN.mapOrbitToArray(initialState.getOrbit(), PositionAngle.MEAN, stateVector, null);
207         final TLEJacobiansMapper mapper = partials.getMapper();
208         mapper.setReferenceState(endState);
209         RealMatrix dYdP = mapper.getParametersJacobian(endState);
210 
211         // compute reference Jacobian using finite differences
212         
213         OrbitType orbitType = OrbitType.CARTESIAN;
214         TLEPropagator propagator2 = TLEPropagator.selectExtrapolator(tle);
215         double[][] dYdPRef = new double[6][1];
216 
217         ParameterDriver selected = bound.getDrivers().get(0);
218         double p0 = selected.getReferenceValue();
219         double h  = selected.getScale();
220         selected.setValue(p0 - 4 * h);
221         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
222         SpacecraftState sM4h = propagator2.propagate(target);
223         selected.setValue(p0 - 3 * h);
224         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
225         SpacecraftState sM3h = propagator2.propagate(target);
226         selected.setValue(p0 - 2 * h);
227         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
228         SpacecraftState sM2h = propagator2.propagate(target);
229         selected.setValue(p0 - 1 * h);
230         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
231         SpacecraftState sM1h = propagator2.propagate(target);
232         selected.setValue(p0 + 1 * h);
233         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
234         SpacecraftState sP1h = propagator2.propagate(target);
235         selected.setValue(p0 + 2 * h);
236         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
237         SpacecraftState sP2h = propagator2.propagate(target);
238         selected.setValue(p0 + 3 * h);
239         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
240         SpacecraftState sP3h = propagator2.propagate(target);
241         selected.setValue(p0 + 4 * h);
242         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
243         SpacecraftState sP4h = propagator2.propagate(target);
244         fillJacobianColumn(dYdPRef, 0, orbitType, h,
245                            sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
246        
247         for (int i = 0; i < 6; ++i) {
248             Assert.assertEquals(dYdPRef[i][0], dYdP.getEntry(i, 0), FastMath.abs(tolerance));
249         }
250 
251     }
252 
253     private void fillJacobianColumn(double[][] jacobian, int column,
254                                     OrbitType orbitType, double h,
255                                     SpacecraftState sM4h, SpacecraftState sM3h,
256                                     SpacecraftState sM2h, SpacecraftState sM1h,
257                                     SpacecraftState sP1h, SpacecraftState sP2h,
258                                     SpacecraftState sP3h, SpacecraftState sP4h) {
259         double[] aM4h = stateToArray(sM4h, orbitType)[0];
260         double[] aM3h = stateToArray(sM3h, orbitType)[0];
261         double[] aM2h = stateToArray(sM2h, orbitType)[0];
262         double[] aM1h = stateToArray(sM1h, orbitType)[0];
263         double[] aP1h = stateToArray(sP1h, orbitType)[0];
264         double[] aP2h = stateToArray(sP2h, orbitType)[0];
265         double[] aP3h = stateToArray(sP3h, orbitType)[0];
266         double[] aP4h = stateToArray(sP4h, orbitType)[0];
267         for (int i = 0; i < jacobian.length; ++i) {
268             jacobian[i][column] = ( -3 * (aP4h[i] - aM4h[i]) +
269                                     32 * (aP3h[i] - aM3h[i]) -
270                                    168 * (aP2h[i] - aM2h[i]) +
271                                    672 * (aP1h[i] - aM1h[i])) / (840 * h);
272         }
273     }
274 
275     private SpacecraftState shiftState(SpacecraftState state, OrbitType orbitType,
276                                        double delta, int column) {
277 
278         double[][] array = stateToArray(state, orbitType);
279         array[0][column] += delta;
280 
281         return arrayToState(array, state.getFrame(), state.getDate(),
282                             state.getMu(), state.getAttitude());
283 
284     }
285     
286     
287     
288     private double[][] stateToArray(SpacecraftState state, OrbitType orbitType) {
289           double[][] array = new double[2][6];
290 
291           orbitType.mapOrbitToArray(state.getOrbit(), PositionAngle.MEAN, array[0], array[1]);
292           return array;
293       }
294 
295 
296     private SpacecraftState arrayToState(double[][] array, 
297                                            Frame frame, AbsoluteDate date, double mu,
298                                            Attitude attitude) {
299         CartesianOrbit orbit = (CartesianOrbit) OrbitType.CARTESIAN.mapArrayToOrbit(array[0], array[1], PositionAngle.MEAN, date, mu, frame);
300         return new SpacecraftState(orbit, attitude);
301     }
302 
303     private TLE newTLE(final TLE template, final double newBStar) {
304         return new TLE(template.getSatelliteNumber(), template.getClassification(),
305                        template.getLaunchYear(), template.getLaunchNumber(), template.getLaunchPiece(),
306                        template.getEphemerisType(), template.getElementNumber(), template.getDate(),
307                        template.getMeanMotion(), template.getMeanMotionFirstDerivative(), template.getMeanMotionSecondDerivative(),
308                        template.getE(), template.getI(), template.getPerigeeArgument(),
309                        template.getRaan(), template.getMeanAnomaly(), template.getRevolutionNumberAtEpoch(),
310                        newBStar);
311     }
312 
313     @Before
314     public void setUp() {
315         Utils.setDataRoot("regular-data");
316         
317         // GPS TLE propagation will use SDP4
318         String line1GPS = "1 11783U 80032A   03300.87313441  .00000062  00000-0  10000-3 0  6416";
319         String line2GPS = "2 11783  62.0472 164.2367 0320924  39.0039 323.3716  2.03455768173530"; 
320         tleGPS = new TLE(line1GPS, line2GPS);
321         
322         // SPOT TLE propagation will use SGP4
323         String line1SPOT = "1 22823U 93061A   03339.49496229  .00000173  00000-0  10336-3 0   133";
324         String line2SPOT = "2 22823  98.4132 359.2998 0017888 100.4310 259.8872 14.18403464527664";
325         tleSPOT = new TLE(line1SPOT, line2SPOT);
326 
327     }
328 
329 }