1   /* Copyright 2002-2025 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.numerical;
18  
19  import org.hipparchus.geometry.euclidean.threed.Vector3D;
20  import org.hipparchus.linear.MatrixUtils;
21  import org.hipparchus.linear.RealMatrix;
22  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
23  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
24  import org.junit.jupiter.api.AfterEach;
25  import org.junit.jupiter.api.Assertions;
26  import org.junit.jupiter.api.BeforeEach;
27  import org.junit.jupiter.api.Test;
28  import org.orekit.errors.OrekitException;
29  import org.orekit.errors.OrekitMessages;
30  import org.orekit.forces.maneuvers.Maneuver;
31  import org.orekit.forces.maneuvers.propulsion.BasicConstantThrustPropulsionModel;
32  import org.orekit.forces.maneuvers.propulsion.PropulsionModel;
33  import org.orekit.forces.maneuvers.trigger.DateBasedManeuverTriggers;
34  import org.orekit.frames.FramesFactory;
35  import org.orekit.orbits.KeplerianOrbit;
36  import org.orekit.orbits.Orbit;
37  import org.orekit.orbits.OrbitType;
38  import org.orekit.orbits.PositionAngleType;
39  import org.orekit.propagation.SpacecraftState;
40  import org.orekit.propagation.ToleranceProvider;
41  import org.orekit.time.AbsoluteDate;
42  import org.orekit.utils.AbsolutePVCoordinates;
43  import org.orekit.utils.Constants;
44  import org.orekit.utils.DoubleArrayDictionary;
45  
46  import java.util.List;
47  
48  public class NumericalPropagationHarvesterTest {
49  
50      @Test
51      public void testNullStmName() {
52           try {
53              propagator.setupMatricesComputation(null, null, null);
54              Assertions.fail("an exception should have been thrown");
55          } catch (OrekitException oe) {
56              Assertions.assertEquals(OrekitMessages.NULL_ARGUMENT, oe.getSpecifier());
57              Assertions.assertEquals("stmName", oe.getParts()[0]);
58          }
59      }
60  
61      @Test
62      public void testUnknownStmName() {
63          NumericalPropagationHarvester harvester =
64                          (NumericalPropagationHarvester) propagator.setupMatricesComputation("stm",
65                                                                                              MatrixUtils.createRealIdentityMatrix(6),
66                                                                                              new DoubleArrayDictionary());
67          Assertions.assertNull(harvester.getStateTransitionMatrix(propagator.getInitialState()));
68      }
69  
70      @Test
71      public void testUnknownColumnName() {
72          NumericalPropagationHarvester harvester =
73                          (NumericalPropagationHarvester) propagator.setupMatricesComputation("stm",
74                                                                                              MatrixUtils.createRealIdentityMatrix(6),
75                                                                                              new DoubleArrayDictionary());
76          Assertions.assertNull(harvester.getParametersJacobian(propagator.getInitialState()));
77      }
78  
79      @Test
80      public void testDefaultNonNullInitialJacobian() {
81          NumericalPropagationHarvester harvester =
82                          (NumericalPropagationHarvester) propagator.setupMatricesComputation("stm",
83                                                                                              MatrixUtils.createRealIdentityMatrix(6),
84                                                                                              new DoubleArrayDictionary());
85          Assertions.assertNotNull(harvester.getInitialJacobianColumn("xyz"));
86      }
87  
88      @Test
89      public void testInitialStmCartesian() {
90          doTestInitialStm(OrbitType.CARTESIAN, 0.0);
91      }
92  
93      @Test
94      public void testInitialStmKeplerian() {
95          doTestInitialStm(OrbitType.KEPLERIAN, 2160.746);
96      }
97  
98      @Test
99      public void testInitialStmAbsPV() {
100         SpacecraftState state = propagator.getInitialState();
101         SpacecraftState absPV =
102                         new SpacecraftState(new AbsolutePVCoordinates(state.getFrame(),
103                                                                       state.getPVCoordinates()));
104         propagator.setInitialState(absPV);
105         doTestInitialStm(null, 0.0);
106     }
107 
108     @Test
109     public void testColumnsNames() {
110 
111         NumericalPropagationHarvester harvester =
112                         (NumericalPropagationHarvester) propagator.setupMatricesComputation("stm",
113                                                                                             MatrixUtils.createRealIdentityMatrix(6),
114                                                                                             new DoubleArrayDictionary());
115         Assertions.assertTrue(harvester.getJacobiansColumnsNames().isEmpty());
116 
117         DateBasedManeuverTriggers triggers = new DateBasedManeuverTriggers("apogee_boost", propagator.getInitialState().getDate().shiftedBy(60.0), 120.0);
118         PropulsionModel propulsion = new BasicConstantThrustPropulsionModel(400.0, 350.0, Vector3D.PLUS_I, "ABM-");
119         propagator.addForceModel(new Maneuver(null, triggers, propulsion));
120         Assertions.assertTrue(harvester.getJacobiansColumnsNames().isEmpty());
121 
122         triggers.getParametersDrivers().get(1).setSelected(true);
123         propulsion.getParametersDrivers().get(0).setSelected(true);
124         List<String> columnsNames = harvester.getJacobiansColumnsNames();
125         Assertions.assertEquals(2, columnsNames.size());
126         Assertions.assertEquals("SpanABM-" + BasicConstantThrustPropulsionModel.THRUST + Integer.toString(0), columnsNames.get(0));
127         Assertions.assertEquals("Spanapogee_boost_STOP" + Integer.toString(0), columnsNames.get(1));
128 
129     }
130 
131     private void doTestInitialStm(OrbitType type, double deltaId) {
132         PositionAngleType angle = PositionAngleType.TRUE;
133         NumericalPropagationHarvester harvester =
134                         (NumericalPropagationHarvester) propagator.setupMatricesComputation("stm", null, null);
135         propagator.setOrbitType(type);
136         propagator.setPositionAngleType(angle);
137         double[] p = new double[36];
138         for (int i = 0; i < p.length; i += 7) {
139             p[i] = 1.0;
140         }
141         SpacecraftState s = propagator.getInitialState().addAdditionalData(harvester.getStmName(), p);
142         RealMatrix stm = harvester.getStateTransitionMatrix(s);
143         Assertions.assertEquals(deltaId, stm.subtract(MatrixUtils.createRealIdentityMatrix(6)).getNorm1(), 1.0e-3);
144         Assertions.assertEquals(type, harvester.getOrbitType());
145         Assertions.assertEquals(angle, harvester.getPositionAngleType());
146     }
147 
148     @BeforeEach
149     public void setUp() {
150         Orbit initialOrbit =
151                         new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngleType.TRUE,
152                                            FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
153                                            Constants.EIGEN5C_EARTH_MU);
154         double minStep = 0.0001;
155         double maxStep = 60;
156         double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(0.001).getTolerances(initialOrbit, initialOrbit.getType());
157         AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tolerances[0], tolerances[1]);
158         integrator.setInitialStepSize(1.0);
159         propagator = new NumericalPropagator(integrator);
160         propagator.setInitialState(new SpacecraftState(initialOrbit));
161     }
162 
163     @AfterEach
164     public void tearDown() {
165         propagator = null;
166     }
167 
168     private NumericalPropagator propagator;
169 
170 }