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.Field;
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.geometry.euclidean.threed.FieldVector3D;
24 import org.hipparchus.geometry.euclidean.threed.Vector3D;
25 import org.hipparchus.util.Binary64;
26 import org.hipparchus.util.Binary64Field;
27 import org.hipparchus.util.FastMath;
28 import org.hipparchus.util.MathUtils;
29 import org.junit.jupiter.api.Assertions;
30 import org.junit.jupiter.api.BeforeEach;
31 import org.junit.jupiter.api.Test;
32 import org.orekit.Utils;
33 import org.orekit.bodies.CelestialBody;
34 import org.orekit.bodies.CelestialBodyFactory;
35 import org.orekit.bodies.GeodeticPoint;
36 import org.orekit.bodies.OneAxisEllipsoid;
37 import org.orekit.errors.OrekitException;
38 import org.orekit.errors.OrekitMessages;
39 import org.orekit.frames.Frame;
40 import org.orekit.frames.FramesFactory;
41 import org.orekit.models.earth.atmosphere.data.JB2008SpaceEnvironmentData;
42 import org.orekit.time.AbsoluteDate;
43 import org.orekit.time.DateComponents;
44 import org.orekit.time.FieldAbsoluteDate;
45 import org.orekit.time.TimeComponents;
46 import org.orekit.time.TimeScale;
47 import org.orekit.time.TimeScalesFactory;
48 import org.orekit.utils.Constants;
49 import org.orekit.utils.IERSConventions;
50 import org.orekit.utils.PVCoordinatesProvider;
51
52 import java.text.ParseException;
53
54 class JB2008Test {
55
56 private static final TimeScale TT = TimeScalesFactory.getTT();
57
58 @Test
59 void testLegacy() {
60 final boolean print = false;
61
62 final CelestialBody sun = CelestialBodyFactory.getSun();
63 final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
64 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
65 Constants.WGS84_EARTH_FLATTENING, itrf);
66 final JB2008 atm = new JB2008(null, sun, earth);
67
68
69 final double[] D1950 = {20035.00454861111, 20035.50362, 20036.00468, 20036.50375,
70 20037.00000, 20037.50556, 20038.00088, 20038.50644,
71 20039.00543, 20039.50450, 20040.00000, 20040.50556};
72 final double[] SUNRA = {3.8826, 3.8914, 3.9001, 3.9089, 3.9176, 3.9265,
73 3.9352, 3.9441, 3.9529, 3.9618, 3.9706, 3.9795};
74 final double[] SUNDEC = {-0.2847, -0.2873, -0.2898, -0.2923, -0.2948, -0.2973,
75 -0.2998, -0.3022, -0.3046, -0.3070, -0.3094, -0.3117};
76 final double[] SATLON = {73.46, 218.68, 34.55, 46.25, 216.56, 32.00,
77 38.83, 213.67, 29.37, 38.12, 211.78, 26.64};
78 final double[] SATLAT = {-85.24, -18.65, 37.71, 74.36, -8.85, -39.64,
79 -51.93, -21.25, 46.43, 65.97, -21.31, -51.87};
80 final double[] SATALT = {398.91, 376.75, 373.45, 380.61, 374.03, 385.05,
81 389.83, 376.98, 374.56, 378.97, 377.76, 390.09};
82 final double[] F10 = {128.80, 128.80, 129.60, 129.60, 124.10, 124.10,
83 140.90, 140.90, 104.60, 104.60, 94.90, 94.90};
84 final double[] F10B = {105.60, 105.60, 105.60, 105.60, 105.60, 105.60,
85 105.60, 105.60, 105.70, 105.70, 105.90, 105.90};
86 final double[] S10 = {103.50, 103.50, 110.20, 110.20, 109.40, 109.40,
87 108.60, 108.60, 107.40, 107.40, 110.90, 110.90};
88 final double[] S10B = {103.80, 103.80, 103.80, 103.80, 103.80, 103.80,
89 103.70, 103.70, 103.70, 103.70, 103.70, 103.70};
90 final double[] XM10 = {110.90, 110.90, 115.60, 115.60, 110.00, 110.00,
91 110.00, 110.00, 106.90, 106.90, 102.20, 102.20};
92 final double[] XM10B = {106.90, 106.90, 106.90, 106.90, 106.90, 106.90,
93 107.00, 107.00, 107.10, 107.10, 107.10, 107.10};
94 final double[] Y10 = {127.90, 127.90, 125.90, 125.90, 127.70, 127.70,
95 125.60, 125.60, 126.60, 126.60, 126.80, 126.80};
96 final double[] Y10B = {112.90, 112.90, 112.90, 112.90, 113.00, 113.00,
97 113.20, 113.20, 113.20, 113.20, 113.30, 113.30};
98 final double[] DSTDTC = { 3., 80., 240., 307., 132., 40.,
99 327., 327., 118., 25., 85., 251.};
100
101
102 for (int i = 0; i < 12; i++) {
103 final double rho = atm.getDensity(MJD(D1950[i]), SUNRA[i], SUNDEC[i],
104 RAP(D1950[i], SATLON[i]),
105 FastMath.toRadians(SATLAT[i]),
106 SATALT[i] * 1000.,
107 F10[i], F10B[i], S10[i], S10B[i],
108 XM10[i], XM10B[i], Y10[i], Y10B[i], DSTDTC[i]);
109 checkLegacy(i, rho, print);
110 }
111
112 }
113
114
115 @Test
116 void testAltitude() {
117 final boolean print = false;
118
119 final InputParams ip = new InputParams();
120
121 final CelestialBody sun = CelestialBodyFactory.getSun();
122
123 final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
124 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
125 Constants.WGS84_EARTH_FLATTENING, itrf);
126 earth.setAngularThreshold(1e-10);
127
128 final JB2008 atm = new JB2008(ip, sun, earth);
129
130
131 final double[][] loc = {{-85.24, 73.46, 91.0e+3},
132 {-18.65, 218.68, 110.0e+3},
133 {-68.05, 145.28, 122.0e+3},
134 { 37.71, 34.55, 150.0e+3},
135 { 74.36, 46.25, 220.0e+3},
136 { -8.85, 216.56, 270.0e+3},
137 {-39.64, 32.00, 400.0e+3},
138 {-51.93, 38.83, 550.0e+3},
139 {-21.25, 213.67, 700.0e+3},
140 { 46.43, 29.37, 900.0e+3},
141 { 65.97, 38.12, 1200.0e+3},
142 {-21.31, 211.78, 1700.0e+3},
143 {-51.87, 26.64, 2300.0e+3}};
144
145
146 for (int i = 0; i < 13; i++) {
147
148 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(loc[i][0]),
149 FastMath.toRadians(loc[i][1]),
150 loc[i][2]);
151
152 final double rho = atm.getDensity(InputParams.TC[i], earth.transform(point), atm.getFrame());
153
154
155 checkAltitude(i, rho, print);
156 }
157 }
158
159 @Test
160 void testException() throws ParseException {
161
162 final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
163 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
164 Constants.WGS84_EARTH_FLATTENING, itrf);
165
166 final JB2008 atm = new JB2008(new InputParams(), CelestialBodyFactory.getSun(), earth);
167
168
169 try {
170 atm.getDensity(0., 0., 0., 0., 0., 89999.0, 0., 0., 0., 0., 0., 0., 0., 0., 0.);
171 Assertions.fail("an exception should have been thrown");
172 } catch (OrekitException oe) {
173 Assertions.assertEquals(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, oe.getSpecifier());
174 Assertions.assertEquals(89999.0, (Double) oe.getParts()[0], 1.0e-15);
175 Assertions.assertEquals(90000.0, (Double) oe.getParts()[1], 1.0e-15);
176 }
177
178 }
179
180 @Test
181 void testDensityField() {
182
183 final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
184 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
185 Constants.WGS84_EARTH_FLATTENING, itrf);
186
187 final JB2008 atm = new JB2008(new InputParams(), CelestialBodyFactory.getSun(), earth);
188
189 final AbsoluteDate date = InputParams.TC[4];
190
191 for (double alt = 100; alt < 1000; alt += 50) {
192 for (double lat = -1.2; lat < 1.2; lat += 0.4) {
193 for (double lon = 0; lon < 6.28; lon += 0.8) {
194
195 final GeodeticPoint point = new GeodeticPoint(lat, lon, alt * 1000.);
196 final Vector3D pos = earth.transform(point);
197 Field<Binary64> field = Binary64Field.getInstance();
198
199
200 final double rho = atm.getDensity(date, pos, itrf);
201 final Binary64 rho64 = atm.getDensity(new FieldAbsoluteDate<>(field, date),
202 new FieldVector3D<>(field, pos),
203 itrf);
204
205 Assertions.assertEquals(rho, rho64.getReal(), rho * 4.0e-10);
206
207 }
208 }
209 }
210
211 }
212
213 @Test
214 void testDensityGradient() {
215
216 final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
217 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
218 Constants.WGS84_EARTH_FLATTENING, itrf);
219
220 final JB2008 atm = new JB2008(new InputParams(), CelestialBodyFactory.getSun(), earth);
221
222 final AbsoluteDate date = InputParams.TC[6];
223
224
225 final double alt = 400.;
226 final double lat = 60.;
227 final double lon = -70.;
228 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(lat),
229 FastMath.toRadians(lon),
230 alt * 1000.);
231 final Vector3D pos = earth.transform(point);
232
233
234 DerivativeStructure zero = new DSFactory(1, 1).variable(0, 0.0);
235 FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(5, 10.0);
236 DerivativeStructure rhoX = differentiator.
237 differentiate((double x) -> {
238 try {
239 return atm.getDensity(date, new Vector3D(1, pos, x, Vector3D.PLUS_I), itrf);
240 } catch (OrekitException oe) {
241 return Double.NaN;
242 }
243 }). value(zero);
244 DerivativeStructure rhoY = differentiator.
245 differentiate((double y) -> {
246 try {
247 return atm.getDensity(date, new Vector3D(1, pos, y, Vector3D.PLUS_J), itrf);
248 } catch (OrekitException oe) {
249 return Double.NaN;
250 }
251 }). value(zero);
252 DerivativeStructure rhoZ = differentiator.
253 differentiate((double z) -> {
254 try {
255 return atm.getDensity(date, new Vector3D(1, pos, z, Vector3D.PLUS_K), itrf);
256 } catch (OrekitException oe) {
257 return Double.NaN;
258 }
259 }). value(zero);
260
261 DSFactory factory3 = new DSFactory(3, 1);
262 Field<DerivativeStructure> field = factory3.getDerivativeField();
263 final DerivativeStructure rhoDS = atm.getDensity(new FieldAbsoluteDate<>(field, date),
264 new FieldVector3D<>(factory3.variable(0, pos.getX()),
265 factory3.variable(1, pos.getY()),
266 factory3.variable(2, pos.getZ())),
267 itrf);
268
269 Assertions.assertEquals(rhoX.getValue(), rhoDS.getReal(), rhoX.getValue() * 2.0e-14);
270 Assertions.assertEquals(rhoY.getValue(), rhoDS.getReal(), rhoY.getValue() * 2.0e-14);
271 Assertions.assertEquals(rhoZ.getValue(), rhoDS.getReal(), rhoZ.getValue() * 2.0e-14);
272 Assertions.assertEquals(rhoX.getPartialDerivative(1),
273 rhoDS.getPartialDerivative(1, 0, 0),
274 FastMath.abs(6.0e-10 * rhoX.getPartialDerivative(1)));
275 Assertions.assertEquals(rhoY.getPartialDerivative(1),
276 rhoDS.getPartialDerivative(0, 1, 0),
277 FastMath.abs(6.0e-10 * rhoY.getPartialDerivative(1)));
278 Assertions.assertEquals(rhoZ.getPartialDerivative(1),
279 rhoDS.getPartialDerivative(0, 0, 1),
280 FastMath.abs(6.0e-10 * rhoY.getPartialDerivative(1)));
281
282 }
283
284 @Test
285 void testComparisonWithReference() {
286
287
288
289
290
291
292 final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
293 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
294 Constants.WGS84_EARTH_FLATTENING, itrf);
295
296
297 final int year = 2004;
298 final int month = 1;
299 final int day = 2;
300 final int hour = 12;
301 final int min = 0;
302 final double sec = 0.0;
303 final double lat = FastMath.toRadians(45.0);
304 final double lon = FastMath.toRadians(45.0);
305 final double alt = 250.0e3;
306
307
308 final JB2008SpaceEnvironmentData JBData = new JB2008SpaceEnvironmentData("SOLFSMY_trunc.txt", "DTCFILE_trunc.TXT");
309 final JB2008 atm = new JB2008(JBData, CelestialBodyFactory.getSun(), earth);
310
311
312 final GeodeticPoint point = new GeodeticPoint(lat, lon, alt);
313 final Vector3D pos = earth.transform(point);
314 final AbsoluteDate date = new AbsoluteDate(year, month, day, hour, min, sec, TimeScalesFactory.getUTC());
315 final double density = atm.getDensity(date, pos, itrf);
316
317
318 final double ref = 6.6862e-11;
319 Assertions.assertEquals(ref, density, 1.0e-15);
320
321 }
322
323
324
325
326
327 private double MJD(final double d1950) {
328 return d1950 + 33281.0;
329 }
330
331
332
333
334
335
336 private double RAP(final double d1950, final double satLon) {
337 double theta;
338 final double nbday = FastMath.floor(d1950);
339 if (nbday < 7305.) {
340 theta = 1.7294446614 + 1.72027915246e-2 * nbday + 6.3003880926 * (d1950 - nbday);
341 } else {
342 final double ts70 = d1950 - 7305.;
343 final double ids70 = FastMath.floor(ts70);
344 final double tfrac = ts70 - ids70;
345 theta = 1.73213438565 + 1.720279169407e-2 * ids70 +
346 (1.720279169407e-2 + MathUtils.TWO_PI) * tfrac +
347 5.0755141943e-15 * ts70 * ts70;
348 }
349 theta = MathUtils.normalizeAngle(theta, FastMath.PI);
350
351 return MathUtils.normalizeAngle(theta + FastMath.toRadians(satLon), 0.);
352 }
353
354
355
356
357
358
359 private void checkLegacy(final int id, final double rho, final boolean print) {
360 final double[] rhoRef = {0.18730056e-11, 0.25650339e-11, 0.57428913e-11,
361 0.83266893e-11, 0.82238726e-11, 0.48686457e-11,
362 0.67210914e-11, 0.74215571e-11, 0.31821075e-11,
363 0.29553578e-11, 0.64122627e-11, 0.79559727e-11};
364 final double dRho = 2.e-5;
365 final int nb = id + 1;
366 if (print) {
367 System.out.printf("Case #%d\n", nb);
368 System.out.printf("Rho: %12.5e %12.5e\n", rhoRef[id], rho);
369 } else {
370 Assertions.assertEquals(rhoRef[id], rho, rhoRef[id] * dRho);
371 }
372
373 }
374
375
376
377
378
379
380 private void checkAltitude(final int id, final double rho, final boolean print) {
381 final double[] rhoRef = {0.27945654e-05, 0.94115202e-07, 0.15025977e-07, 0.21128330e-08,
382 0.15227435e-09, 0.54609767e-10, 0.45899746e-11, 0.14922800e-12,
383 0.17392987e-13, 0.35250121e-14, 0.13482414e-14, 0.77684879e-15,
384 0.19900569e-15};
385 final double dRho = 3.e-4;
386 final int nb = id + 1;
387 if (print) {
388 System.out.printf("Case #%d\n", nb);
389 System.out.printf("Rho: %12.5e %12.5e\n", rhoRef[id], rho);
390 } else {
391 Assertions.assertEquals(rhoRef[id], rho, rhoRef[id] * dRho);
392 }
393
394 }
395
396 @BeforeEach
397 public void setUp() {
398 Utils.setDataRoot("regular-data:atmosphere");
399 }
400
401 private static class InputParams implements JB2008InputParameters {
402
403
404 private static final long serialVersionUID = 5091441542522297257L;
405
406
407 public static final AbsoluteDate[] TC = new AbsoluteDate[] {
408 new AbsoluteDate(new DateComponents(2003, 312), new TimeComponents( 0, 6, 33.0), TT),
409 new AbsoluteDate(new DateComponents(2003, 312), new TimeComponents(12, 5, 13.0), TT),
410 new AbsoluteDate(new DateComponents(2003, 312), new TimeComponents(18, 45, 3.0), TT),
411 new AbsoluteDate(new DateComponents(2003, 313), new TimeComponents( 0, 6, 44.0), TT),
412 new AbsoluteDate(new DateComponents(2003, 313), new TimeComponents(12, 5, 24.0), TT),
413 new AbsoluteDate(new DateComponents(2003, 314), new TimeComponents( 0, 0, 0.0), TT),
414 new AbsoluteDate(new DateComponents(2003, 314), new TimeComponents(12, 8, 0.0), TT),
415 new AbsoluteDate(new DateComponents(2003, 315), new TimeComponents( 0, 1, 16.0), TT),
416 new AbsoluteDate(new DateComponents(2003, 315), new TimeComponents(12, 9, 16.0), TT),
417 new AbsoluteDate(new DateComponents(2003, 316), new TimeComponents( 0, 7, 49.0), TT),
418 new AbsoluteDate(new DateComponents(2003, 316), new TimeComponents(12, 6, 29.0), TT),
419 new AbsoluteDate(new DateComponents(2003, 317), new TimeComponents( 0, 0, 0.0), TT),
420 new AbsoluteDate(new DateComponents(2003, 317), new TimeComponents(12, 8, 0.0), TT)
421 };
422
423
424 private static final double[] F10 = new double[] {
425 91.00, 91.00, 91.00, 92.70, 92.70, 93.00,
426 93.00, 94.60, 94.60, 95.60, 95.60, 98.70, 98.70
427 };
428
429
430 private static final double[] F10B = new double[] {
431 137.10, 137.10, 137.10, 136.90, 136.90, 136.80,
432 136.80, 136.70, 136.70, 136.70, 136.70, 136.90, 136.90
433 };
434
435
436 private static final double[] S10 = new double[] {
437 108.80, 108.80, 108.80, 104.20, 104.20, 102.60,
438 102.60, 100.30, 100.30, 99.50, 99.50, 101.20, 101.20
439 };
440
441
442 private static final double[] S10B = new double[] {
443 123.80, 123.80, 123.80, 123.70, 123.70, 123.60,
444 123.60, 123.50, 123.50, 123.50, 123.50, 123.60, 123.60
445 };
446
447
448 private static final double[] XM10 = new double[] {
449 116.70, 116.70, 116.70, 109.60, 109.60, 100.20,
450 100.20, 97.00, 97.00, 95.40, 95.40, 95.40, 95.40
451 };
452
453
454 private static final double[] XM10B = new double[] {
455 128.50, 128.50, 128.50, 128.00, 128.00, 127.70,
456 127.70, 127.60, 127.60, 127.60, 127.60, 127.70, 127.70
457 };
458
459
460 private static final double[] Y10 = new double[] {
461 168.00, 168.00, 168.00, 147.90, 147.90, 131.60,
462 131.60, 122.60, 122.60, 114.30, 114.30, 112.70, 112.70
463 };
464
465
466 private static final double[] Y10B = new double[] {
467 138.60, 138.60, 138.60, 138.40, 138.40, 138.10,
468 138.10, 137.90, 137.90, 137.90, 137.90, 137.80, 137.80
469 };
470
471
472 private static final double[] DSTDTC = new double[] {
473 43.00, 30.00, 67.00, 90.00, 87.00, 115.00,
474 114.00, 106.00, 148.00, 148.00, 105.00, 146.00, 119.00
475 };
476
477
478 public InputParams() {
479 }
480
481 @Override
482 public AbsoluteDate getMinDate() {
483 return TC[0];
484 }
485
486 @Override
487 public AbsoluteDate getMaxDate() {
488 return TC[TC.length - 1];
489 }
490
491 @Override
492 public double getF10(AbsoluteDate date)
493 {
494 for (int i = 0; i < TC.length; i++) {
495 if (date.equals(TC[i])) {
496 return F10[i];
497 }
498 }
499 return Double.NaN;
500 }
501
502 @Override
503 public double getF10B(AbsoluteDate date)
504 {
505 for (int i = 0; i < TC.length; i++) {
506 if (date.equals(TC[i])) {
507 return F10B[i];
508 }
509 }
510 return Double.NaN;
511 }
512
513 @Override
514 public double getS10(AbsoluteDate date)
515 {
516 for (int i = 0; i < TC.length; i++) {
517 if (date.equals(TC[i])) {
518 return S10[i];
519 }
520 }
521 return Double.NaN;
522 }
523
524 @Override
525 public double getS10B(AbsoluteDate date)
526 {
527 for (int i = 0; i < TC.length; i++) {
528 if (date.equals(TC[i])) {
529 return S10B[i];
530 }
531 }
532 return Double.NaN;
533 }
534
535 @Override
536 public double getXM10(AbsoluteDate date)
537 {
538 for (int i = 0; i < TC.length; i++) {
539 if (date.equals(TC[i])) {
540 return XM10[i];
541 }
542 }
543 return Double.NaN;
544 }
545
546 @Override
547 public double getXM10B(AbsoluteDate date)
548 {
549 for (int i = 0; i < TC.length; i++) {
550 if (date.equals(TC[i])) {
551 return XM10B[i];
552 }
553 }
554 return Double.NaN;
555 }
556
557 @Override
558 public double getY10(AbsoluteDate date)
559 {
560 for (int i = 0; i < TC.length; i++) {
561 if (date.equals(TC[i])) {
562 return Y10[i];
563 }
564 }
565 return Double.NaN;
566 }
567
568 @Override
569 public double getY10B(AbsoluteDate date)
570 {
571 for (int i = 0; i < TC.length; i++) {
572 if (date.equals(TC[i])) {
573 return Y10B[i];
574 }
575 }
576 return Double.NaN;
577 }
578
579 @Override
580 public double getDSTDTC(AbsoluteDate date)
581 {
582 for (int i = 0; i < TC.length; i++) {
583 if (date.equals(TC[i])) {
584 return DSTDTC[i];
585 }
586 }
587 return Double.NaN;
588 }
589
590 }
591
592 }