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