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