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 org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.Field;
21  import org.hipparchus.geometry.euclidean.threed.Vector3D;
22  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
23  import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
24  import org.hipparchus.util.Decimal64Field;
25  import org.junit.After;
26  import org.junit.Assert;
27  import org.junit.Before;
28  import org.junit.Test;
29  import org.orekit.Utils;
30  import org.orekit.forces.gravity.potential.GravityFieldFactory;
31  import org.orekit.forces.gravity.potential.SHMFormatReader;
32  import org.orekit.frames.FramesFactory;
33  import org.orekit.orbits.EquinoctialOrbit;
34  import org.orekit.orbits.Orbit;
35  import org.orekit.orbits.OrbitType;
36  import org.orekit.propagation.FieldSpacecraftState;
37  import org.orekit.propagation.SpacecraftState;
38  import org.orekit.propagation.numerical.FieldNumericalPropagator;
39  import org.orekit.propagation.numerical.NumericalPropagator;
40  import org.orekit.propagation.semianalytical.dsst.FieldDSSTPropagator;
41  import org.orekit.time.AbsoluteDate;
42  import org.orekit.time.FieldAbsoluteDate;
43  import org.orekit.utils.PVCoordinates;
44  
45  @Deprecated
46  public class FieldAdditionalEquationsTest {
47  
48      private double                       mu;
49      private AbsoluteDate                 initDate;
50      private SpacecraftState              initialState;
51      private double[][]                   tolerance;
52  
53      /** Test for issue #401
54       *  with a numerical propagator */
55      @Test
56      public void testInitNumerical() {
57          doTestInitNumerical(Decimal64Field.getInstance());
58      }
59  
60      /** Test for issue #401
61       *  with a DSST propagator */
62      @Test
63      public void testInitDSST() {
64          doTestInitDSST(Decimal64Field.getInstance());
65      }
66  
67      @Test
68      public void testResetState() {
69          doTestResetState(Decimal64Field.getInstance());
70      }
71  
72      private <T extends CalculusFieldElement<T>> void doTestInitNumerical(Field<T> field) {
73          // setup
74          final double reference = 1.25;
75          final double rate      = 1.5;
76          final double dt        = 600.0;
77          Linear<T> linear = new Linear<>("linear", reference, rate);
78          Assert.assertFalse(linear.wasCalled());
79  
80          // action
81          AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200,
82                                                                                                tolerance[0], tolerance[1]);
83          integrator.setInitialStepSize(60);
84          FieldNumericalPropagator<T> propagatorNumerical = new FieldNumericalPropagator<>(field, integrator);
85          propagatorNumerical.setInitialState(new FieldSpacecraftState<>(field, initialState).
86                                              addAdditionalState(linear.getName(), field.getZero().newInstance(reference)));
87          propagatorNumerical.addAdditionalEquations(linear);
88          FieldSpacecraftState<T> finalState = propagatorNumerical.propagate(new FieldAbsoluteDate<>(field, initDate).shiftedBy(dt));
89  
90          // verify
91          Assert.assertTrue(linear.wasCalled());
92          Assert.assertEquals(reference + dt * rate, finalState.getAdditionalState(linear.getName())[0].getReal(), 1.0e-10);
93  
94      }
95  
96      private <T extends CalculusFieldElement<T>> void doTestInitDSST(Field<T> field) {
97          // setup
98          final double reference = 3.5;
99          final double rate      = 1.5;
100         final double dt        = 600.0;
101         Linear<T> linear = new Linear<>("linear", reference, rate);
102         Assert.assertFalse(linear.wasCalled());
103 
104         // action
105         AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200,
106                                                                                               tolerance[0], tolerance[1]);
107         integrator.setInitialStepSize(60);
108         FieldDSSTPropagator<T> propagatorDSST = new FieldDSSTPropagator<>(field, integrator);
109         propagatorDSST.setInitialState(new FieldSpacecraftState<>(field, initialState).
110                                        addAdditionalState(linear.getName(), field.getZero().newInstance(reference)));
111         propagatorDSST.addAdditionalEquations(linear);
112         FieldSpacecraftState<T> finalState = propagatorDSST.propagate(new FieldAbsoluteDate<>(field, initDate).shiftedBy(dt));
113 
114         // verify
115         Assert.assertTrue(linear.wasCalled());
116         Assert.assertEquals(reference + dt * rate, finalState.getAdditionalState(linear.getName())[0].getReal(), 1.0e-10);
117 
118     }
119 
120     private <T extends CalculusFieldElement<T>> void doTestResetState(Field<T> field) {
121         // setup
122         final double reference1 = 3.5;
123         final double rate1      = 1.5;
124         Linear<T> linear1 = new Linear<>("linear-1", reference1, rate1);
125         Assert.assertFalse(linear1.wasCalled());
126         final double reference2 = 4.5;
127         final double rate2      = 1.25;
128         Linear<T> linear2 = new Linear<>("linear-2", reference2, rate2);
129         Assert.assertFalse(linear2.wasCalled());
130         final double dt = 600;
131 
132         // action
133         AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200,
134                                                                                               tolerance[0], tolerance[1]);
135         integrator.setInitialStepSize(60);
136         FieldNumericalPropagator<T> propagatorNumerical = new FieldNumericalPropagator<>(field, integrator);
137         propagatorNumerical.setInitialState(new FieldSpacecraftState<>(field, initialState).
138                                             addAdditionalState(linear1.getName(), field.getZero().newInstance(reference1)).
139                                             addAdditionalState(linear2.getName(), field.getZero().newInstance(reference2)));
140         propagatorNumerical.addAdditionalEquations(linear1);
141         propagatorNumerical.addAdditionalEquations(linear2);
142         FieldSpacecraftState<T> finalState = propagatorNumerical.propagate(new FieldAbsoluteDate<>(field, initDate).shiftedBy(dt));
143 
144         // verify
145         Assert.assertTrue(linear1.wasCalled());
146         Assert.assertTrue(linear2.wasCalled());
147         Assert.assertEquals(reference1 + dt * rate1, finalState.getAdditionalState(linear1.getName())[0].getReal(), 1.0e-10);
148         Assert.assertEquals(reference2 + dt * rate2, finalState.getAdditionalState(linear2.getName())[0].getReal(), 1.0e-10);
149 
150     }
151 
152     @Before
153     public void setUp() {
154         Utils.setDataRoot("regular-data:potential/shm-format");
155         GravityFieldFactory.addPotentialCoefficientsReader(new SHMFormatReader("^eigen_cg03c_coef$", false));
156         mu  = GravityFieldFactory.getUnnormalizedProvider(0, 0).getMu();
157         final Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
158         final Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
159         initDate = AbsoluteDate.J2000_EPOCH;
160         final Orbit orbit = new EquinoctialOrbit(new PVCoordinates(position,  velocity),
161                                                  FramesFactory.getEME2000(), initDate, mu);
162         initialState = new SpacecraftState(orbit);
163         tolerance = NumericalPropagator.tolerances(0.001, orbit, OrbitType.EQUINOCTIAL);
164     }
165 
166     @After
167     public void tearDown() {
168         initDate     = null;
169         initialState = null;
170         tolerance    = null;
171     }
172 
173     private static class Linear<T extends CalculusFieldElement<T>> implements FieldAdditionalEquations<T> {
174 
175         private String  name;
176         private double  expectedAtInit;
177         private double  rate;
178         private boolean called;
179 
180         public Linear(final String name, final double expectedAtInit, final double rate) {
181             this.name           = name;
182             this.expectedAtInit = expectedAtInit;
183             this.rate           = rate;
184             this.called         = false;
185         }
186 
187         @Override
188         public void init(FieldSpacecraftState<T> initiaState, FieldAbsoluteDate<T> target) {
189             Assert.assertEquals(expectedAtInit, initiaState.getAdditionalState(getName())[0].getReal(), 1.0e-15);
190             called = true;
191         }
192 
193         @Override
194         public T[] computeDerivatives(FieldSpacecraftState<T> s, T[] pDot) {
195             pDot[0] = s.getDate().getField().getZero().newInstance(rate);
196             return null;
197         }
198 
199         @Override
200         public String getName() {
201             return name;
202         }
203 
204         public boolean wasCalled() {
205             return called;
206         }
207 
208     }
209 
210 }