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.integration;
18  
19  import java.util.Collections;
20  
21  import org.hipparchus.geometry.euclidean.threed.Vector3D;
22  import org.hipparchus.linear.MatrixUtils;
23  import org.hipparchus.linear.RealMatrix;
24  import org.hipparchus.ode.DenseOutputModel;
25  import org.hipparchus.ode.ExpandableODE;
26  import org.hipparchus.ode.ODEState;
27  import org.hipparchus.ode.OrdinaryDifferentialEquation;
28  import org.hipparchus.ode.SecondaryODE;
29  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
30  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
31  import org.hipparchus.ode.nonstiff.EulerIntegrator;
32  import org.junit.Assert;
33  import org.junit.Before;
34  import org.junit.Test;
35  import org.orekit.Utils;
36  import org.orekit.attitudes.CelestialBodyPointed;
37  import org.orekit.attitudes.InertialProvider;
38  import org.orekit.bodies.CelestialBodyFactory;
39  import org.orekit.errors.OrekitInternalError;
40  import org.orekit.forces.gravity.potential.GravityFieldFactory;
41  import org.orekit.forces.gravity.potential.ICGEMFormatReader;
42  import org.orekit.frames.FramesFactory;
43  import org.orekit.orbits.EquinoctialOrbit;
44  import org.orekit.orbits.Orbit;
45  import org.orekit.orbits.OrbitType;
46  import org.orekit.orbits.PositionAngle;
47  import org.orekit.propagation.BoundedPropagator;
48  import org.orekit.propagation.EphemerisGenerator;
49  import org.orekit.propagation.MatricesHarvester;
50  import org.orekit.propagation.PropagationType;
51  import org.orekit.propagation.SpacecraftState;
52  import org.orekit.propagation.analytical.KeplerianPropagator;
53  import org.orekit.propagation.numerical.NumericalPropagator;
54  import org.orekit.propagation.sampling.OrekitStepHandler;
55  import org.orekit.propagation.sampling.OrekitStepInterpolator;
56  import org.orekit.time.AbsoluteDate;
57  import org.orekit.utils.Constants;
58  import org.orekit.utils.PVCoordinates;
59  import org.orekit.utils.PVCoordinatesProvider;
60  
61  
62  public class IntegratedEphemerisTest {
63  
64      @Test
65      public void testNormalKeplerIntegration() {
66  
67          // Keplerian propagator definition
68          KeplerianPropagator keplerEx = new KeplerianPropagator(initialOrbit);
69  
70          // Integrated ephemeris
71  
72          // Propagation
73          AbsoluteDate finalDate = initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY);
74          final EphemerisGenerator generator = numericalPropagator.getEphemerisGenerator();
75          numericalPropagator.setInitialState(new SpacecraftState(initialOrbit));
76          numericalPropagator.propagate(finalDate);
77          Assert.assertTrue(numericalPropagator.getCalls() < 3200);
78          BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
79  
80          // tests
81          for (int i = 1; i <= Constants.JULIAN_DAY; i++) {
82              AbsoluteDate intermediateDate = initialOrbit.getDate().shiftedBy(i);
83              SpacecraftState keplerIntermediateOrbit = keplerEx.propagate(intermediateDate);
84              SpacecraftState numericIntermediateOrbit = ephemeris.propagate(intermediateDate);
85              Vector3D kepPosition = keplerIntermediateOrbit.getPVCoordinates().getPosition();
86              Vector3D numPosition = numericIntermediateOrbit.getPVCoordinates().getPosition();
87              Assert.assertEquals(0, kepPosition.subtract(numPosition).getNorm(), 0.06);
88          }
89  
90          // test inv
91          AbsoluteDate intermediateDate = initialOrbit.getDate().shiftedBy(41589);
92          SpacecraftState keplerIntermediateOrbit = keplerEx.propagate(intermediateDate);
93          SpacecraftState state = keplerEx.propagate(finalDate);
94          numericalPropagator.setInitialState(state);
95          final EphemerisGenerator generator2 = numericalPropagator.getEphemerisGenerator();
96          numericalPropagator.propagate(initialOrbit.getDate());
97          BoundedPropagator invEphemeris = generator2.getGeneratedEphemeris();
98          SpacecraftState numericIntermediateOrbit = invEphemeris.propagate(intermediateDate);
99          Vector3D kepPosition = keplerIntermediateOrbit.getPVCoordinates().getPosition();
100         Vector3D numPosition = numericIntermediateOrbit.getPVCoordinates().getPosition();
101         Assert.assertEquals(0, kepPosition.subtract(numPosition).getNorm(), 10e-2);
102 
103     }
104 
105     @Test
106     public void testPartialDerivativesIssue16() {
107 
108         final String eqName = "derivatives";
109         final EphemerisGenerator generator = numericalPropagator.getEphemerisGenerator();
110         numericalPropagator.setOrbitType(OrbitType.CARTESIAN);
111         final MatricesHarvester harvester = numericalPropagator.setupMatricesComputation(eqName, null, null);
112         numericalPropagator.setInitialState(new SpacecraftState(initialOrbit));
113         numericalPropagator.propagate(initialOrbit.getDate().shiftedBy(3600.0));
114         BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
115         ephemeris.setStepHandler(new OrekitStepHandler() {
116 
117             public void handleStep(OrekitStepInterpolator interpolator) {
118                 SpacecraftState state = interpolator.getCurrentState();
119                 RealMatrix dYdY0 = harvester.getStateTransitionMatrix(state);
120                 harvester.getParametersJacobian(state); // no parameters, this is a no-op and should work
121                 RealMatrix deltaId = dYdY0.subtract(MatrixUtils.createRealIdentityMatrix(6));
122                 Assert.assertTrue(deltaId.getNorm1() >  100);
123                 Assert.assertTrue(deltaId.getNorm1() < 3100);
124             }
125 
126         });
127 
128         ephemeris.propagate(initialOrbit.getDate().shiftedBy(1800.0));
129 
130     }
131 
132     @Test
133     public void testGetFrame() {
134         // setup
135         AbsoluteDate finalDate = initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY);
136         final EphemerisGenerator generator = numericalPropagator.getEphemerisGenerator();
137         numericalPropagator.setInitialState(new SpacecraftState(initialOrbit));
138         numericalPropagator.propagate(finalDate);
139         Assert.assertTrue(numericalPropagator.getCalls() < 3200);
140         BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
141 
142         //action
143         Assert.assertNotNull(ephemeris.getFrame());
144         Assert.assertSame(ephemeris.getFrame(), numericalPropagator.getFrame());
145     }
146 
147     @Test
148     public void testIssue766() {
149 
150         // setup
151         AbsoluteDate finalDate = initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY);
152         final EphemerisGenerator generator = numericalPropagator.getEphemerisGenerator();
153         numericalPropagator.setInitialState(new SpacecraftState(initialOrbit));
154         numericalPropagator.propagate(finalDate);
155         Assert.assertTrue(numericalPropagator.getCalls() < 3200);
156         BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
157 
158         // verify
159         Assert.assertTrue(ephemeris.getAttitudeProvider() instanceof InertialProvider);
160 
161         // action
162         PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
163         ephemeris.setAttitudeProvider(new CelestialBodyPointed(FramesFactory.getEME2000(), sun, Vector3D.PLUS_K,
164                                                                Vector3D.PLUS_I, Vector3D.PLUS_K));
165         Assert.assertTrue(ephemeris.getAttitudeProvider() instanceof CelestialBodyPointed);
166 
167     }
168 
169     @Deprecated
170     @Test
171     public void testDeprecated() {
172 
173         EulerIntegrator integ = new EulerIntegrator(1.0);
174         DenseOutputModel dom = new DenseOutputModel();
175         integ.addStepHandler(dom);
176         ExpandableODE eode = new ExpandableODE(new OrdinaryDifferentialEquation() {
177             public int getDimension() { return 1; }
178             public double[] computeDerivatives(double t, double[] y) { return y; }
179         });
180         eode.addSecondaryEquations(new SecondaryODE() {
181             public int getDimension() { return 1; }
182             public double[] computeDerivatives(double t, double[] primary,
183                                                double[] primaryDot, double[] secondary) { return secondary; }
184         });
185         integ.integrate(eode, new ODEState(0.0, new double[1], new double[1][1]), 1.0);
186 
187         StateMapper mapper = new StateMapper(AbsoluteDate.ARBITRARY_EPOCH, Constants.EIGEN5C_EARTH_MU,
188                                              OrbitType.CARTESIAN, PositionAngle.TRUE,
189                                              new InertialProvider(FramesFactory.getEME2000()),
190                                              FramesFactory.getEME2000()) {
191             public void mapStateToArray(SpacecraftState state, double[] y, double[] yDot) {}
192             public SpacecraftState mapArrayToState(AbsoluteDate date, double[] y, double[] yDot, PropagationType type) {
193                 return null;
194             }
195         };
196 
197         try {
198             new IntegratedEphemeris(AbsoluteDate.ARBITRARY_EPOCH,
199                                     AbsoluteDate.PAST_INFINITY, AbsoluteDate.FUTURE_INFINITY,
200                                     mapper, PropagationType.OSCULATING,
201                                     dom, Collections.emptyMap(), Collections.emptyList(),
202                                     new String[] { "equation-1", "equation-2" });
203             Assert.fail("an exception should have been thrown");
204         } catch (OrekitInternalError oie) {
205             // expected as only one equation could be handled properly by this deprecated constructor
206         }
207 
208         IntegratedEphemeris ie = new IntegratedEphemeris(AbsoluteDate.ARBITRARY_EPOCH,
209                                                          AbsoluteDate.PAST_INFINITY, AbsoluteDate.FUTURE_INFINITY,
210                                                          mapper, PropagationType.OSCULATING,
211                                                          dom, Collections.emptyMap(), Collections.emptyList(),
212                                                          new String[] { "equation-1" });
213         Assert.assertNotNull(ie);
214 
215     }
216 
217     @Before
218     public void setUp() {
219 
220         Utils.setDataRoot("regular-data:potential/icgem-format");
221         GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", true));
222 
223         // Definition of initial conditions with position and velocity
224         Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
225         Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
226         double mu = 3.9860047e14;
227 
228         AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
229         initialOrbit =
230             new EquinoctialOrbit(new PVCoordinates(position, velocity),
231                                  FramesFactory.getEME2000(), initDate, mu);
232 
233         // Numerical propagator definition
234         double[] absTolerance = {
235             0.0001, 1.0e-11, 1.0e-11, 1.0e-8, 1.0e-8, 1.0e-8, 0.001
236         };
237         double[] relTolerance = {
238             1.0e-8, 1.0e-8, 1.0e-8, 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-7
239         };
240         AdaptiveStepsizeIntegrator integrator =
241             new DormandPrince853Integrator(0.001, 500, absTolerance, relTolerance);
242         integrator.setInitialStepSize(100);
243         numericalPropagator = new NumericalPropagator(integrator);
244 
245     }
246 
247     private Orbit initialOrbit;
248     private NumericalPropagator numericalPropagator;
249 
250 }