1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.models.earth.ionosphere.nequick;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.Field;
21 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
22 import org.hipparchus.geometry.euclidean.threed.Vector3D;
23 import org.hipparchus.util.Binary64;
24 import org.hipparchus.util.Binary64Field;
25 import org.hipparchus.util.FastMath;
26 import org.junit.jupiter.api.Assertions;
27 import org.junit.jupiter.api.BeforeEach;
28 import org.junit.jupiter.api.Test;
29 import org.orekit.Utils;
30 import org.orekit.bodies.FieldGeodeticPoint;
31 import org.orekit.bodies.GeodeticPoint;
32 import org.orekit.bodies.OneAxisEllipsoid;
33 import org.orekit.frames.FramesFactory;
34 import org.orekit.frames.TopocentricFrame;
35 import org.orekit.gnss.PredefinedGnssSignal;
36 import org.orekit.orbits.CartesianOrbit;
37 import org.orekit.orbits.FieldCartesianOrbit;
38 import org.orekit.orbits.FieldOrbit;
39 import org.orekit.orbits.Orbit;
40 import org.orekit.propagation.FieldSpacecraftState;
41 import org.orekit.propagation.SpacecraftState;
42 import org.orekit.time.AbsoluteDate;
43 import org.orekit.time.FieldAbsoluteDate;
44 import org.orekit.time.TimeScalesFactory;
45 import org.orekit.utils.Constants;
46 import org.orekit.utils.FieldPVCoordinates;
47 import org.orekit.utils.IERSConventions;
48 import org.orekit.utils.PVCoordinates;
49
50
51
52
53
54 public class NeQuickItuTest {
55
56 @BeforeEach
57 public void setUp() {
58 Utils.setDataRoot("regular-data");
59 }
60
61 @Test
62 public void testMediumSolarActivityItu() {
63
64
65 final NeQuickItu model = new NeQuickItu(137.568737, TimeScalesFactory.getUTC());
66
67
68 Assertions.assertEquals(137.568737, model.getF107(), 1.0e-6);
69 Assertions.assertEquals("UTC", model.getUtc().getName());
70
71
72 final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(-31.80), FastMath.toRadians(115.89), 12.78);
73 final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(-14.31), FastMath.toRadians(124.09), 20100697.90);
74
75
76 final AbsoluteDate date = new AbsoluteDate(2018, 1, 2, 16, 0, 0, TimeScalesFactory.getUTC());
77
78
79 final double stec = model.stec(date, recP, satP);
80 Assertions.assertEquals(14.607, stec, 1.0e-3);
81 }
82
83 @Test
84 public void testFieldMediumSolarActivityItu() {
85 doTestFieldMediumSolarActivityItu(Binary64Field.getInstance());
86 }
87
88 private <T extends CalculusFieldElement<T>> void doTestFieldMediumSolarActivityItu(final Field<T> field) {
89
90
91 final T zero = field.getZero();
92
93
94 final NeQuickItu model = new NeQuickItu(137.568737, TimeScalesFactory.getUTC());
95
96
97 final FieldGeodeticPoint<T> recP = new FieldGeodeticPoint<>(zero.newInstance(FastMath.toRadians(-31.80)),
98 zero.newInstance(FastMath.toRadians(115.89)),
99 zero.newInstance(12.78));
100 final FieldGeodeticPoint<T> satP = new FieldGeodeticPoint<>(zero.newInstance(FastMath.toRadians(-14.31)),
101 zero.newInstance(FastMath.toRadians(124.09)),
102 zero.newInstance(20100697.90));
103
104
105 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2018, 1, 2, 16, 0, 0, TimeScalesFactory.getUTC());
106
107
108 final T stec = model.stec(date, recP, satP);
109 Assertions.assertEquals(14.607, stec.getReal(), 1.0e-3);
110 }
111
112
113 @Test
114 public void testValidationItu63() {
115 doTestValidationItu(63, 4.36189);
116 }
117
118
119 @Test
120 public void testValidationItu128() {
121 doTestValidationItu(128, 11.83090);
122 }
123
124
125 @Test
126 public void testValidationItu193() {
127 doTestValidationItu(193, 20.85232);
128 }
129
130 private void doTestValidationItu(final double flux, final double expected) {
131
132
133 final NeQuickItu model = new NeQuickItu(flux, TimeScalesFactory.getUTC());
134
135
136 final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(82.494293510492980204),
137 FastMath.toRadians(-62.340460202287275138),
138 78.107446296427042398);
139 final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(54.445029415916160076),
140 FastMath.toRadians(-118.47006897550868132),
141 20370730.845002099872);
142
143
144 final AbsoluteDate date = new AbsoluteDate(2007, 4, 1, TimeScalesFactory.getUTC());
145
146
147 final double stec = model.stec(date, recP, satP);
148 Assertions.assertEquals(expected, stec, 1.0e-3);
149 }
150
151
152 @Test
153 public void testValidationItu63Field() {
154 doTestValidationItu(Binary64Field.getInstance(), 63, 4.36189);
155 }
156
157
158 @Test
159 public void testValidationItu128Field() {
160 doTestValidationItu(Binary64Field.getInstance(), 128, 11.83090);
161 }
162
163
164 @Test
165 public void testValidationItu193Field() {
166 doTestValidationItu(Binary64Field.getInstance(), 193, 20.85232);
167 }
168
169 private <T extends CalculusFieldElement<T>> void doTestValidationItu(final Field<T> field,
170 final double flux, final double expected) {
171
172
173 final NeQuickItu model = new NeQuickItu(flux, TimeScalesFactory.getUTC());
174
175
176 final FieldGeodeticPoint<T> recP =
177 new FieldGeodeticPoint<>(field,
178 new GeodeticPoint(FastMath.toRadians(82.494293510492980204),
179 FastMath.toRadians(-62.340460202287275138),
180 78.1074462964270423981));
181 final FieldGeodeticPoint<T> satP =
182 new FieldGeodeticPoint<>(field,
183 new GeodeticPoint(FastMath.toRadians(54.4450294159161600763),
184 FastMath.toRadians(-118.47006897550868132),
185 20370730.8450020998725));
186
187
188 final FieldAbsoluteDate<T> date =
189 new FieldAbsoluteDate<>(field, new AbsoluteDate(2007, 4, 1, TimeScalesFactory.getUTC()));
190
191
192 final T stec = model.stec(date, recP, satP);
193 Assertions.assertEquals(expected, stec.getReal(), 1.0e-3);
194 }
195
196 @Test
197 public void testHeights1() {
198 doTestHeights(30.0, 205000000.0, 10.730);
199 }
200
201 @Test
202 public void testHeights2() {
203 doTestHeights(1300000.0, 205000000.0, 0.943);
204 }
205
206 @Test
207 public void testHeights3() {
208 doTestHeights(2300000.0, 205000000.0, 0.580);
209 }
210
211 @Test
212 public void testHeights4() {
213 doTestHeights(30.0, 1950000.0, 29.670);
214 }
215
216 @Test
217 public void testHeights5() {
218 doTestHeights(1030000.0, 1950000.0, 2.594);
219 }
220
221 private void doTestHeights(final double hRec, final double hSat, final double expected) {
222
223
224 final NeQuickItu model = new NeQuickItu(128.0, TimeScalesFactory.getUTC());
225
226
227 final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians( 82.494),
228 FastMath.toRadians(-62.340),
229 hRec);
230 final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(54.445),
231 FastMath.toRadians(-118.470),
232 hSat);
233
234
235 final AbsoluteDate date = new AbsoluteDate(2007, 4, 1, TimeScalesFactory.getUTC());
236
237
238 final double stec = model.stec(date, recP, satP);
239 Assertions.assertEquals(expected, stec, 1.0e-3);
240
241 }
242
243 @Test
244 public void testMeridian() {
245 final NeQuickItu model = new NeQuickItu(128.0, TimeScalesFactory.getUTC());
246 final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(1.0e-3), FastMath.toRadians(0), 0);
247 final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(0.761e-3), FastMath.toRadians(0), 2.0e6);
248 final AbsoluteDate date = new AbsoluteDate(2007, 4, 1, TimeScalesFactory.getUTC());
249 Assertions.assertEquals(38.031, model.stec(date, recP, satP), 1.0e-3);
250 }
251
252 @Test
253 public void testMeridianField() {
254 final Field<Binary64> field = Binary64Field.getInstance();
255 final Binary64 zero = field.getZero();
256 final NeQuickItu model = new NeQuickItu(128.0, TimeScalesFactory.getUTC());
257 final FieldGeodeticPoint<Binary64> recP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(1.0e-3)),
258 FastMath.toRadians(zero.newInstance(0)),
259 zero.newInstance(0));
260 final FieldGeodeticPoint<Binary64> satP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(0.761e-3)),
261 FastMath.toRadians(zero.newInstance(0)),
262 zero.newInstance(2.0e6));
263 final FieldAbsoluteDate<Binary64> date = new FieldAbsoluteDate<>(field,
264 new AbsoluteDate(2007, 4, 1, TimeScalesFactory.getUTC()));
265 Assertions.assertEquals(38.031, model.stec(date, recP, satP).getReal(), 1.0e-3);
266 }
267
268 @Test
269 public void testHeights1Field() {
270 doTestHeights(Binary64Field.getInstance(), 30.0, 205000000.0, 10.730);
271 }
272
273 @Test
274 public void testHeights2Field() {
275 doTestHeights(Binary64Field.getInstance(), 1300000.0, 205000000.0, 0.943);
276 }
277
278 @Test
279 public void testHeights3Field() {
280 doTestHeights(Binary64Field.getInstance(), 2300000.0, 205000000.0, 0.580);
281 }
282
283 @Test
284 public void testHeights4Field() {
285 doTestHeights(Binary64Field.getInstance(), 30.0, 1950000.0, 29.670);
286 }
287
288 @Test
289 public void testHeights5Field() {
290 doTestHeights(Binary64Field.getInstance(), 1030000.0, 1950000.0, 2.594);
291 }
292
293 private <T extends CalculusFieldElement<T>> void doTestHeights(final Field<T> field,
294 final double hRec, final double hSat, final double expected) {
295
296 final T zero = field.getZero();
297
298
299 final NeQuickItu model = new NeQuickItu(128.0, TimeScalesFactory.getUTC());
300
301
302 final FieldGeodeticPoint<T> recP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance( 82.494)),
303 FastMath.toRadians(zero.newInstance(-62.340)),
304 zero.newInstance(hRec));
305 final FieldGeodeticPoint<T> satP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(54.445)),
306 FastMath.toRadians(zero.newInstance(-118.470)),
307 zero.newInstance(hSat));
308
309
310 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2007, 4, 1, TimeScalesFactory.getUTC());
311
312
313 final T stec = model.stec(date, recP, satP);
314 Assertions.assertEquals(expected, stec.getReal(), 1.0e-3);
315
316 }
317
318 @Test
319 public void testDelay() {
320
321
322 final NeQuickItu model = new NeQuickItu(128.0, TimeScalesFactory.getUTC());
323
324
325 final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(-31.80), FastMath.toRadians(115.89), 12.78);
326 final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(-14.31), FastMath.toRadians(124.09), 20100697.90);
327
328
329 final AbsoluteDate date = new AbsoluteDate(2018, 4, 2, 16, 0, 0, TimeScalesFactory.getUTC());
330
331
332 final OneAxisEllipsoid ellipsoid = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
333 Constants.WGS84_EARTH_FLATTENING,
334 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
335
336 final Vector3D satPosInITRF = ellipsoid.transform(satP);
337 final Vector3D satPosInEME2000 = ellipsoid.getBodyFrame().getStaticTransformTo(FramesFactory.getEME2000(), date).transformPosition(satPosInITRF);
338
339
340 final PVCoordinates pv = new PVCoordinates(satPosInEME2000, new Vector3D(1.0, 1.0, 1.0));
341 final Orbit orbit = new CartesianOrbit(pv, FramesFactory.getEME2000(), date, Constants.WGS84_EARTH_MU);
342 final SpacecraftState state = new SpacecraftState(orbit);
343
344 final double delay = model.pathDelay(state, new TopocentricFrame(ellipsoid, recP, null),
345 PredefinedGnssSignal.G01.getFrequency(), model.getParameters());
346
347
348 Assertions.assertEquals(1.32, delay, 0.01);
349 }
350
351 @Test
352 public void testFieldDelay() {
353 doTestFieldDelay(Binary64Field.getInstance());
354 }
355
356 private <T extends CalculusFieldElement<T>> void doTestFieldDelay(final Field<T> field) {
357
358
359 final T zero = field.getZero();
360 final T one = field.getOne();
361
362
363 final NeQuickItu model = new NeQuickItu(128.0, TimeScalesFactory.getUTC());
364
365
366 final double recLat = FastMath.toRadians(-31.80);
367 final double recLon = FastMath.toRadians(115.89);
368 final double recAlt = 12.78;
369 final GeodeticPoint recP = new GeodeticPoint(recLat, recLon, recAlt);
370 final FieldGeodeticPoint<T> satP = new FieldGeodeticPoint<>(zero.newInstance(FastMath.toRadians(-14.31)),
371 zero.newInstance(FastMath.toRadians(124.09)), zero.newInstance(20100697.90));
372
373
374 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2018, 4, 2, 16, 0, 0, TimeScalesFactory.getUTC());
375
376
377 final OneAxisEllipsoid ellipsoid = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
378 Constants.WGS84_EARTH_FLATTENING,
379 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
380
381 final FieldVector3D<T> satPosInITRF = ellipsoid.transform(satP);
382 final FieldVector3D<T> satPosInEME2000 = ellipsoid.getBodyFrame().getStaticTransformTo(FramesFactory.getEME2000(), date).transformPosition(satPosInITRF);
383
384
385 final FieldPVCoordinates<T> pv = new FieldPVCoordinates<>(satPosInEME2000, new FieldVector3D<>(one, one, one));
386 final FieldOrbit<T> orbit = new FieldCartesianOrbit<>(pv, FramesFactory.getEME2000(), date, zero.newInstance(Constants.WGS84_EARTH_MU));
387 final FieldSpacecraftState<T> state = new FieldSpacecraftState<>(orbit);
388
389 final T delay = model.pathDelay(state, new TopocentricFrame(ellipsoid, recP, null),
390 PredefinedGnssSignal.G01.getFrequency(), model.getParameters(field));
391
392
393 Assertions.assertEquals(1.32, delay.getReal(), 0.01);
394 }
395
396 @Test
397 public void testAntiMeridian() {
398
399
400 final NeQuickItu model = new NeQuickItu(128.0, TimeScalesFactory.getUTC());
401
402
403 final AbsoluteDate date = new AbsoluteDate(2018, 11, 2, 16, 0, 0, TimeScalesFactory.getUTC());
404
405
406 final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(-31.80), FastMath.toRadians(-179.99), 12.78);
407 final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(-14.31), FastMath.toRadians(-177.43), 20100697.90);
408 double stec = model.stec(date, recP, satP);
409 Assertions.assertEquals(13.269, stec, 0.001);
410
411 }
412
413 @Test
414 public void testFieldAntiMeridian() {
415 doTestFieldAntiMeridian(Binary64Field.getInstance());
416 }
417
418 private <T extends CalculusFieldElement<T>> void doTestFieldAntiMeridian(final Field<T> field) {
419
420 final T zero = field.getZero();
421
422
423 final NeQuickItu model = new NeQuickItu(128.0, TimeScalesFactory.getUTC());
424
425
426 final FieldAbsoluteDate<T> date =
427 new FieldAbsoluteDate<>(field,
428 new AbsoluteDate(2018, 11, 2, 16, 0, 0, TimeScalesFactory.getUTC()));
429
430
431 final FieldGeodeticPoint<T> recP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(-31.80)),
432 FastMath.toRadians(zero.newInstance(-179.99)),
433 zero.newInstance(12.78));
434 final FieldGeodeticPoint<T> satP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(-14.31)),
435 FastMath.toRadians(zero.newInstance(-177.43)),
436 zero.newInstance(20100697.90));
437 T stec = model.stec(date, recP, satP);
438 Assertions.assertEquals(13.269, stec.getReal(), 0.001);
439
440 }
441
442 @Test
443 public void testPole() {
444
445
446 final NeQuickItu model = new NeQuickItu(128.0, TimeScalesFactory.getUTC());
447
448
449 final AbsoluteDate date = new AbsoluteDate(2018, 11, 2, 16, 0, 0, TimeScalesFactory.getUTC());
450
451
452 final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(90), FastMath.toRadians(0), 100);
453 final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(-14.31), FastMath.toRadians(-177.43), 20100697.90);
454 double stec = model.stec(date, recP, satP);
455 Assertions.assertEquals(70.372, stec, 0.001);
456
457 }
458
459 @Test
460 public void testFieldPole() {
461 doTestFieldPole(Binary64Field.getInstance());
462 }
463
464 private <T extends CalculusFieldElement<T>> void doTestFieldPole(final Field<T> field) {
465
466 final T zero = field.getZero();
467
468
469 final NeQuickItu model = new NeQuickItu(128.0, TimeScalesFactory.getUTC());
470
471
472 final FieldAbsoluteDate<T> date =
473 new FieldAbsoluteDate<>(field,
474 new AbsoluteDate(2018, 11, 2, 16, 0, 0, TimeScalesFactory.getUTC()));
475
476
477 final FieldGeodeticPoint<T> recP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(90)),
478 FastMath.toRadians(zero.newInstance(0)),
479 zero.newInstance(12.78));
480 final FieldGeodeticPoint<T> satP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(-14.31)),
481 FastMath.toRadians(zero.newInstance(-177.43)),
482 zero.newInstance(20100697.90));
483 T stec = model.stec(date, recP, satP);
484 Assertions.assertEquals(70.373, stec.getReal(), 0.001);
485
486 }
487
488 @Test
489 public void testZenith() {
490
491
492 final NeQuickItu model = new NeQuickItu(128.0, TimeScalesFactory.getUTC());
493
494
495 final AbsoluteDate date = new AbsoluteDate(2018, 4, 2, 16, 0, 0, TimeScalesFactory.getUTC());
496
497
498 final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(51.678), FastMath.toRadians(-9.448), 0.0);
499 final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(51.678), FastMath.toRadians(-9.448), 600000.0);
500 double stec = model.stec(date, recP, satP);
501 Assertions.assertEquals(16.618, stec, 0.001);
502
503 }
504
505 @Test
506 public void testFieldZenith() {
507 doTestFieldZenith(Binary64Field.getInstance());
508 }
509
510 private <T extends CalculusFieldElement<T>> void doTestFieldZenith(final Field<T> field) {
511
512 final T zero = field.getZero();
513
514
515 final NeQuickItu model = new NeQuickItu(128.0, TimeScalesFactory.getUTC());
516
517
518 final FieldAbsoluteDate<T> date =
519 new FieldAbsoluteDate<>(field,
520 new AbsoluteDate(2018, 4, 2, 16, 0, 0, TimeScalesFactory.getUTC()));
521
522
523 final FieldGeodeticPoint<T> recP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(51.678)),
524 FastMath.toRadians(zero.newInstance(-9.448)),
525 zero.newInstance(0.0));
526 final FieldGeodeticPoint<T> satP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(51.678)),
527 FastMath.toRadians(zero.newInstance(-9.448)),
528 zero.newInstance(600000.0));
529 T stec = model.stec(date, recP, satP);
530 Assertions.assertEquals(16.618, stec.getReal(), 0.001);
531
532 }
533
534 }