1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.models.earth.troposphere;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.Field;
21 import org.hipparchus.util.FastMath;
22 import org.hipparchus.util.FieldSinCos;
23 import org.hipparchus.util.MathArrays;
24 import org.hipparchus.util.SinCos;
25 import org.orekit.annotation.DefaultDataContext;
26 import org.orekit.bodies.FieldGeodeticPoint;
27 import org.orekit.bodies.GeodeticPoint;
28 import org.orekit.data.DataContext;
29 import org.orekit.models.earth.weather.FieldPressureTemperatureHumidity;
30 import org.orekit.models.earth.weather.PressureTemperatureHumidity;
31 import org.orekit.time.AbsoluteDate;
32 import org.orekit.time.DateTimeComponents;
33 import org.orekit.time.FieldAbsoluteDate;
34 import org.orekit.time.TimeScale;
35 import org.orekit.utils.FieldLegendrePolynomials;
36 import org.orekit.utils.FieldTrackingCoordinates;
37 import org.orekit.utils.LegendrePolynomials;
38 import org.orekit.utils.TrackingCoordinates;
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61 @SuppressWarnings("deprecation")
62 public class GlobalMappingFunctionModel implements MappingFunction, TroposphereMappingFunction {
63
64
65 private static final double FACTOR = 1.0e-5;
66
67
68 private final TimeScale utc;
69
70
71
72
73
74
75
76 @DefaultDataContext
77 public GlobalMappingFunctionModel() {
78 this(DataContext.getDefault().getTimeScales().getUTC());
79 }
80
81
82
83
84
85 public GlobalMappingFunctionModel(final TimeScale utc) {
86 this.utc = utc;
87 }
88
89
90 @Override
91 @Deprecated
92 public double[] mappingFactors(final double elevation, final GeodeticPoint point,
93 final AbsoluteDate date) {
94 return mappingFactors(new TrackingCoordinates(0.0, elevation, 0.0), point,
95 TroposphericModelUtils.STANDARD_ATMOSPHERE,
96 date);
97 }
98
99
100 @Override
101 public double[] mappingFactors(final TrackingCoordinates trackingCoordinates,
102 final GeodeticPoint point,
103 final PressureTemperatureHumidity weather,
104 final AbsoluteDate date) {
105
106 final DateTimeComponents dtc = date.getComponents(utc);
107 final int dofyear = dtc.getDate().getDayOfYear();
108
109
110 final double bh = 0.0029;
111 final double c0h = 0.062;
112 final double c10h;
113 final double c11h;
114 final double psi;
115
116
117 final double latitude = point.getLatitude();
118 final double longitude = point.getLongitude();
119
120 if (FastMath.sin(latitude) > 0) {
121
122 c10h = 0.001;
123 c11h = 0.005;
124 psi = 0.0;
125 } else {
126
127 c10h = 0.002;
128 c11h = 0.007;
129 psi = FastMath.PI;
130 }
131
132 double t0 = 28;
133 if (latitude < 0) {
134
135 t0 += 183;
136 }
137 final double coef = ((dofyear + 1 - t0) / 365.25) * 2 * FastMath.PI + psi;
138 final double ch = c0h + ((FastMath.cos(coef) + 1) * (c11h / 2.0) + c10h) * (1.0 - FastMath.cos(latitude));
139
140
141 final double bw = 0.00146;
142 final double cw = 0.04391;
143
144
145
146
147 final int degree = 9;
148 final int order = 9;
149 final LegendrePolynomials p = new LegendrePolynomials(degree, order, FastMath.sin(latitude));
150
151 double a0Hydro = 0.;
152 double amplHydro = 0.;
153 double a0Wet = 0.;
154 double amplWet = 0.;
155 final ABCoefficients abCoef = new ABCoefficients();
156 int j = 0;
157 for (int n = 0; n <= 9; n++) {
158 for (int m = 0; m <= n; m++) {
159
160 final SinCos sc = FastMath.sinCos(m * longitude);
161
162 a0Hydro = a0Hydro + (abCoef.getAHMean(j) * p.getPnm(n, m) * sc.cos() +
163 abCoef.getBHMean(j) * p.getPnm(n, m) * sc.sin()) * FACTOR;
164
165 a0Wet = a0Wet + (abCoef.getAWMean(j) * p.getPnm(n, m) * sc.cos() +
166 abCoef.getBWMean(j) * p.getPnm(n, m) * sc.sin()) * FACTOR;
167
168 amplHydro = amplHydro + (abCoef.getAHAmplitude(j) * p.getPnm(n, m) * sc.cos() +
169 abCoef.getBHAmplitude(j) * p.getPnm(n, m) * sc.sin()) * FACTOR;
170
171 amplWet = amplWet + (abCoef.getAWAmplitude(j) * p.getPnm(n, m) * sc.cos() +
172 abCoef.getBWAmplitude(j) * p.getPnm(n, m) * sc.sin()) * FACTOR;
173
174 j = j + 1;
175 }
176 }
177
178
179 final double ah = a0Hydro + amplHydro * FastMath.cos(coef - psi);
180 final double aw = a0Wet + amplWet * FastMath.cos(coef - psi);
181
182 final double[] function = new double[2];
183 function[0] = TroposphericModelUtils.mappingFunction(ah, bh, ch,
184 trackingCoordinates.getElevation());
185 function[1] = TroposphericModelUtils.mappingFunction(aw, bw, cw,
186 trackingCoordinates.getElevation());
187
188
189 final double correction = TroposphericModelUtils.computeHeightCorrection(trackingCoordinates.getElevation(),
190 point.getAltitude());
191 function[0] = function[0] + correction;
192
193 return function;
194 }
195
196
197 @Override
198 @Deprecated
199 public <T extends CalculusFieldElement<T>> T[] mappingFactors(final T elevation, final FieldGeodeticPoint<T> point,
200 final FieldAbsoluteDate<T> date) {
201 return mappingFactors(new FieldTrackingCoordinates<>(date.getField().getZero(), elevation, date.getField().getZero()),
202 point,
203 new FieldPressureTemperatureHumidity<>(date.getField(),
204 TroposphericModelUtils.STANDARD_ATMOSPHERE),
205 date);
206 }
207
208
209 @Override
210 public <T extends CalculusFieldElement<T>> T[] mappingFactors(final FieldTrackingCoordinates<T> trackingCoordinates,
211 final FieldGeodeticPoint<T> point,
212 final FieldPressureTemperatureHumidity<T> weather,
213 final FieldAbsoluteDate<T> date) {
214
215 final DateTimeComponents dtc = date.getComponents(utc);
216 final int dofyear = dtc.getDate().getDayOfYear();
217
218 final Field<T> field = date.getField();
219 final T zero = field.getZero();
220
221
222 final T bh = zero.newInstance(0.0029);
223 final T c0h = zero.newInstance(0.062);
224 final T c10h;
225 final T c11h;
226 final T psi;
227
228
229 final T latitude = point.getLatitude();
230 final T longitude = point.getLongitude();
231
232
233 if (FastMath.sin(latitude.getReal()) > 0) {
234 c10h = zero.newInstance(0.001);
235 c11h = zero.newInstance(0.005);
236 psi = zero;
237 } else {
238 c10h = zero.newInstance(0.002);
239 c11h = zero.newInstance(0.007);
240 psi = zero.getPi();
241 }
242
243 double t0 = 28;
244 if (latitude.getReal() < 0) {
245
246 t0 += 183;
247 }
248 final T coef = psi.add(zero.getPi().multiply(2.0).multiply((dofyear + 1 - t0) / 365.25));
249 final T ch = c11h.divide(2.0).multiply(FastMath.cos(coef).add(1.0)).add(c10h).multiply(FastMath.cos(latitude).negate().add(1.0)).add(c0h);
250
251
252 final T bw = zero.newInstance(0.00146);
253 final T cw = zero.newInstance(0.04391);
254
255
256
257
258 final int degree = 9;
259 final int order = 9;
260 final FieldLegendrePolynomials<T> p = new FieldLegendrePolynomials<>(degree, order, FastMath.sin(latitude));
261
262 T a0Hydro = zero;
263 T amplHydro = zero;
264 T a0Wet = zero;
265 T amplWet = zero;
266 final ABCoefficients abCoef = new ABCoefficients();
267 int j = 0;
268 for (int n = 0; n <= 9; n++) {
269 for (int m = 0; m <= n; m++) {
270
271 final FieldSinCos<T> sc = FastMath.sinCos(longitude.multiply(m));
272
273
274 a0Hydro = a0Hydro.add(p.getPnm(n, m).multiply(abCoef.getAHMean(j)).multiply(sc.cos()).
275 add(p.getPnm(n, m).multiply(abCoef.getBHMean(j)).multiply(sc.sin())).
276 multiply(FACTOR));
277
278 a0Wet = a0Wet.add(p.getPnm(n, m).multiply(abCoef.getAWMean(j)).multiply(sc.cos()).
279 add(p.getPnm(n, m).multiply(abCoef.getBWMean(j)).multiply(sc.sin())).
280 multiply(FACTOR));
281
282 amplHydro = amplHydro.add(p.getPnm(n, m).multiply(abCoef.getAHAmplitude(j)).multiply(sc.cos()).
283 add(p.getPnm(n, m).multiply(abCoef.getBHAmplitude(j)).multiply(sc.sin())).
284 multiply(FACTOR));
285
286 amplWet = amplWet.add(p.getPnm(n, m).multiply(abCoef.getAWAmplitude(j)).multiply(sc.cos()).
287 add(p.getPnm(n, m).multiply(abCoef.getBWAmplitude(j)).multiply(sc.sin())).
288 multiply(FACTOR));
289
290 j = j + 1;
291 }
292 }
293
294
295 final T ah = a0Hydro.add(amplHydro.multiply(FastMath.cos(coef.subtract(psi))));
296 final T aw = a0Wet.add(amplWet.multiply(FastMath.cos(coef.subtract(psi))));
297
298 final T[] function = MathArrays.buildArray(field, 2);
299 function[0] = TroposphericModelUtils.mappingFunction(ah, bh, ch,
300 trackingCoordinates.getElevation());
301 function[1] = TroposphericModelUtils.mappingFunction(aw, bw, cw,
302 trackingCoordinates.getElevation());
303
304
305 final T correction = TroposphericModelUtils.computeHeightCorrection(trackingCoordinates.getElevation(),
306 point.getAltitude(),
307 field);
308 function[0] = function[0].add(correction);
309
310 return function;
311 }
312
313 private static class ABCoefficients {
314
315
316 private static final double[] AH_MEAN = {
317 1.2517e02,
318 8.503e-01,
319 6.936e-02,
320 -6.760e+00,
321 1.771e-01,
322 1.130e-02,
323 5.963e-01,
324 1.808e-02,
325 2.801e-03,
326 -1.414e-03,
327 -1.212e+00,
328 9.300e-02,
329 3.683e-03,
330 1.095e-03,
331 4.671e-05,
332 3.959e-01,
333 -3.867e-02,
334 5.413e-03,
335 -5.289e-04,
336 3.229e-04,
337 2.067e-05,
338 3.000e-01,
339 2.031e-02,
340 5.900e-03,
341 4.573e-04,
342 -7.619e-05,
343 2.327e-06,
344 3.845e-06,
345 1.182e-01,
346 1.158e-02,
347 5.445e-03,
348 6.219e-05,
349 4.204e-06,
350 -2.093e-06,
351 1.540e-07,
352 -4.280e-08,
353 -4.751e-01,
354 -3.490e-02,
355 1.758e-03,
356 4.019e-04,
357 -2.799e-06,
358 -1.287e-06,
359 5.468e-07,
360 7.580e-08,
361 -6.300e-09,
362 -1.160e-01,
363 8.301e-03,
364 8.771e-04,
365 9.955e-05,
366 -1.718e-06,
367 -2.012e-06,
368 1.170e-08,
369 1.790e-08,
370 -1.300e-09,
371 1.000e-10
372 };
373
374
375 private static final double[] BH_MEAN = {
376 0.000e+00,
377 0.000e+00,
378 3.249e-02,
379 0.000e+00,
380 3.324e-02,
381 1.850e-02,
382 0.000e+00,
383 -1.115e-01,
384 2.519e-02,
385 4.923e-03,
386 0.000e+00,
387 2.737e-02,
388 1.595e-02,
389 -7.332e-04,
390 1.933e-04,
391 0.000e+00,
392 -4.796e-02,
393 6.381e-03,
394 -1.599e-04,
395 -3.685e-04,
396 1.815e-05,
397 0.000e+00,
398 7.033e-02,
399 2.426e-03,
400 -1.111e-03,
401 -1.357e-04,
402 -7.828e-06,
403 2.547e-06,
404 0.000e+00,
405 5.779e-03,
406 3.133e-03,
407 -5.312e-04,
408 -2.028e-05,
409 2.323e-07,
410 -9.100e-08,
411 -1.650e-08,
412 0.000e+00,
413 3.688e-02,
414 -8.638e-04,
415 -8.514e-05,
416 -2.828e-05,
417 5.403e-07,
418 4.390e-07,
419 1.350e-08,
420 1.800e-09,
421 0.000e+00,
422 -2.736e-02,
423 -2.977e-04,
424 8.113e-05,
425 2.329e-07,
426 8.451e-07,
427 4.490e-08,
428 -8.100e-09,
429 -1.500e-09,
430 2.000e-10
431 };
432
433
434 private static final double[] AH_AMPL = {
435 -2.738e-01,
436 -2.837e+00,
437 1.298e-02,
438 -3.588e-01,
439 2.413e-02,
440 3.427e-02,
441 -7.624e-01,
442 7.272e-02,
443 2.160e-02,
444 -3.385e-03,
445 4.424e-01,
446 3.722e-02,
447 2.195e-02,
448 -1.503e-03,
449 2.426e-04,
450 3.013e-01,
451 5.762e-02,
452 1.019e-02,
453 -4.476e-04,
454 6.790e-05,
455 3.227e-05,
456 3.123e-01,
457 -3.535e-02,
458 4.840e-03,
459 3.025e-06,
460 -4.363e-05,
461 2.854e-07,
462 -1.286e-06,
463 -6.725e-01,
464 -3.730e-02,
465 8.964e-04,
466 1.399e-04,
467 -3.990e-06,
468 7.431e-06,
469 -2.796e-07,
470 -1.601e-07,
471 4.068e-02,
472 -1.352e-02,
473 7.282e-04,
474 9.594e-05,
475 2.070e-06,
476 -9.620e-08,
477 -2.742e-07,
478 -6.370e-08,
479 -6.300e-09,
480 8.625e-02,
481 -5.971e-03,
482 4.705e-04,
483 2.335e-05,
484 4.226e-06,
485 2.475e-07,
486 -8.850e-08,
487 -3.600e-08,
488 -2.900e-09,
489 0.000e+00
490 };
491
492
493 private static final double[] BH_AMPL = {
494 0.000e+00,
495 0.000e+00,
496 -1.136e-01,
497 0.000e+00,
498 -1.868e-01,
499 -1.399e-02,
500 0.000e+00,
501 -1.043e-01,
502 1.175e-02,
503 -2.240e-03,
504 0.000e+00,
505 -3.222e-02,
506 1.333e-02,
507 -2.647e-03,
508 -2.316e-05,
509 0.000e+00,
510 5.339e-02,
511 1.107e-02,
512 -3.116e-03,
513 -1.079e-04,
514 -1.299e-05,
515 0.000e+00,
516 4.861e-03,
517 8.891e-03,
518 -6.448e-04,
519 -1.279e-05,
520 6.358e-06,
521 -1.417e-07,
522 0.000e+00,
523 3.041e-02,
524 1.150e-03,
525 -8.743e-04,
526 -2.781e-05,
527 6.367e-07,
528 -1.140e-08,
529 -4.200e-08,
530 0.000e+00,
531 -2.982e-02,
532 -3.000e-03,
533 1.394e-05,
534 -3.290e-05,
535 -1.705e-07,
536 7.440e-08,
537 2.720e-08,
538 -6.600e-09,
539 0.000e+00,
540 1.236e-02,
541 -9.981e-04,
542 -3.792e-05,
543 -1.355e-05,
544 1.162e-06,
545 -1.789e-07,
546 1.470e-08,
547 -2.400e-09,
548 -4.000e-10
549 };
550
551
552 private static final double[] AW_MEAN = {
553 5.640e+01,
554 1.555e+00,
555 -1.011e+00,
556 -3.975e+00,
557 3.171e-02,
558 1.065e-01,
559 6.175e-01,
560 1.376e-01,
561 4.229e-02,
562 3.028e-03,
563 1.688e+00,
564 -1.692e-01,
565 5.478e-02,
566 2.473e-02,
567 6.059e-04,
568 2.278e+00,
569 6.614e-03,
570 -3.505e-04,
571 -6.697e-03,
572 8.402e-04,
573 7.033e-04,
574 -3.236e+00,
575 2.184e-01,
576 -4.611e-02,
577 -1.613e-02,
578 -1.604e-03,
579 5.420e-05,
580 7.922e-05,
581 -2.711e-01,
582 -4.406e-01,
583 -3.376e-02,
584 -2.801e-03,
585 -4.090e-04,
586 -2.056e-05,
587 6.894e-06,
588 2.317e-06,
589 1.941e+00,
590 -2.562e-01,
591 1.598e-02,
592 5.449e-03,
593 3.544e-04,
594 1.148e-05,
595 7.503e-06,
596 -5.667e-07,
597 -3.660e-08,
598 8.683e-01,
599 -5.931e-02,
600 -1.864e-03,
601 -1.277e-04,
602 2.029e-04,
603 1.269e-05,
604 1.629e-06,
605 9.660e-08,
606 -1.015e-07,
607 -5.000e-10
608 };
609
610
611 private static final double[] BW_MEAN = {
612 0.000e+00,
613 0.000e+00,
614 2.592e-01,
615 0.000e+00,
616 2.974e-02,
617 -5.471e-01,
618 0.000e+00,
619 -5.926e-01,
620 -1.030e-01,
621 -1.567e-02,
622 0.000e+00,
623 1.710e-01,
624 9.025e-02,
625 2.689e-02,
626 2.243e-03,
627 0.000e+00,
628 3.439e-01,
629 2.402e-02,
630 5.410e-03,
631 1.601e-03,
632 9.669e-05,
633 0.000e+00,
634 9.502e-02,
635 -3.063e-02,
636 -1.055e-03,
637 -1.067e-04,
638 -1.130e-04,
639 2.124e-05,
640 0.000e+00,
641 -3.129e-01,
642 8.463e-03,
643 2.253e-04,
644 7.413e-05,
645 -9.376e-05,
646 -1.606e-06,
647 2.060e-06,
648 0.000e+00,
649 2.739e-01,
650 1.167e-03,
651 -2.246e-05,
652 -1.287e-04,
653 -2.438e-05,
654 -7.561e-07,
655 1.158e-06,
656 4.950e-08,
657 0.000e+00,
658 -1.344e-01,
659 5.342e-03,
660 3.775e-04,
661 -6.756e-05,
662 -1.686e-06,
663 -1.184e-06,
664 2.768e-07,
665 2.730e-08,
666 5.700e-09
667 };
668
669
670 private static final double[] AW_AMPL = {
671 1.023e-01,
672 -2.695e+00,
673 3.417e-01,
674 -1.405e-01,
675 3.175e-01,
676 2.116e-01,
677 3.536e+00,
678 -1.505e-01,
679 -1.660e-02,
680 2.967e-02,
681 3.819e-01,
682 -1.695e-01,
683 -7.444e-02,
684 7.409e-03,
685 -6.262e-03,
686 -1.836e+00,
687 -1.759e-02,
688 -6.256e-02,
689 -2.371e-03,
690 7.947e-04,
691 1.501e-04,
692 -8.603e-01,
693 -1.360e-01,
694 -3.629e-02,
695 -3.706e-03,
696 -2.976e-04,
697 1.857e-05,
698 3.021e-05,
699 2.248e+00,
700 -1.178e-01,
701 1.255e-02,
702 1.134e-03,
703 -2.161e-04,
704 -5.817e-06,
705 8.836e-07,
706 -1.769e-07,
707 7.313e-01,
708 -1.188e-01,
709 1.145e-02,
710 1.011e-03,
711 1.083e-04,
712 2.570e-06,
713 -2.140e-06,
714 -5.710e-08,
715 2.000e-08,
716 -1.632e+00,
717 -6.948e-03,
718 -3.893e-03,
719 8.592e-04,
720 7.577e-05,
721 4.539e-06,
722 -3.852e-07,
723 -2.213e-07,
724 -1.370e-08,
725 5.800e-09
726 };
727
728
729 private static final double[] BW_AMPL = {
730 0.000e+00,
731 0.000e+00,
732 -8.865e-02,
733 0.000e+00,
734 -4.309e-01,
735 6.340e-02,
736 0.000e+00,
737 1.162e-01,
738 6.176e-02,
739 -4.234e-03,
740 0.000e+00,
741 2.530e-01,
742 4.017e-02,
743 -6.204e-03,
744 4.977e-03,
745 0.000e+00,
746 -1.737e-01,
747 -5.638e-03,
748 1.488e-04,
749 4.857e-04,
750 -1.809e-04,
751 0.000e+00,
752 -1.514e-01,
753 -1.685e-02,
754 5.333e-03,
755 -7.611e-05,
756 2.394e-05,
757 8.195e-06,
758 0.000e+00,
759 9.326e-02,
760 -1.275e-02,
761 -3.071e-04,
762 5.374e-05,
763 -3.391e-05,
764 -7.436e-06,
765 6.747e-07,
766 0.000e+00,
767 -8.637e-02,
768 -3.807e-03,
769 -6.833e-04,
770 -3.861e-05,
771 -2.268e-05,
772 1.454e-06,
773 3.860e-07,
774 -1.068e-07,
775 0.000e+00,
776 -2.658e-02,
777 -1.947e-03,
778 7.131e-04,
779 -3.506e-05,
780 1.885e-07,
781 5.792e-07,
782 3.990e-08,
783 2.000e-08,
784 -5.700e-09
785 };
786
787
788 ABCoefficients() {
789
790 }
791
792
793
794
795
796 public double getAHMean(final int index) {
797 return AH_MEAN[index];
798 }
799
800
801
802
803
804 public double getBHMean(final int index) {
805 return BH_MEAN[index];
806 }
807
808
809
810
811
812 public double getAWMean(final int index) {
813 return AW_MEAN[index];
814 }
815
816
817
818
819
820 public double getBWMean(final int index) {
821 return BW_MEAN[index];
822 }
823
824
825
826
827
828 public double getAHAmplitude(final int index) {
829 return AH_AMPL[index];
830 }
831
832
833
834
835
836 public double getBHAmplitude(final int index) {
837 return BH_AMPL[index];
838 }
839
840
841
842
843
844 public double getAWAmplitude(final int index) {
845 return AW_AMPL[index];
846 }
847
848
849
850
851
852 public double getBWAmplitude(final int index) {
853 return BW_AMPL[index];
854 }
855 }
856 }