1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.models.earth.atmosphere;
18
19 import org.hipparchus.analysis.UnivariateFunction;
20 import org.hipparchus.analysis.differentiation.DSFactory;
21 import org.hipparchus.analysis.differentiation.DerivativeStructure;
22 import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
23 import org.hipparchus.complex.Complex;
24 import org.hipparchus.complex.ComplexField;
25 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
26 import org.hipparchus.geometry.euclidean.threed.Vector3D;
27 import org.junit.jupiter.api.AfterEach;
28 import org.junit.jupiter.api.Assertions;
29 import org.junit.jupiter.api.BeforeEach;
30 import org.junit.jupiter.api.Test;
31 import org.orekit.Utils;
32 import org.orekit.bodies.CelestialBody;
33 import org.orekit.bodies.CelestialBodyFactory;
34 import org.orekit.bodies.GeodeticPoint;
35 import org.orekit.bodies.OneAxisEllipsoid;
36 import org.orekit.errors.OrekitException;
37 import org.orekit.frames.Frame;
38 import org.orekit.frames.FramesFactory;
39 import org.orekit.time.AbsoluteDate;
40 import org.orekit.time.DateComponents;
41 import org.orekit.time.FieldAbsoluteDate;
42 import org.orekit.time.TimeComponents;
43 import org.orekit.time.TimeScale;
44 import org.orekit.time.TimeScalesFactory;
45 import org.orekit.utils.Constants;
46 import org.orekit.utils.PVCoordinatesProvider;
47
48 class HarrisPriesterTest {
49
50
51 private CelestialBody sun;
52
53
54 private OneAxisEllipsoid earth;
55
56
57 private Frame earthFrame;
58
59
60 private TimeScale utc;
61
62
63 private AbsoluteDate date;
64
65 @Test
66 void testStandard() {
67
68 final HarrisPriester hp = new HarrisPriester(sun, earth);
69
70
71 final GeodeticPoint point = new GeodeticPoint(0, 0, 500000.);
72 final Vector3D pos = earth.transform(point);
73
74
75 final double rho = hp.getDensity(date, pos, earthFrame);
76
77 Assertions.assertEquals(3.9237E-13, rho, 1.0e-17);
78
79 }
80
81 @Test
82 void testParameterN() {
83
84 final HarrisPriester hp = new HarrisPriester(sun, earth);
85
86
87 final GeodeticPoint point = new GeodeticPoint(0, 0, 500000.);
88 final Vector3D pos = earth.transform(point);
89
90
91 final double rho4 = hp.getDensity(date, pos, earthFrame);
92
93
94 final HarrisPriester hp2 = new HarrisPriester(sun, earth, 2);
95
96
97 final double rho2 = hp2.getDensity(date, pos, earthFrame);
98
99 final HarrisPriester hp6 = new HarrisPriester(sun, earth, 6);
100
101
102 final double rho6 = hp6.getDensity(date, pos, earthFrame);
103
104 final double c2Psi2 = 0.02163787;
105
106 Assertions.assertEquals(c2Psi2, (rho6-rho2)/(rho4-rho2) - 1., 1.e-8);
107
108 }
109
110 @Test
111 void testMaxAlt() {
112
113 final HarrisPriester hp = new HarrisPriester(sun, earth);
114
115
116 final GeodeticPoint point = new GeodeticPoint(0, 0, 1500000.);
117 final Vector3D pos = earth.transform(point);
118
119
120 final double rho = hp.getDensity(date, pos, earthFrame);
121
122 Assertions.assertEquals(0.0, rho, 0.0);
123 }
124
125 @Test
126 void testUserTab() {
127
128 final double[][] userTab = {
129 {100000., 4.974e+02, 4.974e+02},
130 {110000., 7.800e+01, 7.800e+01},
131 {120000., 2.490e+01, 2.400e+01},
132 {130000., 8.377e+00, 8.710e+00},
133 {140000., 3.899e+00, 4.059e+00},
134 {150000., 2.122e+00, 2.215e+00},
135 {160000., 1.263e+00, 1.344e+00},
136 {170000., 8.008e-01, 8.758e-01},
137 {180000., 5.283e-01, 6.010e-01},
138 {190000., 3.618e-01, 4.297e-01},
139 {200000., 2.557e-01, 3.162e-01},
140 {210000., 1.839e-01, 2.396e-01},
141 {220000., 1.341e-01, 1.853e-01},
142 {230000., 9.949e-02, 1.455e-01},
143 {240000., 7.488e-02, 1.157e-01},
144 {250000., 5.709e-02, 9.308e-02},
145 {260000., 4.403e-02, 7.555e-02},
146 {270000., 3.430e-02, 6.182e-02},
147 {280000., 2.697e-02, 5.095e-02},
148 {290000., 2.139e-02, 4.226e-02},
149 {300000., 1.708e-02, 3.526e-02},
150 {320000., 1.099e-02, 2.511e-02},
151 {340000., 7.214e-03, 1.819e-02},
152 {360000., 4.824e-03, 1.337e-02},
153 {380000., 3.274e-03, 9.955e-03},
154 {400000., 2.249e-03, 7.492e-03},
155 {420000., 1.558e-03, 5.684e-03},
156 {440000., 1.091e-03, 4.355e-03},
157 {460000., 7.701e-04, 3.362e-03},
158 {480000., 5.474e-04, 2.612e-03},
159 {500000., 3.916e-04, 2.042e-03},
160 {520000., 2.819e-04, 1.605e-03},
161 {540000., 2.042e-04, 1.267e-03},
162 {560000., 1.488e-04, 1.005e-03},
163 {580000., 1.092e-04, 7.997e-04},
164 {600000., 8.070e-05, 6.390e-04},
165 {620000., 6.012e-05, 5.123e-04},
166 {640000., 4.519e-05, 4.121e-04},
167 {660000., 3.430e-05, 3.325e-04},
168 {680000., 2.632e-05, 2.691e-04},
169 {700000., 2.043e-05, 2.185e-04},
170 {720000., 1.607e-05, 1.779e-04},
171 {740000., 1.281e-05, 1.452e-04},
172 {760000., 1.036e-05, 1.190e-04},
173 {780000., 8.496e-06, 9.776e-05},
174 {800000., 7.069e-06, 8.059e-05},
175 {850000., 4.800e-06, 5.500e-05},
176 {900000., 3.300e-06, 3.700e-05},
177 {950000., 2.450e-06, 2.400e-05},
178 {1000000., 1.900e-06, 1.700e-05},
179 {1100000., 1.180e-06, 8.700e-06},
180 {1200000., 7.500e-07, 4.800e-06},
181 {1300000., 5.300e-07, 3.200e-06},
182 {1400000., 4.100e-07, 2.000e-06},
183 {1500000., 2.900e-07, 1.350e-06},
184 {1600000., 2.000e-07, 9.500e-07},
185 {1700000., 1.600e-07, 7.700e-07},
186 {1800000., 1.200e-07, 6.300e-07},
187 {1900000., 9.600e-08, 5.200e-07},
188 {2000000., 7.300e-08, 4.400e-07}
189 };
190
191
192 final GeodeticPoint point = new GeodeticPoint(0, 0, 1500000.);
193 final Vector3D pos = earth.transform(point);
194
195 final HarrisPriester hp = new HarrisPriester(sun, earth, userTab);
196
197
198 final double rho = hp.getDensity(date, pos, earthFrame);
199
200 Assertions.assertEquals(2.9049E-7, rho, 1.0e-11);
201
202 final HarrisPriester hp6 = new HarrisPriester(sun, earth, userTab, 6);
203 final double rho6 = hp6.getDensity(date, pos, earthFrame);
204
205 final HarrisPriester hp2 = new HarrisPriester(sun, earth, userTab, 2);
206 final double rho2 = hp2.getDensity(date, pos, earthFrame);
207
208 final double c2Psi2 = 0.02163787;
209
210 Assertions.assertEquals(c2Psi2, (rho6-rho2)/(rho-rho2) - 1., 1.0e-8);
211
212 }
213
214 @Test
215 void testOutOfRange() {
216 Assertions.assertThrows(OrekitException.class, () -> {
217 final HarrisPriester hp = new HarrisPriester(sun, earth);
218
219
220 final GeodeticPoint point = new GeodeticPoint(0, 0, 50000.);
221 final Vector3D pos = earth.transform(point);
222
223
224 hp.getDensity(date, pos, earthFrame);
225 });
226 }
227
228 @Test
229 void testGetDensityField() {
230
231 final ComplexField field = ComplexField.getInstance();
232 final HarrisPriester hp = new HarrisPriester(sun, earth);
233 final Frame frame = FramesFactory.getGCRF();
234 final FieldAbsoluteDate<Complex> fieldDate = new FieldAbsoluteDate<>(field, date);
235
236 for (double height = 300e3; height <= 2000e3; height += 100) {
237 final Vector3D position = new Vector3D(Constants.EGM96_EARTH_EQUATORIAL_RADIUS + height, 0., 0.);
238 final FieldVector3D<Complex> fieldPosition = new FieldVector3D<>(field, position);
239 final Complex fieldDensity = hp.getDensity(fieldDate, fieldPosition, frame);
240 final double density = hp.getDensity(date, position, frame);
241 Assertions.assertEquals(density, fieldDensity.getReal());
242 }
243 }
244
245 @Test
246 void testVelocityDerivative() {
247 final Frame eme2000 = FramesFactory.getEME2000();
248 final HarrisPriester hp = new HarrisPriester(sun, earth);
249 final Vector3D pos = earth.getBodyFrame().getStaticTransformTo(eme2000, date).
250 transformPosition(earth.transform(new GeodeticPoint(-1.7, 4.2, 987654.321)));
251 double dP = 100.0;
252 double dVxdX = gradientComponent(hp, pos, Vector3D.PLUS_I, dP, eme2000, v -> v.getX());
253 double dVxdY = gradientComponent(hp, pos, Vector3D.PLUS_J, dP, eme2000, v -> v.getX());
254 double dVxdZ = gradientComponent(hp, pos, Vector3D.PLUS_K, dP, eme2000, v -> v.getX());
255 double dVydX = gradientComponent(hp, pos, Vector3D.PLUS_I, dP, eme2000, v -> v.getY());
256 double dVydY = gradientComponent(hp, pos, Vector3D.PLUS_J, dP, eme2000, v -> v.getY());
257 double dVydZ = gradientComponent(hp, pos, Vector3D.PLUS_K, dP, eme2000, v -> v.getY());
258 double dVzdX = gradientComponent(hp, pos, Vector3D.PLUS_I, dP, eme2000, v -> v.getZ());
259 double dVzdY = gradientComponent(hp, pos, Vector3D.PLUS_J, dP, eme2000, v -> v.getZ());
260 double dVzdZ = gradientComponent(hp, pos, Vector3D.PLUS_K, dP, eme2000, v -> v.getZ());
261
262 DSFactory factory = new DSFactory(3, 1);
263 FieldVector3D<DerivativeStructure> dsPos = new FieldVector3D<>(factory.variable(0, pos.getX()),
264 factory.variable(1, pos.getY()),
265 factory.variable(2, pos.getZ()));
266 FieldVector3D<DerivativeStructure> dsVel = hp.getVelocity(new FieldAbsoluteDate<>(factory.getDerivativeField(),
267 date),
268 dsPos, eme2000);
269 Assertions.assertEquals(dVxdX, dsVel.getX().getPartialDerivative(1, 0, 0), 1.0e-16);
270 Assertions.assertEquals(dVxdY, dsVel.getX().getPartialDerivative(0, 1, 0), 1.0e-16);
271 Assertions.assertEquals(dVxdZ, dsVel.getX().getPartialDerivative(0, 0, 1), 1.0e-16);
272 Assertions.assertEquals(dVydX, dsVel.getY().getPartialDerivative(1, 0, 0), 1.0e-16);
273 Assertions.assertEquals(dVydY, dsVel.getY().getPartialDerivative(0, 1, 0), 1.0e-16);
274 Assertions.assertEquals(dVydZ, dsVel.getY().getPartialDerivative(0, 0, 1), 1.0e-16);
275 Assertions.assertEquals(dVzdX, dsVel.getZ().getPartialDerivative(1, 0, 0), 1.0e-16);
276 Assertions.assertEquals(dVzdY, dsVel.getZ().getPartialDerivative(0, 1, 0), 1.0e-16);
277 Assertions.assertEquals(dVzdZ, dsVel.getZ().getPartialDerivative(0, 0, 1), 1.0e-16);
278
279 }
280
281 private double gradientComponent(final Atmosphere atm, final Vector3D position, final Vector3D direction,
282 final double dP, final Frame frame, final ComponentGetter getter) {
283 FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(5, dP);
284 UnivariateFunction f = delta -> {
285 try {
286 return getter.get(atm.getVelocity(date, new Vector3D(1, position, delta, direction), frame));
287 } catch (OrekitException oe) {
288 return Double.NaN;
289 }
290 };
291 return differentiator.differentiate(f).value(new DSFactory(1, 1).variable(0, 0.0)).getPartialDerivative(1);
292 }
293
294 private interface ComponentGetter {
295 double get(Vector3D v);
296 }
297
298 @BeforeEach
299 public void setUp() {
300 Utils.setDataRoot("regular-data");
301 sun = CelestialBodyFactory.getSun();
302
303 earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
304 earth = new OneAxisEllipsoid(6378136.460, 1.0 / 298.257222101, earthFrame);
305
306
307 utc = TimeScalesFactory.getUTC();
308 date = new AbsoluteDate(new DateComponents(2003, 03, 21), new TimeComponents(1, 0, 0.), utc);
309 }
310
311 @AfterEach
312 public void tearDown() {
313 utc = null;
314 }
315
316 }