1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.conversion;
18
19 import org.hipparchus.analysis.MultivariateVectorFunction;
20 import org.hipparchus.exception.LocalizedCoreFormats;
21 import org.hipparchus.exception.MathIllegalArgumentException;
22 import org.hipparchus.geometry.euclidean.threed.Vector3D;
23 import org.hipparchus.linear.ArrayRealVector;
24 import org.hipparchus.linear.RealMatrix;
25 import org.hipparchus.linear.RealVector;
26 import org.hipparchus.random.RandomGenerator;
27 import org.hipparchus.random.Well19937a;
28 import org.hipparchus.util.FastMath;
29 import org.hipparchus.util.Pair;
30 import org.junit.jupiter.api.Assertions;
31 import org.junit.jupiter.api.BeforeEach;
32 import org.junit.jupiter.api.Test;
33 import org.orekit.Utils;
34 import org.orekit.bodies.OneAxisEllipsoid;
35 import org.orekit.forces.ForceModel;
36 import org.orekit.forces.drag.DragForce;
37 import org.orekit.forces.drag.DragSensitive;
38 import org.orekit.forces.drag.IsotropicDrag;
39 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
40 import org.orekit.forces.gravity.NewtonianAttraction;
41 import org.orekit.forces.gravity.potential.GravityFieldFactory;
42 import org.orekit.frames.FramesFactory;
43 import org.orekit.models.earth.atmosphere.Atmosphere;
44 import org.orekit.models.earth.atmosphere.SimpleExponentialAtmosphere;
45 import org.orekit.orbits.EquinoctialOrbit;
46 import org.orekit.orbits.Orbit;
47 import org.orekit.orbits.OrbitType;
48 import org.orekit.orbits.PositionAngleType;
49 import org.orekit.propagation.Propagator;
50 import org.orekit.propagation.SpacecraftState;
51 import org.orekit.time.AbsoluteDate;
52 import org.orekit.time.TimeScalesFactory;
53 import org.orekit.utils.Constants;
54 import org.orekit.utils.IERSConventions;
55 import org.orekit.utils.PVCoordinates;
56 import org.orekit.utils.ParameterDriver;
57
58 import java.io.IOException;
59 import java.lang.reflect.InvocationTargetException;
60 import java.lang.reflect.Method;
61 import java.text.ParseException;
62 import java.util.ArrayList;
63 import java.util.List;
64
65 public class JacobianPropagatorConverterTest {
66
67 private double mu;
68 private double dP;
69
70 private Orbit orbit;
71 private ForceModel gravity;
72 private ForceModel drag;
73 private Atmosphere atmosphere;
74 private double crossSection;
75
76 @Test
77 public void testDerivativesNothing() {
78 try {
79 doTestDerivatives(1.0, 1.0);
80 Assertions.fail("an exception should have been thrown");
81 } catch (MathIllegalArgumentException miae) {
82 Assertions.assertEquals(LocalizedCoreFormats.AT_LEAST_ONE_COLUMN, miae.getSpecifier());
83 }
84 }
85
86 @Test
87 public void testDerivativesOrbitOnly() {
88 doTestDerivatives(4.8e-9, 3.5e-12,
89 "Px", "Py", "Pz", "Vx", "Vy", "Vz");
90 }
91
92 @Test
93 public void testDerivativesPositionAndDrag() {
94 doTestDerivatives(5.1e-9, 4.8e-12,
95 "Px", "Py", "Pz", DragSensitive.DRAG_COEFFICIENT);
96 }
97
98 @Test
99 public void testDerivativesDrag() {
100 doTestDerivatives(3.2e-9, 3.2e-12,
101 DragSensitive.DRAG_COEFFICIENT);
102 }
103
104 @Test
105 public void testDerivativesCentralAttraction() {
106 doTestDerivatives(3.6e-9, 4.0e-12,
107 NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT);
108 }
109
110 @Test
111 public void testDerivativesAllParameters() {
112 doTestDerivatives(1.1e-8, 1.1e-11,
113 "Px", "Py", "Pz", "Vx", "Vy", "Vz",
114 DragSensitive.DRAG_COEFFICIENT,
115 NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT);
116 }
117
118 private void doTestDerivatives(double tolP, double tolV, String... names) {
119
120
121
122
123 NumericalPropagatorBuilder builder =
124 new NumericalPropagatorBuilder(OrbitType.CARTESIAN.convertType(orbit),
125 new LutherIntegratorBuilder(10.0),
126 PositionAngleType.TRUE, dP);
127 builder.setMass(200.0);
128 builder.addForceModel(drag);
129 builder.addForceModel(gravity);
130
131
132
133 RandomGenerator random = new Well19937a(0xe67f19c1a678d037l);
134 List<ParameterDriver> all = new ArrayList<ParameterDriver>();
135 for (final ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
136 all.add(driver);
137 }
138 for (final ParameterDriver driver : builder.getPropagationParametersDrivers().getDrivers()) {
139 all.add(driver);
140 }
141 double[] normalized = new double[names.length];
142 List<ParameterDriver> selected = new ArrayList<ParameterDriver>(names.length);
143 int index = 0;
144 for (final ParameterDriver driver : all) {
145 boolean found = false;
146 for (final String name : names) {
147 if (name.equals(driver.getName())) {
148 found = true;
149 normalized[index++] = driver.getNormalizedValue(new AbsoluteDate()) + (2 * random.nextDouble() - 1);
150 selected.add(driver);
151 }
152 }
153 driver.setSelected(found);
154 }
155
156
157
158 final List<SpacecraftState> sample = new ArrayList<SpacecraftState>();
159 Propagator propagator = builder.buildPropagator(normalized);
160 propagator.setStepHandler(60.0, currentState -> sample.add(currentState));
161 propagator.propagate(orbit.getDate().shiftedBy(600.0), orbit.getDate().shiftedBy(4200.0));
162
163 JacobianPropagatorConverter fitter = new JacobianPropagatorConverter(builder, 1.0e-3, 5000);
164 try {
165 Method setSample = AbstractPropagatorConverter.class.getDeclaredMethod("setSample", List.class);
166 setSample.setAccessible(true);
167 setSample.invoke(fitter, sample);
168 } catch (NoSuchMethodException | SecurityException | IllegalAccessException |
169 IllegalArgumentException | InvocationTargetException e) {
170 Assertions.fail(e.getLocalizedMessage());
171 }
172
173 MultivariateVectorFunction f = fitter.getObjectiveFunction();
174 Pair<RealVector, RealMatrix> p = fitter.getModel().value(new ArrayRealVector(normalized));
175
176
177
178 RealMatrix m = p.getSecond();
179 double h = 10.0;
180 double[] shifted = normalized.clone();
181 double maxErrorP = 0;
182 double maxErrorV = 0;
183 for (int j = 0; j < selected.size(); ++j) {
184 shifted[j] = normalized[j] + 2.0 * h;
185 double[] valueP2 = f.value(shifted);
186 shifted[j] = normalized[j] + 1.0 * h;
187 double[] valueP1 = f.value(shifted);
188 shifted[j] = normalized[j] - 1.0 * h;
189 double[] valueM1 = f.value(shifted);
190 shifted[j] = normalized[j] - 2.0 * h;
191 double[] valueM2 = f.value(shifted);
192 shifted[j] = normalized[j];
193 for (int i = 0; i < valueP2.length; ++i) {
194 double d = (8 * (valueP1[i] - valueM1[i]) - (valueP2[i] - valueM2[i])) / (12 * h);
195 if (i % 6 < 3) {
196
197 maxErrorP = FastMath.max(maxErrorP, FastMath.abs(m.getEntry(i, j) - d));
198 } else {
199
200 maxErrorV = FastMath.max(maxErrorV, FastMath.abs(m.getEntry(i, j) - d));
201 }
202 }
203 }
204 Assertions.assertEquals(0.0, maxErrorP, tolP);
205 Assertions.assertEquals(0.0, maxErrorV, tolV);
206
207 }
208
209 @BeforeEach
210 public void setUp() throws IOException, ParseException {
211
212 Utils.setDataRoot("regular-data:potential/shm-format");
213 gravity = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true),
214 GravityFieldFactory.getNormalizedProvider(2, 0));
215 mu = gravity.getParameterDriver(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT).getValue();
216 dP = 1.0;
217
218
219 final Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6).normalize()
220 .scalarMultiply(Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 300e3);
221 final Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
222 final AbsoluteDate initDate = new AbsoluteDate(2010, 10, 10, 10, 10, 10.0, TimeScalesFactory.getUTC());
223 orbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
224 FramesFactory.getEME2000(), initDate, mu);
225
226 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
227 Constants.WGS84_EARTH_FLATTENING,
228 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
229 earth.setAngularThreshold(1.e-7);
230 atmosphere = new SimpleExponentialAtmosphere(earth, 0.0004, 42000.0, 7500.0);
231 final double dragCoef = 2.0;
232 crossSection = 25.0;
233 drag = new DragForce(atmosphere, new IsotropicDrag(crossSection, dragCoef));
234
235 }
236
237 }
238