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 java.lang.reflect.Constructor;
20 import java.lang.reflect.InvocationTargetException;
21 import java.lang.reflect.Method;
22
23 import org.hipparchus.Field;
24 import org.hipparchus.CalculusFieldElement;
25 import org.hipparchus.analysis.differentiation.DSFactory;
26 import org.hipparchus.analysis.differentiation.DerivativeStructure;
27 import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
28 import org.hipparchus.exception.LocalizedCoreFormats;
29 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
30 import org.hipparchus.geometry.euclidean.threed.Vector3D;
31 import org.hipparchus.random.RandomGenerator;
32 import org.hipparchus.random.Well19937a;
33 import org.hipparchus.util.Decimal64;
34 import org.hipparchus.util.Decimal64Field;
35 import org.hipparchus.util.FastMath;
36 import org.hipparchus.util.Precision;
37 import org.junit.Assert;
38 import org.junit.Before;
39 import org.junit.Test;
40 import org.orekit.Utils;
41 import org.orekit.bodies.CelestialBodyFactory;
42 import org.orekit.bodies.GeodeticPoint;
43 import org.orekit.bodies.OneAxisEllipsoid;
44 import org.orekit.errors.OrekitException;
45 import org.orekit.errors.OrekitMessages;
46 import org.orekit.frames.Frame;
47 import org.orekit.frames.FramesFactory;
48 import org.orekit.time.AbsoluteDate;
49 import org.orekit.time.DateComponents;
50 import org.orekit.time.FieldAbsoluteDate;
51 import org.orekit.time.TimeComponents;
52 import org.orekit.time.TimeScalesFactory;
53 import org.orekit.utils.Constants;
54 import org.orekit.utils.IERSConventions;
55 import org.orekit.utils.PVCoordinatesProvider;
56
57
58 public class NRLMSISE00Test {
59
60 @Test
61 public void testLegacy() throws
62 NoSuchMethodException, SecurityException, InstantiationException,
63 IllegalAccessException, IllegalArgumentException, InvocationTargetException {
64
65
66 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
67 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
68 Constants.WGS84_EARTH_FLATTENING, itrf);
69 NRLMSISE00 atm = new NRLMSISE00(null, null, earth);
70
71
72 final int doy = 172;
73 final double sec = 29000.;
74 final double alt = 400.;
75 final double lat = 60.;
76 final double lon = -70.;
77 final double hl = 16.;
78 final double f107a = 150.;
79 final double f107 = 150.;
80 double[] ap = {4., 100., 100., 100., 100., 100., 100.};
81 final boolean print = false;
82
83 Method gtd7 = getOutputClass().getDeclaredMethod("gtd7", Double.TYPE);
84 gtd7.setAccessible(true);
85
86
87 final Object out1 = createOutput(atm, doy, sec, lat, lon, hl, f107a, f107, ap);
88 gtd7.invoke(out1, alt);
89 checkLegacy(1, out1, print);
90
91
92 final int doy2 = 81;
93 final Object out2 = createOutput(atm, doy2, sec, lat, lon, hl, f107a, f107, ap);
94 gtd7.invoke(out2, alt);
95 checkLegacy(2, out2, print);
96
97
98 final double sec3 = 75000.;
99 final double alt3 = 1000.;
100 final Object out3 = createOutput(atm, doy, sec3, lat, lon, hl, f107a, f107, ap);
101 gtd7.invoke(out3, alt3);
102 checkLegacy(3, out3, print);
103
104
105 final double alt4 = 100.;
106 final Object out4 = createOutput(atm, doy, sec, lat, lon, hl, f107a, f107, ap);
107 gtd7.invoke(out4, alt4);
108 checkLegacy(4, out4, print);
109
110
111 final double lat5 = 0.;
112 final Object out5 = createOutput(atm, doy, sec, lat5, lon, hl, f107a, f107, ap);
113 gtd7.invoke(out5, alt);
114 checkLegacy(5, out5, print);
115
116
117 final double lon6 = 0.;
118 final Object out6 = createOutput(atm, doy, sec, lat, lon6, hl, f107a, f107, ap);
119 gtd7.invoke(out6, alt);
120 checkLegacy(6, out6, print);
121
122
123 final double hl7 = 4.;
124 final Object out7 = createOutput(atm, doy, sec, lat, lon, hl7, f107a, f107, ap);
125 gtd7.invoke(out7, alt);
126 checkLegacy(7, out7, print);
127
128
129 final double f107a8 = 70.;
130 final Object out8 = createOutput(atm, doy, sec, lat, lon, hl, f107a8, f107, ap);
131 gtd7.invoke(out8, alt);
132 checkLegacy(8, out8, print);
133
134
135 final double f1079 = 180.;
136 final Object out9 = createOutput(atm, doy, sec, lat, lon, hl, f107a, f1079, ap);
137 gtd7.invoke(out9, alt);
138 checkLegacy(9, out9, print);
139
140
141 ap[0] = 40.;
142 final Object out10 = createOutput(atm, doy, sec, lat, lon, hl, f107a, f107, ap);
143 gtd7.invoke(out10, alt);
144 checkLegacy(10, out10, print);
145 ap[0] = 4.;
146
147
148 final double alt11 = 0.;
149 final Object out11 = createOutput(atm, doy, sec, lat, lon, hl, f107a, f107, ap);
150 gtd7.invoke(out11, alt11);
151 checkLegacy(11, out11, print);
152
153
154 final double alt12 = 10.;
155 final Object out12 = createOutput(atm, doy, sec, lat, lon, hl, f107a, f107, ap);
156 gtd7.invoke(out12, alt12);
157 checkLegacy(12, out12, print);
158
159
160 final double alt13 = 30.;
161 final Object out13 = createOutput(atm, doy, sec, lat, lon, hl, f107a, f107, ap);
162 gtd7.invoke(out13, alt13);
163 checkLegacy(13, out13, print);
164
165
166 final double alt14 = 50.;
167 final Object out14 = createOutput(atm, doy, sec, lat, lon, hl, f107a, f107, ap);
168 gtd7.invoke(out14, alt14);
169 checkLegacy(14, out14, print);
170
171
172 final double alt15 = 70.;
173 final Object out15 = createOutput(atm, doy, sec, lat, lon, hl, f107a, f107, ap);
174 gtd7.invoke(out15, alt15);
175 checkLegacy(15, out15, print);
176
177
178 NRLMSISE00 otherAtm = atm.withSwitch(9, -1);
179 final Object out16 = createOutput(otherAtm, doy, sec, lat, lon, hl, f107a, f107, ap);
180 gtd7.invoke(out16, alt);
181 checkLegacy(16, out16, print);
182
183
184 final double alt17 = 100.;
185 final Object out17 = createOutput(otherAtm, doy, sec, lat, lon, hl, f107a, f107, ap);
186 gtd7.invoke(out17, alt17);
187 checkLegacy(17, out17, print);
188
189 }
190
191 @Test
192 public void testDensity() throws
193 InstantiationException, IllegalAccessException,
194 IllegalArgumentException, InvocationTargetException,
195 NoSuchMethodException, SecurityException {
196
197 final InputParams ip = new InputParams();
198
199 final PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
200
201 final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
202 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
203 Constants.WGS84_EARTH_FLATTENING, itrf);
204
205 final NRLMSISE00 atm = new NRLMSISE00(ip, sun, earth).withSwitch(9, -1);
206
207 final AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 172),
208 new TimeComponents(29000.),
209 TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true));
210
211 final double alt = 400.;
212 final double lat = 60.;
213 final double lon = -70.;
214 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(lat),
215 FastMath.toRadians(lon),
216 alt * 1000.);
217 final Vector3D pos = earth.transform(point);
218
219
220 try {
221 atm.getDensity(date.shiftedBy(2 * Constants.JULIAN_YEAR), pos, itrf);
222 Assert.fail("an exception should have been thrown");
223 } catch (OrekitException oe) {
224 Assert.assertEquals(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE, oe.getSpecifier());
225 }
226
227 final double rho = atm.getDensity(date, pos, itrf);
228 final double lst = 29000. / 3600. - 70. / 15.;
229 final double[] ap = {4., 100., 100., 100., 100., 100., 100.};
230
231 Class<?> outputClass = getOutputClass();
232 Constructor<?> cons = outputClass.getDeclaredConstructor(NRLMSISE00.class,
233 Integer.TYPE,
234 Double.TYPE,
235 Double.TYPE,
236 Double.TYPE,
237 Double.TYPE,
238 Double.TYPE,
239 Double.TYPE,
240 double[].class);
241 cons.setAccessible(true);
242 Method gtd7d = outputClass.getDeclaredMethod("gtd7d", Double.TYPE);
243 gtd7d.setAccessible(true);
244 Method getDensity = outputClass.getDeclaredMethod("getDensity", Integer.TYPE);
245 getDensity.setAccessible(true);
246
247 final Object out = createOutput(atm, 172, 29000., 60., -70, lst, 150., 150., ap);
248 gtd7d.invoke(out, 400.0);
249 Assert.assertEquals(rho, ((Double) getDensity.invoke(out, 5)).doubleValue(), rho * 1.e-3);
250
251 }
252
253 @Test
254 public void testDensityField() {
255
256 final InputParams ip = new InputParams();
257
258 final PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
259
260 final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
261 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
262 Constants.WGS84_EARTH_FLATTENING, itrf);
263
264 final NRLMSISE00 atm = new NRLMSISE00(ip, sun, earth);
265
266 final AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 172),
267 new TimeComponents(29000.),
268 TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true));
269
270 final double alt = 400.;
271 final double lat = 60.;
272 final double lon = -70.;
273 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(lat),
274 FastMath.toRadians(lon),
275 alt * 1000.);
276 final Vector3D pos = earth.transform(point);
277 Field<Decimal64> field = Decimal64Field.getInstance();
278
279
280 try {
281 atm.getDensity(new FieldAbsoluteDate<>(field, date).shiftedBy(2 * Constants.JULIAN_YEAR),
282 new FieldVector3D<>(field.getOne(), pos),
283 itrf);
284 Assert.fail("an exception should have been thrown");
285 } catch (OrekitException oe) {
286 Assert.assertEquals(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE, oe.getSpecifier());
287 }
288
289 final double rho = atm.getDensity(date, pos, itrf);
290 final Decimal64 rho64 = atm.getDensity(new FieldAbsoluteDate<>(field, date),
291 new FieldVector3D<>(field.getOne(), pos),
292 itrf);
293
294 Assert.assertEquals(rho, rho64.getReal(), rho * 2.0e-13);
295
296 }
297
298 @Test
299 public void testDensityGradient() {
300
301 final InputParams ip = new InputParams();
302
303 final PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
304
305 final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
306 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
307 Constants.WGS84_EARTH_FLATTENING, itrf);
308
309 final NRLMSISE00 atm = new NRLMSISE00(ip, sun, earth);
310
311 final AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 172),
312 new TimeComponents(29000.),
313 TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true));
314
315 final double alt = 400.;
316 final double lat = 60.;
317 final double lon = -70.;
318 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(lat),
319 FastMath.toRadians(lon),
320 alt * 1000.);
321 final Vector3D pos = earth.transform(point);
322
323
324 DerivativeStructure zero = new DSFactory(1, 1).variable(0, 0.0);
325 FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(5, 10.0);
326 DerivativeStructure rhoX = differentiator.
327 differentiate((double x) -> {
328 try {
329 return atm.getDensity(date, new Vector3D(1, pos, x, Vector3D.PLUS_I), itrf);
330 } catch (OrekitException oe) {
331 return Double.NaN;
332 }
333 }). value(zero);
334 DerivativeStructure rhoY = differentiator.
335 differentiate((double y) -> {
336 try {
337 return atm.getDensity(date, new Vector3D(1, pos, y, Vector3D.PLUS_J), itrf);
338 } catch (OrekitException oe) {
339 return Double.NaN;
340 }
341 }). value(zero);
342 DerivativeStructure rhoZ = differentiator.
343 differentiate((double z) -> {
344 try {
345 return atm.getDensity(date, new Vector3D(1, pos, z, Vector3D.PLUS_K), itrf);
346 } catch (OrekitException oe) {
347 return Double.NaN;
348 }
349 }). value(zero);
350
351 DSFactory factory3 = new DSFactory(3, 1);
352 Field<DerivativeStructure> field = factory3.getDerivativeField();
353 final DerivativeStructure rhoDS = atm.getDensity(new FieldAbsoluteDate<>(field, date),
354 new FieldVector3D<>(factory3.variable(0, pos.getX()),
355 factory3.variable(1, pos.getY()),
356 factory3.variable(2, pos.getZ())),
357 itrf);
358
359 Assert.assertEquals(rhoX.getValue(), rhoDS.getReal(), rhoX.getValue() * 2.0e-13);
360 Assert.assertEquals(rhoY.getValue(), rhoDS.getReal(), rhoY.getValue() * 2.0e-13);
361 Assert.assertEquals(rhoZ.getValue(), rhoDS.getReal(), rhoZ.getValue() * 2.0e-13);
362 Assert.assertEquals(rhoX.getPartialDerivative(1),
363 rhoDS.getPartialDerivative(1, 0, 0),
364 FastMath.abs(2.0e-10 * rhoX.getPartialDerivative(1)));
365 Assert.assertEquals(rhoY.getPartialDerivative(1),
366 rhoDS.getPartialDerivative(0, 1, 0),
367 FastMath.abs(2.0e-10 * rhoY.getPartialDerivative(1)));
368 Assert.assertEquals(rhoZ.getPartialDerivative(1),
369 rhoDS.getPartialDerivative(0, 0, 1),
370 FastMath.abs(2.0e-10 * rhoY.getPartialDerivative(1)));
371
372 }
373
374 @Test
375 public void testWrongNumberLow() {
376 try {
377 new NRLMSISE00(null, null, null).withSwitch(0, 17);
378 Assert.fail("an exception should have been thrown");
379 } catch (OrekitException oe) {
380 Assert.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
381 Assert.assertEquals( 0, oe.getParts()[0]);
382 Assert.assertEquals( 1, oe.getParts()[1]);
383 Assert.assertEquals(23, oe.getParts()[2]);
384 }
385 }
386
387 @Test
388 public void testWrongNumberHigh() {
389 try {
390 new NRLMSISE00(null, null, null).withSwitch(24, 17);
391 Assert.fail("an exception should have been thrown");
392 } catch (OrekitException oe) {
393 Assert.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
394 Assert.assertEquals(24, oe.getParts()[0]);
395 Assert.assertEquals( 1, oe.getParts()[1]);
396 Assert.assertEquals(23, oe.getParts()[2]);
397 }
398 }
399
400 @Test
401 public void testGlobe7SwitchesOn() {
402 RandomGenerator random = new Well19937a(0xb9d06451353d23cbl);
403 NRLMSISE00 atm = new NRLMSISE00(null, null, null);
404 for (int i = 1; i <= 23; ++i) {
405 atm = atm.withSwitch(i, 1);
406 }
407 doTestDoubleMethod(atm, random, "globe7", 2.0e-14, 2.2e-16);
408 }
409
410 @Test
411 public void testGlobe7SwitchesOff() {
412 RandomGenerator random = new Well19937a(0x778b486a40464b8fl);
413 NRLMSISE00 atm = new NRLMSISE00(null, null, null);
414 for (int i = 1; i <= 23; ++i) {
415 atm = atm.withSwitch(i, 0);
416 }
417 doTestDoubleMethod(atm, random, "globe7", 1.0e-50, 1.0e-50);
418 }
419
420 @Test
421 public void testGlobe7SwitchesRandom() {
422 RandomGenerator random = new Well19937a(0xe20a69235cc9583dl);
423 NRLMSISE00 atm = new NRLMSISE00(null, null, null);
424 for (int i = 1; i <= 23; ++i) {
425 atm = atm.withSwitch(i, random.nextInt(3) - 2);
426 }
427 doTestDoubleMethod(atm, random, "globe7", 3.0e-14, 4.0e-16);
428 }
429
430 @Test
431 public void testGlob7sSwitchesOn() {
432 RandomGenerator random = new Well19937a(0xc7c218fabec5e98cl);
433 NRLMSISE00 atm = new NRLMSISE00(null, null, null);
434 for (int i = 1; i <= 23; ++i) {
435 atm = atm.withSwitch(i, 1);
436 }
437 doTestDoubleMethod(atm, random, "glob7s", 4.0e-15, 9.0e-16);
438 }
439
440 @Test
441 public void testGlob7sSwitchesOff() {
442 RandomGenerator random = new Well19937a(0x141f7aa933299a83l);
443 NRLMSISE00 atm = new NRLMSISE00(null, null, null);
444 for (int i = 1; i <= 23; ++i) {
445 atm = atm.withSwitch(i, 0);
446 }
447 doTestDoubleMethod(atm, random, "glob7s", 1.0e-50, 1.0e-50);
448 }
449
450 @Test
451 public void testGlob7sSwitchesRandom() {
452 RandomGenerator random = new Well19937a(0x3671893ce741fc5cl);
453 NRLMSISE00 atm = new NRLMSISE00(null, null, null);
454 for (int i = 1; i <= 23; ++i) {
455 atm = atm.withSwitch(i, random.nextInt(3) - 2);
456 }
457 doTestDoubleMethod(atm, random, "glob7s", 1.0e-50, 1.0e-50);
458 }
459
460 @Test
461 public void testgts7SwitchesOn() {
462 RandomGenerator random = new Well19937a(0xb6dcf73ed5e5d985l);
463 NRLMSISE00 atm = new NRLMSISE00(null, null, null);
464 for (int i = 1; i <= 23; ++i) {
465 atm = atm.withSwitch(i, 1);
466 }
467 doTestVoidMethod(atm, random, "gts7", 1.0e-50, 6.0e-14);
468 }
469
470 @Test
471 public void testgts7SwitchesOff() {
472 RandomGenerator random = new Well19937a(0x0c953641bea0f6d2l);
473 NRLMSISE00 atm = new NRLMSISE00(null, null, null);
474 for (int i = 1; i <= 23; ++i) {
475 atm = atm.withSwitch(i, 0);
476 }
477 doTestVoidMethod(atm, random, "gts7", 1.0e-50, 6.0e-14);
478 }
479
480 @Test
481 public void testgts7SwitchesRandom() {
482 RandomGenerator random = new Well19937a(0x7347cacb946cb93bl);
483 NRLMSISE00 atm = new NRLMSISE00(null, null, null);
484 for (int i = 1; i <= 23; ++i) {
485 atm = atm.withSwitch(i, random.nextInt(3) - 2);
486 }
487 doTestVoidMethod(atm, random, "gts7", 1.0e-50, 6.0e-14);
488 }
489
490 @Test
491 public void testgtd7SwitchesOn() {
492 RandomGenerator random = new Well19937a(0x3439206bdd4dff5dl);
493 NRLMSISE00 atm = new NRLMSISE00(null, null, null);
494 for (int i = 1; i <= 23; ++i) {
495 atm = atm.withSwitch(i, 1);
496 }
497 doTestVoidMethod(atm, random, "gtd7", 5.0e-16, 3.0e-14);
498 }
499
500 @Test
501 public void testgtd7SwitchesOff() {
502 RandomGenerator random = new Well19937a(0x3dc1f824e1033d1bl);
503 NRLMSISE00 atm = new NRLMSISE00(null, null, null);
504 for (int i = 1; i <= 23; ++i) {
505 atm = atm.withSwitch(i, 0);
506 }
507 doTestVoidMethod(atm, random, "gtd7", 1.0e-50, 3.0e-14);
508 }
509
510 @Test
511 public void testgtd7SwitchesRandom() {
512 RandomGenerator random = new Well19937a(0xa12175ef0b689b04l);
513 NRLMSISE00 atm = new NRLMSISE00(null, null, null);
514 for (int i = 1; i <= 23; ++i) {
515 atm = atm.withSwitch(i, random.nextInt(3) - 2);
516 }
517 doTestVoidMethod(atm, random, "gtd7", 1.0e-50, 3.0e-14);
518 }
519
520 @Test
521 public void testgtd7dSwitchesOn() {
522 RandomGenerator random = new Well19937a(0x4bbb424422a1b909l);
523 NRLMSISE00 atm = new NRLMSISE00(null, null, null);
524 for (int i = 1; i <= 23; ++i) {
525 atm = atm.withSwitch(i, 1);
526 }
527 doTestVoidMethod(atm, random, "gtd7d", 3.0e-16, 3.0e-14);
528 }
529
530 @Test
531 public void testgtd7dSwitchesOff() {
532 RandomGenerator random = new Well19937a(0x7f6da37655e30103l);
533 NRLMSISE00 atm = new NRLMSISE00(null, null, null);
534 for (int i = 1; i <= 23; ++i) {
535 atm = atm.withSwitch(i, 0);
536 }
537 doTestVoidMethod(atm, random, "gtd7d", 1.0e-50, 2.0e-14);
538 }
539
540 @Test
541 public void testgtd7dSwitchesRandom() {
542 RandomGenerator random = new Well19937a(0x4a75e29ddf23ccd7l);
543 NRLMSISE00 atm = new NRLMSISE00(null, null, null);
544 for (int i = 1; i <= 23; ++i) {
545 atm = atm.withSwitch(i, random.nextInt(3) - 2);
546 }
547 doTestVoidMethod(atm, random, "gtd7d", 1.0e-50, 3.0e-14);
548 }
549
550 private void doTestDoubleMethod(NRLMSISE00 atm, RandomGenerator random, String methodName,
551 double absTolerance, double relTolerance)
552 {
553 try {
554
555 final int doy = 172;
556 final double sec = 29000.;
557 final double lat = 60.;
558 final double lon = -70.;
559 final double hl = 16.;
560 final double f107a = 149.;
561 final double f107 = 150.;
562 double[] ap = {4., 100., 100., 100., 100., 100., 100.};
563
564 Method methodD = getOutputClass().getDeclaredMethod(methodName, double[].class);
565 methodD.setAccessible(true);
566
567 Method methodF = getFieldOutputClass().getDeclaredMethod(methodName, double[].class);
568 methodF.setAccessible(true);
569
570 double[] p = new double[150];
571
572 Object output = createOutput(atm, doy, sec, lat, lon, hl, f107a, f107, ap);
573 Object fieldOutput = createFieldOutput(Decimal64Field.getInstance(),
574 atm, doy, sec, lat, lon, hl, f107a, f107, ap);
575 double maxAbsoluteError = 0;
576 double maxRelativeError = 0;
577 for (int i = 0; i < 100; ++i) {
578 for (int k = 0; k < p.length; ++k) {
579 p[k] = random.nextDouble();
580 }
581 double resDouble = ((Double) methodD.invoke(output, p)).doubleValue();
582 double resField = ((CalculusFieldElement<?>) methodF.invoke(fieldOutput, p)).getReal();
583 maxAbsoluteError = FastMath.max(maxAbsoluteError, FastMath.abs(resDouble - resField));
584 maxRelativeError = FastMath.max(maxRelativeError, FastMath.abs((resDouble - resField) / resDouble));
585 }
586 Assert.assertEquals(0.0, maxAbsoluteError, absTolerance);
587 if (maxAbsoluteError != 0.0) {
588 Assert.assertEquals(0.0, maxRelativeError, relTolerance);
589 }
590
591 } catch (NoSuchMethodException | SecurityException |
592 IllegalAccessException | IllegalArgumentException | InvocationTargetException e) {
593 Assert.fail(e.getLocalizedMessage());
594 }
595 }
596
597 private void doTestVoidMethod(NRLMSISE00 atm, RandomGenerator random, String methodName,
598 double temperatureRelativeTolerance, double densityRelativeTolerance)
599 {
600 try {
601
602
603 final int doy = 172;
604 final double sec = 29000.;
605 final double lat = 60.;
606 final double lon = -70.;
607 final double hl = 16.;
608 final double f107a = 149.;
609 final double f107 = 150.;
610 double[] ap = {4., 100., 100., 100., 100., 100., 100.};
611
612 Method methodD = getOutputClass().getDeclaredMethod(methodName, Double.TYPE);
613 methodD.setAccessible(true);
614
615 Method methodF = getFieldOutputClass().getDeclaredMethod(methodName, CalculusFieldElement.class);
616 methodF.setAccessible(true);
617
618 Object output = createOutput(atm, doy, sec, lat, lon, hl, f107a, f107, ap);
619 Object fieldOutput = createFieldOutput(Decimal64Field.getInstance(),
620 atm, doy, sec, lat, lon, hl, f107a, f107, ap);
621 double maxTemperatureError = 0;
622 double maxDensityError = 0;
623 for (int i = 0; i < 100; ++i) {
624 double alt = 500.0 * random.nextDouble();
625 methodD.invoke(output, alt);
626 methodF.invoke(fieldOutput, new Decimal64(alt));
627 for (int index = 0; index < 2; ++index) {
628 double tD = getOutputTemperature(output, index);
629 double tF = getFieldOutputTemperature(fieldOutput, index);
630 maxTemperatureError = FastMath.max(maxTemperatureError, FastMath.abs((tD - tF) / tF));
631 }
632 for (int index = 0; index < 9; ++index) {
633 double dD = getOutputDensity(output, index);
634 double dF = getFieldOutputDensity(fieldOutput, index);
635 if (Double.isNaN(dD)) {
636
637
638
639 Assert.assertTrue(Double.isNaN(dF));
640 } else if (dD == 0) {
641
642 Assert.assertEquals(dD, dF, Precision.SAFE_MIN);
643 } else {
644 maxDensityError = FastMath.max(maxDensityError, FastMath.abs((dD - dF) / dD));
645 }
646 }
647 }
648 Assert.assertEquals(0.0, maxTemperatureError, temperatureRelativeTolerance);
649 Assert.assertEquals(0.0, maxDensityError, densityRelativeTolerance);
650
651 } catch (NoSuchMethodException | SecurityException |
652 IllegalAccessException | IllegalArgumentException | InvocationTargetException e) {
653 Assert.fail(e.getLocalizedMessage());
654 }
655 }
656
657 private Class<?> getOutputClass() {
658 for (final Class<?> c : NRLMSISE00.class.getDeclaredClasses()) {
659 if (c.getName().endsWith("$Output")) {
660 return c;
661 }
662 }
663 return null;
664 }
665
666 private Object createOutput(final NRLMSISE00 atm,
667 final int doy, final double sec,
668 final double lat, final double lon, final double hl,
669 final double f107a, final double f107, final double[] ap) {
670 try {
671 Class<?> outputClass = getOutputClass();
672 Constructor<?> cons = outputClass.getDeclaredConstructor(NRLMSISE00.class,
673 Integer.TYPE,
674 Double.TYPE,
675 Double.TYPE,
676 Double.TYPE,
677 Double.TYPE,
678 Double.TYPE,
679 Double.TYPE,
680 double[].class);
681 cons.setAccessible(true);
682
683 return cons.newInstance(atm, doy, sec, lat, lon, hl, f107a, f107, ap);
684 } catch (NoSuchMethodException | SecurityException | InstantiationException |
685 IllegalAccessException | IllegalArgumentException | InvocationTargetException e) {
686 Assert.fail(e.getLocalizedMessage());
687 return null;
688 }
689
690 }
691
692 private double getOutputDensity(Object o, int index) {
693 try {
694 Method getDensity = getOutputClass().
695 getDeclaredMethod("getDensity", Integer.TYPE);
696 getDensity.setAccessible(true);
697 return ((Double) getDensity.invoke(o, index)).doubleValue();
698 } catch (NoSuchMethodException | SecurityException |
699 IllegalAccessException | IllegalArgumentException |
700 InvocationTargetException e) {
701 Assert.fail(e.getLocalizedMessage());
702 return Double.NaN;
703 }
704 }
705
706 private double getOutputTemperature(Object o, int index) {
707 try {
708 java.lang.reflect.Field temperaturesField = getOutputClass().getDeclaredField("temperatures");
709 temperaturesField.setAccessible(true);
710 return ((double[]) temperaturesField.get(o))[index];
711 } catch (NoSuchFieldException | SecurityException |
712 IllegalAccessException | IllegalArgumentException e) {
713 Assert.fail(e.getLocalizedMessage());
714 return Double.NaN;
715 }
716
717 }
718
719 private Class<?> getFieldOutputClass() {
720 for (final Class<?> c : NRLMSISE00.class.getDeclaredClasses()) {
721 if (c.getName().endsWith("$FieldOutput")) {
722 return c;
723 }
724 }
725 return null;
726 }
727
728 private <T extends CalculusFieldElement<T>> Object createFieldOutput(Field<T>field,
729 final NRLMSISE00 atm,
730 final int doy, final double sec,
731 final double lat, final double lon,
732 final double hl,
733 final double f107a, final double f107,
734 final double[] ap) {
735 try {
736 Class<?> fieldOutputClass = getFieldOutputClass();
737 Constructor<?> cons = fieldOutputClass.getDeclaredConstructor(NRLMSISE00.class,
738 Integer.TYPE,
739 CalculusFieldElement.class,
740 CalculusFieldElement.class,
741 CalculusFieldElement.class,
742 CalculusFieldElement.class,
743 Double.TYPE,
744 Double.TYPE,
745 double[].class);
746 cons.setAccessible(true);
747 return cons.newInstance(atm, doy,
748 field.getZero().add(sec),
749 field.getZero().add(lat),
750 field.getZero().add(lon),
751 field.getZero().add(hl),
752 f107a, f107, ap);
753 } catch (NoSuchMethodException | SecurityException | InstantiationException |
754 IllegalAccessException | IllegalArgumentException | InvocationTargetException e) {
755 Assert.fail(e.getLocalizedMessage());
756 return null;
757 }
758
759 }
760
761 private double getFieldOutputDensity(Object o, int index) {
762 try {
763 Method getDensity = getFieldOutputClass().
764 getDeclaredMethod("getDensity", Integer.TYPE);
765 getDensity.setAccessible(true);
766 return ((CalculusFieldElement<?>) getDensity.invoke(o, index)).getReal();
767 } catch (NoSuchMethodException | SecurityException |
768 IllegalAccessException | IllegalArgumentException |
769 InvocationTargetException e) {
770 Assert.fail(e.getLocalizedMessage());
771 return Double.NaN;
772 }
773 }
774
775 private double getFieldOutputTemperature(Object o, int index) {
776 try {
777 java.lang.reflect.Field temperaturesField = getFieldOutputClass().getDeclaredField("temperatures");
778 temperaturesField.setAccessible(true);
779 return ((CalculusFieldElement[]) temperaturesField.get(o))[index].getReal();
780 } catch (NoSuchFieldException | SecurityException |
781 IllegalAccessException | IllegalArgumentException e) {
782 Assert.fail(e.getLocalizedMessage());
783 return Double.NaN;
784 }
785
786 }
787
788 private void checkLegacy(final int nb, final Object out, final boolean print)
789 throws IllegalAccessException, IllegalArgumentException, InvocationTargetException {
790 final double[] tInfRef = {1.250540E+03, 1.166754E+03, 1.239892E+03, 1.027318E+03,
791 1.212396E+03, 1.220146E+03, 1.116385E+03, 1.031247E+03,
792 1.306052E+03, 1.361868E+03, 1.027318E+03, 1.027318E+03,
793 1.027318E+03, 1.027318E+03, 1.027318E+03, 1.426412E+03,
794 1.027318E+03};
795 final double[] tAltRef = {1.241416E+03, 1.161710E+03, 1.239891E+03, 2.068878E+02,
796 1.208135E+03, 1.212712E+03, 1.112999E+03, 1.024848E+03,
797 1.293374E+03, 1.347389E+03, 2.814648E+02, 2.274180E+02,
798 2.374389E+02, 2.795551E+02, 2.190732E+02, 1.408608E+03,
799 1.934071E+02};
800 final double[] dHeRef = {6.665177E+05, 3.407293E+06, 1.123767E+05, 5.411554E+07,
801 1.851122E+06, 8.673095E+05, 5.776251E+05, 3.740304E+05,
802 6.748339E+05, 5.528601E+05, 1.375488E+14, 4.427443E+13,
803 2.127829E+12, 1.412184E+11, 1.254884E+10, 5.196477E+05,
804 4.260860E+07};
805 final double[] dORef = {1.138806E+08, 1.586333E+08, 6.934130E+04, 1.918893E+11,
806 1.476555E+08, 1.278862E+08, 6.979139E+07, 4.782720E+07,
807 1.245315E+08, 1.198041E+08, 0.000000E+00, 0.000000E+00,
808 0.000000E+00, 0.000000E+00, 0.000000E+00, 1.274494E+08,
809 1.241342E+11};
810 final double[] dN2Ref = {1.998211E+07, 1.391117E+07, 4.247105E+01, 6.115826E+12,
811 1.579356E+07, 1.822577E+07, 1.236814E+07, 5.240380E+06,
812 2.369010E+07, 3.495798E+07, 2.049687E+19, 6.597567E+18,
813 3.170791E+17, 2.104370E+16, 1.874533E+15, 4.850450E+07,
814 4.929562E+12};
815 final double[] dO2Ref = {4.022764E+05, 3.262560E+05, 1.322750E-01, 1.225201E+12,
816 2.633795E+05, 2.922214E+05, 2.492868E+05, 1.759875E+05,
817 4.911583E+05, 9.339618E+05, 5.498695E+18, 1.769929E+18,
818 8.506280E+16, 5.645392E+15, 4.923051E+14, 1.720838E+06,
819 1.048407E+12};
820 final double[] dARRef = {3.557465E+03, 1.559618E+03, 2.618848E-05, 6.023212E+10,
821 1.588781E+03, 2.402962E+03, 1.405739E+03, 5.501649E+02,
822 4.578781E+03, 1.096255E+04, 2.451733E+17, 7.891680E+16,
823 3.792741E+15, 2.517142E+14, 2.239685E+13, 2.354487E+04,
824 4.993465E+10};
825 final double[] dHRef = {3.475312E+04, 4.854208E+04, 2.016750E+04, 1.059880E+07,
826 5.816167E+04, 3.686389E+04, 5.291986E+04, 8.896776E+04,
827 3.244595E+04, 2.686428E+04, 0.000000E+00, 0.000000E+00,
828 0.000000E+00, 0.000000E+00, 0.000000E+00, 2.500078E+04,
829 8.831229E+06};
830 final double[] dNRef = {4.095913E+06, 4.380967E+06, 5.741256E+03, 2.615737E+05,
831 5.478984E+06, 3.897276E+06, 1.069814E+06, 1.979741E+06,
832 5.370833E+06, 4.889974E+06, 0.000000E+00, 0.000000E+00,
833 0.000000E+00, 0.000000E+00, 0.000000E+00, 6.279210E+06,
834 2.252516E+05};
835 final double[] dAnORef = {2.667273E+04, 6.956682E+03, 2.374394E+04, 2.819879E-42,
836 1.264446E+03, 2.667273E+04, 2.667273E+04, 9.121815E+03,
837 2.667273E+04, 2.805445E+04, 0.000000E+00, 0.000000E+00,
838 0.000000E+00, 0.000000E+00, 0.000000E+00, 2.667273E+04,
839 2.415246E-42};
840 final double[] rhoRef = {4.074714E-15, 5.001846E-15, 2.756772E-18, 3.584426E-10,
841 4.809630E-15, 4.355866E-15, 2.470651E-15, 1.571889E-15,
842 4.564420E-15, 4.974543E-15, 1.261066E-03, 4.059139E-04,
843 1.950822E-05, 1.294709E-06, 1.147668E-07, 5.881940E-15,
844 2.914304E-10};
845 final double deltaT = 1.e-2;
846 final double deltaD = 5.e-7;
847 final int id = nb - 1;
848 if (print) {
849 System.out.printf("Case #%d\n", nb);
850 System.out.printf("Tinf: %E %E\n", tInfRef[id], getOutputTemperature(out, 0));
851 System.out.printf("Talt: %E %E\n", tAltRef[id], getOutputTemperature(out, 1));
852 System.out.printf("He: %E %E\n", dHeRef[id], getOutputDensity(out, 0) * 1e-6);
853 System.out.printf("O: %E %E\n", dORef[id], getOutputDensity(out, 1) * 1e-6);
854 System.out.printf("N2: %E %E\n", dN2Ref[id], getOutputDensity(out, 2) * 1e-6);
855 System.out.printf("O2: %E %E\n", dO2Ref[id], getOutputDensity(out, 3) * 1e-6);
856 System.out.printf("Ar: %E %E\n", dARRef[id], getOutputDensity(out, 4) * 1e-6);
857 System.out.printf("H: %E %E\n", dHRef[id], getOutputDensity(out, 6) * 1e-6);
858 System.out.printf("N: %E %E\n", dNRef[id], getOutputDensity(out, 7) * 1e-6);
859 System.out.printf("AnO: %E %E\n", dAnORef[id], getOutputDensity(out, 8) * 1e-6);
860 System.out.printf("Rho: %E %E\n\n", rhoRef[id], getOutputDensity(out, 5) * 1e-3);
861 } else {
862 Assert.assertEquals(tInfRef[id], getOutputTemperature(out, 0), deltaT);
863 Assert.assertEquals(tAltRef[id], getOutputTemperature(out, 1), deltaT);
864 Assert.assertEquals(dHeRef[id], getOutputDensity(out, 0) * 1e-6, dHeRef[id] * deltaD);
865 Assert.assertEquals(dORef[id], getOutputDensity(out, 1) * 1e-6, dORef[id] * deltaD);
866 Assert.assertEquals(dN2Ref[id], getOutputDensity(out, 2) * 1e-6, dN2Ref[id] * deltaD);
867 Assert.assertEquals(dO2Ref[id], getOutputDensity(out, 3) * 1e-6, dO2Ref[id] * deltaD);
868 Assert.assertEquals(dARRef[id], getOutputDensity(out, 4) * 1e-6, dARRef[id] * deltaD);
869 Assert.assertEquals(dHRef[id], getOutputDensity(out, 6) * 1e-6, dHRef[id] * deltaD);
870 Assert.assertEquals(dNRef[id], getOutputDensity(out, 7) * 1e-6, dNRef[id] * deltaD);
871 Assert.assertEquals(dAnORef[id], getOutputDensity(out, 8) * 1e-6, dAnORef[id] * deltaD);
872 Assert.assertEquals(rhoRef[id], getOutputDensity(out, 5) * 1e-3, rhoRef[id] * deltaD);
873 }
874
875 }
876
877 @Before
878 public void setUp() {
879 Utils.setDataRoot("regular-data");
880 }
881
882 private static class InputParams implements NRLMSISE00InputParameters {
883
884
885 private static final long serialVersionUID = 1L;
886
887
888 public InputParams() {
889
890 }
891
892 @Override
893 public AbsoluteDate getMinDate() {
894 return new AbsoluteDate(2003, 1, 1, TimeScalesFactory.getUTC());
895 }
896
897 @Override
898 public AbsoluteDate getMaxDate() {
899 return new AbsoluteDate(2003, 12, 31, TimeScalesFactory.getUTC());
900 }
901
902 @Override
903 public double getDailyFlux(AbsoluteDate date) {
904 return 150.;
905 }
906
907 @Override
908 public double getAverageFlux(AbsoluteDate date) {
909 return 150.;
910 }
911
912 @Override
913 public double[] getAp(AbsoluteDate date) {
914 return new double[] {4., 100., 100., 100., 100., 100., 100.};
915 }
916 }
917 }