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 org.hipparchus.Field;
20 import org.hipparchus.CalculusFieldElement;
21 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
22 import org.hipparchus.geometry.euclidean.threed.Vector3D;
23 import org.hipparchus.util.FastMath;
24 import org.hipparchus.util.FieldSinCos;
25 import org.hipparchus.util.MathArrays;
26 import org.hipparchus.util.MathUtils;
27 import org.hipparchus.util.SinCos;
28 import org.orekit.annotation.DefaultDataContext;
29 import org.orekit.bodies.BodyShape;
30 import org.orekit.bodies.FieldGeodeticPoint;
31 import org.orekit.bodies.GeodeticPoint;
32 import org.orekit.data.DataContext;
33 import org.orekit.errors.OrekitException;
34 import org.orekit.errors.OrekitMessages;
35 import org.orekit.frames.Frame;
36 import org.orekit.time.AbsoluteDate;
37 import org.orekit.time.DateTimeComponents;
38 import org.orekit.time.FieldAbsoluteDate;
39 import org.orekit.time.TimeScale;
40 import org.orekit.utils.Constants;
41 import org.orekit.utils.ExtendedPositionProvider;
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78 public class JB2008 extends AbstractSunInfluencedAtmosphere {
79
80
81 private static final double ALT_MIN = 90000.;
82
83
84 private static final double EARTH_RADIUS = 6356.766;
85
86
87 private static final double LOG10 = FastMath.log(10.);
88
89
90 private static final double AVOGAD = 6.02257e26;
91
92
93 private static final double RSTAR = 8.31432;
94
95
96 private static final double[] ALPHA = {
97 0, 0, 0, 0, -0.38
98 };
99
100
101 private static final double[] AMW = {
102 28.0134, 31.9988, 15.9994, 39.9480, 4.0026, 1.00797
103 };
104
105
106 private static final double[] FRAC = {
107 0.78110, 0.20955, 9.3400e-3, 1.2890e-5
108 };
109
110
111 private static final double R1 = 0.010;
112
113
114 private static final double R2 = 0.025;
115
116
117 private static final double R3 = 0.075;
118
119
120 private static final double[] WT = {
121 0.311111111111111, 1.422222222222222, 0.533333333333333, 1.422222222222222, 0.311111111111111
122 };
123
124
125 private static final double[] CHT = {
126 0.22, -0.20e-02, 0.115e-02, -0.211e-05
127 };
128
129
130 private static final double[] FZM = {
131 0.2689e+00, -0.1176e-01, 0.2782e-01, -0.2782e-01, 0.3470e-03
132 };
133
134
135 private static final double[] GTM = {
136 -0.3633e+00, 0.8506e-01, 0.2401e+00, -0.1897e+00, -0.2554e+00,
137 -0.1790e-01, 0.5650e-03, -0.6407e-03, -0.3418e-02, -0.1252e-02
138 };
139
140
141 private static final double[] CMB = {
142 28.15204, -8.5586e-2, +1.2840e-4, -1.0056e-5, -1.0210e-5, +1.5044e-6, +9.9826e-8
143 };
144
145
146 private static final double[] BDTC = {
147 -0.457512297e+01, -0.512114909e+01, -0.693003609e+02,
148 0.203716701e+03, 0.703316291e+03, -0.194349234e+04,
149 0.110651308e+04, -0.174378996e+03, 0.188594601e+04,
150 -0.709371517e+04, 0.922454523e+04, -0.384508073e+04,
151 -0.645841789e+01, 0.409703319e+02, -0.482006560e+03,
152 0.181870931e+04, -0.237389204e+04, 0.996703815e+03,
153 0.361416936e+02
154 };
155
156
157 private static final double[] CDTC = {
158 -0.155986211e+02, -0.512114909e+01, -0.693003609e+02,
159 0.203716701e+03, 0.703316291e+03, -0.194349234e+04,
160 0.110651308e+04, -0.220835117e+03, 0.143256989e+04,
161 -0.318481844e+04, 0.328981513e+04, -0.135332119e+04,
162 0.199956489e+02, -0.127093998e+02, 0.212825156e+02,
163 -0.275555432e+01, 0.110234982e+02, 0.148881951e+03,
164 -0.751640284e+03, 0.637876542e+03, 0.127093998e+02,
165 -0.212825156e+02, 0.275555432e+01
166 };
167
168
169 private JB2008InputParameters inputParams;
170
171
172 private final BodyShape earth;
173
174
175 private final TimeScale utc;
176
177
178
179
180
181
182
183
184
185
186 @DefaultDataContext
187 public JB2008(final JB2008InputParameters parameters,
188 final ExtendedPositionProvider sun, final BodyShape earth) {
189 this(parameters, sun, earth, DataContext.getDefault().getTimeScales().getUTC());
190 }
191
192
193
194
195
196
197
198
199
200
201 public JB2008(final JB2008InputParameters parameters,
202 final ExtendedPositionProvider sun,
203 final BodyShape earth,
204 final TimeScale utc) {
205 super(sun);
206 this.earth = earth;
207 this.inputParams = parameters;
208 this.utc = utc;
209 }
210
211
212 public Frame getFrame() {
213 return earth.getBodyFrame();
214 }
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242 public double getDensity(final double dateMJD, final double sunRA, final double sunDecli,
243 final double satLon, final double satLat, final double satAlt,
244 final double f10, final double f10B, final double s10,
245 final double s10B, final double xm10, final double xm10B,
246 final double y10, final double y10B, final double dstdtc) {
247
248 if (satAlt < ALT_MIN) {
249 throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, satAlt, ALT_MIN);
250 }
251
252 final double altKm = satAlt / 1000.0;
253
254
255 final double fn = FastMath.min(1.0, FastMath.pow(f10B / 240., 0.25));
256 final double fsb = f10B * fn + s10B * (1. - fn);
257 final double tsubc = 392.4 + 3.227 * fsb + 0.298 * (f10 - f10B) + 2.259 * (s10 - s10B) +
258 0.312 * (xm10 - xm10B) + 0.178 * (y10 - y10B);
259
260
261 final double eta = 0.5 * FastMath.abs(satLat - sunDecli);
262 final double theta = 0.5 * FastMath.abs(satLat + sunDecli);
263
264
265 final double h = satLon - sunRA;
266 final double tau = h - 0.64577182 + 0.10471976 * FastMath.sin(h + 0.75049158);
267 double solarTime = FastMath.toDegrees(h + FastMath.PI) / 15.0;
268 while (solarTime >= 24) {
269 solarTime -= 24.;
270 }
271 while (solarTime < 0) {
272 solarTime += 24.;
273 }
274
275
276 final double cosEta = FastMath.pow(FastMath.cos(eta), 2.5);
277 final double sinTeta = FastMath.pow(FastMath.sin(theta), 2.5);
278 final double cosTau = FastMath.abs(FastMath.cos(0.5 * tau));
279 final double df = sinTeta + (cosEta - sinTeta) * cosTau * cosTau * cosTau;
280 final double tsubl = tsubc * (1. + 0.31 * df);
281
282
283 final double dtclst = dTc(f10, solarTime, satLat, altKm);
284
285
286
287 final double[] temp = new double[2];
288 temp[0] = tsubl + dstdtc;
289 final double tinf = temp[0] + dtclst;
290
291
292 final double tsubx = 444.3807 + 0.02385 * tinf - 392.8292 * FastMath.exp(-0.0021357 * tinf);
293
294
295 final double gsubx = 0.054285714 * (tsubx - 183.);
296
297
298
299 final double[] tc = new double[4];
300 tc[0] = tsubx;
301 tc[1] = gsubx;
302
303 tc[2] = (tinf - tsubx) * 2. / FastMath.PI;
304 tc[3] = gsubx / tc[2];
305
306
307 final double z1 = 90.;
308 final double z2 = FastMath.min(altKm, 105.0);
309 double al = FastMath.log(z2 / z1);
310 int n = 1 + (int) FastMath.floor(al / R1);
311 double zr = FastMath.exp(al / n);
312 final double mb1 = mBar(z1);
313 final double tloc1 = localTemp(z1, tc);
314 double zend = z1;
315 double sub2 = 0.;
316 double ain = mb1 * gravity(z1) / tloc1;
317 double mb2 = 0;
318 double tloc2 = 0;
319 double z = 0;
320 double gravl = 0;
321
322 for (int i = 0; i < n; ++i) {
323 z = zend;
324 zend = zr * z;
325 final double dz = 0.25 * (zend - z);
326 double sum1 = WT[0] * ain;
327 for (int j = 1; j < 5; ++j) {
328 z += dz;
329 mb2 = mBar(z);
330 tloc2 = localTemp(z, tc);
331 gravl = gravity(z);
332 ain = mb2 * gravl / tloc2;
333 sum1 += WT[j] * ain;
334 }
335 sub2 += dz * sum1;
336 }
337
338 double rho = 3.46e-6 * mb2 * tloc1 / FastMath.exp(sub2 / RSTAR) / (mb1 * tloc2);
339
340
341 final double anm = AVOGAD * rho;
342 final double an = anm / mb2;
343
344
345 double fact2 = anm / 28.960;
346 final double[] aln = new double[6];
347 aln[0] = FastMath.log(FRAC[0] * fact2);
348 aln[3] = FastMath.log(FRAC[2] * fact2);
349 aln[4] = FastMath.log(FRAC[3] * fact2);
350
351 aln[1] = FastMath.log(fact2 * (1. + FRAC[1]) - an);
352 aln[2] = FastMath.log(2. * (an - fact2));
353
354 if (altKm <= 105.0) {
355 temp[1] = tloc2;
356
357 aln[5] = aln[4] - 25.0;
358 } else {
359
360 al = FastMath.log(FastMath.min(altKm, 500.0) / z);
361 n = 1 + (int) FastMath.floor(al / R2);
362 zr = FastMath.exp(al / n);
363 sub2 = 0.;
364 ain = gravl / tloc2;
365
366 double tloc3 = 0;
367 for (int i = 0; i < n; ++i) {
368 z = zend;
369 zend = zr * z;
370 final double dz = 0.25 * (zend - z);
371 double sum1 = WT[0] * ain;
372 for (int j = 1; j < 5; ++j) {
373 z += dz;
374 tloc3 = localTemp(z, tc);
375 gravl = gravity(z);
376 ain = gravl / tloc3;
377 sum1 += WT[j] * ain;
378 }
379 sub2 += dz * sum1;
380 }
381
382 al = FastMath.log(FastMath.max(altKm, 500.0) / z);
383 final double r = (altKm > 500.0) ? R3 : R2;
384 n = 1 + (int) FastMath.floor(al / r);
385 zr = FastMath.exp(al / n);
386 double sum3 = 0.;
387 double tloc4 = 0;
388 for (int i = 0; i < n; ++i) {
389 z = zend;
390 zend = zr * z;
391 final double dz = 0.25 * (zend - z);
392 double sum1 = WT[0] * ain;
393 for (int j = 1; j < 5; ++j) {
394 z += dz;
395 tloc4 = localTemp(z, tc);
396 gravl = gravity(z);
397 ain = gravl / tloc4;
398 sum1 += WT[j] * ain;
399 }
400 sum3 += dz * sum1;
401 }
402 final double altr;
403 final double hSign;
404 if (altKm <= 500.) {
405 temp[1] = tloc3;
406 altr = FastMath.log(tloc3 / tloc2);
407 fact2 = sub2 / RSTAR;
408 hSign = 1.0;
409 } else {
410 temp[1] = tloc4;
411 altr = FastMath.log(tloc4 / tloc2);
412 fact2 = (sub2 + sum3) / RSTAR;
413 hSign = -1.0;
414 }
415 for (int i = 0; i < 5; ++i) {
416 aln[i] = aln[i] - (1.0 + ALPHA[i]) * altr - fact2 * AMW[i];
417 }
418
419
420 final double al10t5 = FastMath.log10(tinf);
421 final double alnh5 = (5.5 * al10t5 - 39.40) * al10t5 + 73.13;
422 aln[5] = LOG10 * (alnh5 + 6.) + hSign * (FastMath.log(tloc4 / tloc3) + sum3 * AMW[5] / RSTAR);
423 }
424
425
426 final double capPhi = ((dateMJD - 36204.0) / 365.2422) % 1;
427 final int signum = (satLat >= 0.) ? 1 : -1;
428 final double sinLat = FastMath.sin(satLat);
429 final double hm90 = altKm - 90.;
430 final double dlrsl = 0.02 * hm90 * FastMath.exp(-0.045 * hm90) *
431 signum * FastMath.sin(MathUtils.TWO_PI * capPhi + 1.72) * sinLat * sinLat;
432
433
434 double dlrsa = 0;
435 if (z < 2000.0) {
436
437 dlrsa = semian08(dayOfYear(dateMJD), altKm, f10B, s10B, xm10B);
438 }
439
440
441
442
443
444 final double dlr = LOG10 * (dlrsl + dlrsa);
445 for (int i = 0; i < 6; ++i) {
446 aln[i] += dlr;
447 }
448
449
450
451 double sumnm = 0.0;
452 for (int i = 0; i < 6; ++i) {
453 sumnm += FastMath.exp(aln[i]) * AMW[i];
454 }
455 rho = sumnm / AVOGAD;
456
457
458 double fex = 1.0;
459 if (altKm >= 1000.0 && altKm < 1500.0) {
460 final double zeta = (altKm - 1000.) * 0.002;
461 final double f15c = CHT[0] + CHT[1] * f10B + (CHT[2] + CHT[3] * f10B) * 1500.0;
462 final double f15cZeta = (CHT[2] + CHT[3] * f10B) * 500.0;
463 final double fex2 = 3.0 * f15c - f15cZeta - 3.0;
464 final double fex3 = f15cZeta - 2.0 * f15c + 2.0;
465 fex += zeta * zeta * (fex2 + fex3 * zeta);
466 } else if (altKm >= 1500.0) {
467 fex = CHT[0] + CHT[1] * f10B + CHT[2] * altKm + CHT[3] * f10B * altKm;
468 }
469
470
471 rho *= fex;
472
473 return rho;
474 }
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503 public <T extends CalculusFieldElement<T>> T getDensity(final T dateMJD, final T sunRA, final T sunDecli,
504 final T satLon, final T satLat, final T satAlt,
505 final double f10, final double f10B, final double s10,
506 final double s10B, final double xm10, final double xm10B,
507 final double y10, final double y10B, final double dstdtc) {
508
509 if (satAlt.getReal() < ALT_MIN) {
510 throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD,
511 satAlt.getReal(), ALT_MIN);
512 }
513
514 final Field<T> field = satAlt.getField();
515 final T pi = field.getOne().getPi();
516 final T altKm = satAlt.divide(1000.0);
517
518
519 final double fn = FastMath.min(1.0, FastMath.pow(f10B / 240., 0.25));
520 final double fsb = f10B * fn + s10B * (1. - fn);
521 final double tsubc = 392.4 + 3.227 * fsb + 0.298 * (f10 - f10B) + 2.259 * (s10 - s10B) +
522 0.312 * (xm10 - xm10B) + 0.178 * (y10 - y10B);
523
524
525 final T eta = satLat.subtract(sunDecli).abs().multiply(0.5);
526 final T theta = satLat.add(sunDecli).abs().multiply(0.5);
527
528
529 final T h = satLon.subtract(sunRA);
530 final T tau = h.subtract(0.64577182).add(h.add(0.75049158).sin().multiply(0.10471976));
531 T solarTime = FastMath.toDegrees(h.add(pi)).divide(15.0);
532 while (solarTime.getReal() >= 24) {
533 solarTime = solarTime.subtract(24);
534 }
535 while (solarTime.getReal() < 0) {
536 solarTime = solarTime.add(24);
537 }
538
539
540 final T cos = eta.cos();
541 final T cosEta = cos.square().multiply(cos.sqrt());
542 final T sin = theta.sin();
543 final T sinTeta = sin.square().multiply(sin.sqrt());
544 final T cosTau = tau.multiply(0.5).cos().abs();
545 final T df = sinTeta.add(cosEta.subtract(sinTeta).multiply(cosTau).multiply(cosTau).multiply(cosTau));
546 final T tsubl = df.multiply(0.31).add(1).multiply(tsubc);
547
548
549 final T dtclst = dTc(f10, solarTime, satLat, altKm);
550
551
552
553 final T[] temp = MathArrays.buildArray(field, 2);
554 temp[0] = tsubl.add(dstdtc);
555 final T tinf = temp[0].add(dtclst);
556
557
558 final T tsubx = tinf.multiply(0.02385).add(444.3807).subtract(tinf.multiply(-0.0021357).exp().multiply(392.8292));
559
560
561 final T gsubx = tsubx.subtract(183.).multiply(0.054285714);
562
563
564
565 final T[] tc = MathArrays.buildArray(field, 4);
566 tc[0] = tsubx;
567 tc[1] = gsubx;
568
569 tc[2] = tinf.subtract(tsubx).multiply(pi.reciprocal().multiply(2.0));
570 tc[3] = gsubx.divide(tc[2]);
571
572
573 final T z1 = field.getZero().newInstance(90.);
574 final T z2 = min(105.0, altKm);
575 T al = z2.divide(z1).log();
576 int n = 1 + (int) FastMath.floor(al.getReal() / R1);
577 T zr = al.divide(n).exp();
578 final T mb1 = mBar(z1);
579 final T tloc1 = localTemp(z1, tc);
580 T zend = z1;
581 T sub2 = field.getZero();
582 T ain = mb1.multiply(gravity(z1)).divide(tloc1);
583 T mb2 = field.getZero();
584 T tloc2 = field.getZero();
585 T z = field.getZero();
586 T gravl = field.getZero();
587 for (int i = 0; i < n; ++i) {
588 z = zend;
589 zend = zr.multiply(z);
590 final T dz = zend.subtract(z).multiply(0.25);
591 T sum1 = ain.multiply(WT[0]);
592 for (int j = 1; j < 5; ++j) {
593 z = z.add(dz);
594 mb2 = mBar(z);
595 tloc2 = localTemp(z, tc);
596 gravl = gravity(z);
597 ain = mb2.multiply(gravl).divide(tloc2);
598 sum1 = sum1.add(ain.multiply(WT[j]));
599 }
600 sub2 = sub2.add(dz.multiply(sum1));
601 }
602
603 T rho = mb2.multiply(3.46e-6).multiply(tloc1).divide(sub2.divide(RSTAR).exp().multiply(mb1.multiply(tloc2)));
604
605
606 final T anm = rho.multiply(AVOGAD);
607 final T an = anm.divide(mb2);
608
609
610 T fact2 = anm.divide(28.960);
611 final T[] aln = MathArrays.buildArray(field, 6);
612 aln[0] = fact2.multiply(FRAC[0]).log();
613 aln[3] = fact2.multiply(FRAC[2]).log();
614 aln[4] = fact2.multiply(FRAC[3]).log();
615
616 aln[1] = fact2.multiply(1. + FRAC[1]).subtract(an).log();
617 aln[2] = an.subtract(fact2).multiply(2).log();
618
619 if (altKm.getReal() <= 105.0) {
620 temp[1] = tloc2;
621
622 aln[5] = aln[4].subtract(25.0);
623 } else {
624
625 al = min(500.0, altKm).divide(z).log();
626 n = 1 + (int) FastMath.floor(al.getReal() / R2);
627 zr = al.divide(n).exp();
628 sub2 = field.getZero();
629 ain = gravl.divide(tloc2);
630
631 T tloc3 = field.getZero();
632 for (int i = 0; i < n; ++i) {
633 z = zend;
634 zend = zr.multiply(z);
635 final T dz = zend.subtract(z).multiply(0.25);
636 T sum1 = ain.multiply(WT[0]);
637 for (int j = 1; j < 5; ++j) {
638 z = z.add(dz);
639 tloc3 = localTemp(z, tc);
640 gravl = gravity(z);
641 ain = gravl.divide(tloc3);
642 sum1 = sum1.add(ain.multiply(WT[j]));
643 }
644 sub2 = sub2.add(dz.multiply(sum1));
645 }
646
647 al = max(500.0, altKm).divide(z).log();
648 final double r = (altKm.getReal() > 500.0) ? R3 : R2;
649 n = 1 + (int) FastMath.floor(al.getReal() / r);
650 zr = al.divide(n).exp();
651 T sum3 = field.getZero();
652 T tloc4 = field.getZero();
653 for (int i = 0; i < n; ++i) {
654 z = zend;
655 zend = zr.multiply(z);
656 final T dz = zend.subtract(z).multiply(0.25);
657 T sum1 = ain.multiply(WT[0]);
658 for (int j = 1; j < 5; ++j) {
659 z = z.add(dz);
660 tloc4 = localTemp(z, tc);
661 gravl = gravity(z);
662 ain = gravl.divide(tloc4);
663 sum1 = sum1.add(ain.multiply(WT[j]));
664 }
665 sum3 = sum3.add(dz.multiply(sum1));
666 }
667 final T altr;
668 final double hSign;
669 if (altKm.getReal() <= 500.) {
670 temp[1] = tloc3;
671 altr = tloc3.divide(tloc2).log();
672 fact2 = sub2.divide(RSTAR);
673 hSign = 1.0;
674 } else {
675 temp[1] = tloc4;
676 altr = tloc4.divide(tloc2).log();
677 fact2 = sub2.add(sum3).divide(RSTAR);
678 hSign = -1.0;
679 }
680 for (int i = 0; i < 5; ++i) {
681 aln[i] = aln[i].subtract(altr.multiply(1.0 + ALPHA[i])).subtract(fact2.multiply(AMW[i]));
682 }
683
684
685 final T al10t5 = tinf.log10();
686 final T alnh5 = al10t5.multiply(5.5).subtract(39.40).multiply(al10t5).add(73.13);
687 aln[5] = alnh5.add(6.).multiply(LOG10).
688 add(tloc4.divide(tloc3).log().add(sum3.multiply(AMW[5] / RSTAR)).multiply(hSign));
689 }
690
691
692 T capPhi = dateMJD.subtract(36204.0).divide(365.2422);
693 capPhi = capPhi.subtract(FastMath.floor(capPhi.getReal()));
694 final int signum = (satLat.getReal() >= 0.) ? 1 : -1;
695 final T sinLat = satLat.sin();
696 final T hm90 = altKm.subtract(90.);
697 final T dlrsl = hm90.multiply(0.02).multiply(hm90.multiply(-0.045).exp()).
698 multiply(capPhi.multiply(MathUtils.TWO_PI).add(1.72).sin()).
699 multiply(signum).multiply(sinLat).multiply(sinLat);
700
701
702 T dlrsa = field.getZero();
703 if (z.getReal() < 2000.0) {
704
705 dlrsa = semian08(dayOfYear(dateMJD), altKm, f10B, s10B, xm10B);
706 }
707
708
709
710
711
712 final T dlr = dlrsl.add(dlrsa).multiply(LOG10);
713 for (int i = 0; i < 6; ++i) {
714 aln[i] = aln[i].add(dlr);
715 }
716
717
718
719 T sumnm = field.getZero();
720 for (int i = 0; i < 6; ++i) {
721 sumnm = sumnm.add(aln[i].exp().multiply(AMW[i]));
722 }
723 rho = sumnm.divide(AVOGAD);
724
725
726 T fex = field.getOne();
727 if (altKm.getReal() >= 1000.0 && altKm.getReal() < 1500.0) {
728 final T zeta = altKm.subtract(1000.).multiply(0.002);
729 final double f15c = CHT[0] + CHT[1] * f10B + (CHT[2] + CHT[3] * f10B) * 1500.0;
730 final double f15cZeta = (CHT[2] + CHT[3] * f10B) * 500.0;
731 final double fex2 = 3.0 * f15c - f15cZeta - 3.0;
732 final double fex3 = f15cZeta - 2.0 * f15c + 2.0;
733 fex = fex.add(zeta.multiply(zeta).multiply(zeta.multiply(fex3).add(fex2)));
734 } else if (altKm.getReal() >= 1500.0) {
735 fex = altKm.multiply(CHT[3] * f10B).add(altKm.multiply(CHT[2])).add(CHT[0] + CHT[1] * f10B);
736 }
737
738
739 rho = rho.multiply(fex);
740
741 return rho;
742 }
743
744
745
746
747
748
749
750
751 private static double dTc(final double f10, final double solarTime,
752 final double satLat, final double satAlt) {
753 double dTc = 0.;
754 final double st = solarTime / 24.0;
755 final double cs = FastMath.cos(satLat);
756 final double fs = (f10 - 100.0) / 100.0;
757
758
759 if (satAlt >= 120 && satAlt <= 200) {
760 final double dtc200 = poly2CDTC(fs, st, cs);
761 final double dtc200dz = poly1CDTC(fs, st, cs);
762 final double cc = 3.0 * dtc200 - dtc200dz;
763 final double dd = dtc200 - cc;
764 final double zp = (satAlt - 120.0) / 80.0;
765 dTc = zp * zp * (cc + dd * zp);
766 } else if (satAlt > 200.0 && satAlt <= 240.0) {
767 final double h = (satAlt - 200.0) / 50.0;
768 dTc = poly1CDTC(fs, st, cs) * h + poly2CDTC(fs, st, cs);
769 } else if (satAlt > 240.0 && satAlt <= 300.0) {
770 final double h = 0.8;
771 final double bb = poly1CDTC(fs, st, cs);
772 final double aa = bb * h + poly2CDTC(fs, st, cs);
773 final double p2BDT = poly2BDTC(st);
774 final double dtc300 = poly1BDTC(fs, st, cs, 3 * p2BDT);
775 final double dtc300dz = cs * p2BDT;
776 final double cc = 3.0 * dtc300 - dtc300dz - 3.0 * aa - 2.0 * bb;
777 final double dd = dtc300 - aa - bb - cc;
778 final double zp = (satAlt - 240.0) / 60.0;
779 dTc = aa + zp * (bb + zp * (cc + zp * dd));
780 } else if (satAlt > 300.0 && satAlt <= 600.0) {
781 final double h = satAlt / 100.0;
782 dTc = poly1BDTC(fs, st, cs, h * poly2BDTC(st));
783 } else if (satAlt > 600.0 && satAlt <= 800.0) {
784 final double poly2 = poly2BDTC(st);
785 final double aa = poly1BDTC(fs, st, cs, 6 * poly2);
786 final double bb = cs * poly2;
787 final double cc = -(3.0 * aa + 4.0 * bb) / 4.0;
788 final double dd = (aa + bb) / 4.0;
789 final double zp = (satAlt - 600.0) / 100.0;
790 dTc = aa + zp * (bb + zp * (cc + zp * dd));
791 }
792
793 return dTc;
794 }
795
796
797
798
799
800
801
802
803
804 private static <T extends CalculusFieldElement<T>> T dTc(final double f10, final T solarTime,
805 final T satLat, final T satAlt) {
806 T dTc = solarTime.getField().getZero();
807 final T st = solarTime.divide(24.0);
808 final T cs = satLat.cos();
809 final double fs = (f10 - 100.0) / 100.0;
810
811
812 if (satAlt.getReal() >= 120 && satAlt.getReal() <= 200) {
813 final T dtc200 = poly2CDTC(fs, st, cs);
814 final T dtc200dz = poly1CDTC(fs, st, cs);
815 final T cc = dtc200.multiply(3).subtract(dtc200dz);
816 final T dd = dtc200.subtract(cc);
817 final T zp = satAlt.subtract(120.0).divide(80.0);
818 dTc = zp.multiply(zp).multiply(cc.add(dd.multiply(zp)));
819 } else if (satAlt.getReal() > 200.0 && satAlt.getReal() <= 240.0) {
820 final T h = satAlt.subtract(200.0).divide(50.0);
821 dTc = poly1CDTC(fs, st, cs).multiply(h).add(poly2CDTC(fs, st, cs));
822 } else if (satAlt.getReal() > 240.0 && satAlt.getReal() <= 300.0) {
823 final T h = solarTime.getField().getZero().newInstance(0.8);
824 final T bb = poly1CDTC(fs, st, cs);
825 final T aa = bb.multiply(h).add(poly2CDTC(fs, st, cs));
826 final T p2BDT = poly2BDTC(st);
827 final T dtc300 = poly1BDTC(fs, st, cs, p2BDT.multiply(3));
828 final T dtc300dz = cs.multiply(p2BDT);
829 final T cc = dtc300.multiply(3).subtract(dtc300dz).subtract(aa.multiply(3)).subtract(bb.multiply(2));
830 final T dd = dtc300.subtract(aa).subtract(bb).subtract(cc);
831 final T zp = satAlt.subtract(240.0).divide(60.0);
832 dTc = aa.add(zp.multiply(bb.add(zp.multiply(cc.add(zp.multiply(dd))))));
833 } else if (satAlt.getReal() > 300.0 && satAlt.getReal() <= 600.0) {
834 final T h = satAlt.divide(100.0);
835 dTc = poly1BDTC(fs, st, cs, h.multiply(poly2BDTC(st)));
836 } else if (satAlt.getReal() > 600.0 && satAlt.getReal() <= 800.0) {
837 final T poly2 = poly2BDTC(st);
838 final T aa = poly1BDTC(fs, st, cs, poly2.multiply(6));
839 final T bb = cs.multiply(poly2);
840 final T cc = aa.multiply(3).add(bb.multiply(4)).divide(-4.0);
841 final T dd = aa.add(bb).divide(4.0);
842 final T zp = satAlt.subtract(600.0).divide(100.0);
843 dTc = aa.add(zp.multiply(bb.add(zp.multiply(cc.add(zp.multiply(dd))))));
844 }
845
846 return dTc;
847 }
848
849
850
851
852
853
854
855 private static double poly1CDTC(final double fs, final double st, final double cs) {
856 return CDTC[0] +
857 fs * (CDTC[1] + st * (CDTC[2] + st * (CDTC[3] + st * (CDTC[4] + st * (CDTC[5] + st * CDTC[6]))))) +
858 cs * st * (CDTC[7] + st * (CDTC[8] + st * (CDTC[9] + st * (CDTC[10] + st * CDTC[11])))) +
859 cs * (CDTC[12] + fs * (CDTC[13] + st * (CDTC[14] + st * CDTC[15])));
860 }
861
862
863
864
865
866
867
868
869 private static <T extends CalculusFieldElement<T>> T poly1CDTC(final double fs, final T st, final T cs) {
870 return st.multiply(CDTC[6]).
871 add(CDTC[5]).multiply(st).
872 add(CDTC[4]).multiply(st).
873 add(CDTC[3]).multiply(st).
874 add(CDTC[2]).multiply(st).
875 add(CDTC[1]).multiply(fs).
876 add(st.multiply(CDTC[11]).
877 add(CDTC[10]).multiply(st).
878 add(CDTC[ 9]).multiply(st).
879 add(CDTC[ 8]).multiply(st).
880 add(CDTC[7]).multiply(st).multiply(cs)).
881 add(st.multiply(CDTC[15]).
882 add(CDTC[14]).multiply(st).
883 add(CDTC[13]).multiply(fs).
884 add(CDTC[12]).multiply(cs)).
885 add(CDTC[0]);
886 }
887
888
889
890
891
892
893
894 private static double poly2CDTC(final double fs, final double st, final double cs) {
895 return CDTC[16] + st * cs * (CDTC[17] + st * (CDTC[18] + st * CDTC[19])) +
896 fs * cs * (CDTC[20] + st * (CDTC[21] + st * CDTC[22]));
897 }
898
899
900
901
902
903
904
905
906 private static <T extends CalculusFieldElement<T>> T poly2CDTC(final double fs, final T st, final T cs) {
907 return st.multiply(CDTC[19]).
908 add(CDTC[18]).multiply(st).
909 add(CDTC[17]).multiply(cs).multiply(st).
910 add( st.multiply(CDTC[22]).
911 add(CDTC[21]).multiply(st).
912 add(CDTC[20]).multiply(cs).multiply(fs)).
913 add(CDTC[16]);
914 }
915
916
917
918
919
920
921
922
923 private static double poly1BDTC(final double fs, final double st, final double cs, final double hp) {
924 return BDTC[0] +
925 fs * (BDTC[1] + st * (BDTC[2] + st * (BDTC[3] + st * (BDTC[4] + st * (BDTC[5] + st * BDTC[6]))))) +
926 cs * (st * (BDTC[7] + st * (BDTC[8] + st * (BDTC[9] + st * (BDTC[10] + st * BDTC[11])))) + hp + BDTC[18]);
927 }
928
929
930
931
932
933
934
935
936
937 private static <T extends CalculusFieldElement<T>> T poly1BDTC(final double fs, final T st, final T cs, final T hp) {
938 return st.multiply(BDTC[6]).
939 add(BDTC[5]).multiply(st).
940 add(BDTC[4]).multiply(st).
941 add(BDTC[3]).multiply(st).
942 add(BDTC[2]).multiply(st).
943 add(BDTC[1]).multiply(fs).
944 add( st.multiply(BDTC[11]).
945 add(BDTC[10]).multiply(st).
946 add(BDTC[ 9]).multiply(st).
947 add(BDTC[ 8]).multiply(st).
948 add(BDTC[ 7]).multiply(st).
949 add(hp).add(BDTC[18]).multiply(cs)).
950 add(BDTC[0]);
951 }
952
953
954
955
956
957 private static double poly2BDTC(final double st) {
958 return BDTC[12] + st * (BDTC[13] + st * (BDTC[14] + st * (BDTC[15] + st * (BDTC[16] + st * BDTC[17]))));
959 }
960
961
962
963
964
965
966 private static <T extends CalculusFieldElement<T>> T poly2BDTC(final T st) {
967 return st.multiply(BDTC[17]).
968 add(BDTC[16]).multiply(st).
969 add(BDTC[15]).multiply(st).
970 add(BDTC[14]).multiply(st).
971 add(BDTC[13]).multiply(st).
972 add(BDTC[12]);
973 }
974
975
976
977
978
979 private static double mBar(final double z) {
980 final double dz = z - 100.;
981 double amb = CMB[6];
982 for (int i = 5; i >= 0; --i) {
983 amb = dz * amb + CMB[i];
984 }
985 return amb;
986 }
987
988
989
990
991
992
993 private static <T extends CalculusFieldElement<T>> T mBar(final T z) {
994 final T dz = z.subtract(100.);
995 T amb = z.getField().getZero().newInstance(CMB[6]);
996 for (int i = 5; i >= 0; --i) {
997 amb = dz.multiply(amb).add(CMB[i]);
998 }
999 return amb;
1000 }
1001
1002
1003
1004
1005
1006
1007 private static double localTemp(final double z, final double[] tc) {
1008 final double dz = z - 125.;
1009 if (dz <= 0.) {
1010 return ((-9.8204695e-6 * dz - 7.3039742e-4) * dz * dz + 1.0) * dz * tc[1] + tc[0];
1011 } else {
1012 return tc[0] + tc[2] * FastMath.atan(tc[3] * dz * (1 + 4.5e-6 * FastMath.pow(dz, 2.5)));
1013 }
1014 }
1015
1016
1017
1018
1019
1020
1021
1022 private static <T extends CalculusFieldElement<T>> T localTemp(final T z, final T[] tc) {
1023 final T dz = z.subtract(125.);
1024 if (dz.getReal() <= 0.) {
1025 return dz.multiply(-9.8204695e-6).subtract(7.3039742e-4).multiply(dz).multiply(dz).add(1.0).multiply(dz).multiply(tc[1]).add(tc[0]);
1026 } else {
1027 return dz.multiply(dz).multiply(dz.sqrt()).multiply(4.5e-6).add(1).multiply(dz).multiply(tc[3]).atan().multiply(tc[2]).add(tc[0]);
1028 }
1029 }
1030
1031
1032
1033
1034
1035 private static double gravity(final double z) {
1036 final double tmp = 1.0 + z / EARTH_RADIUS;
1037 return Constants.G0_STANDARD_GRAVITY / (tmp * tmp);
1038 }
1039
1040
1041
1042
1043
1044
1045 private static <T extends CalculusFieldElement<T>> T gravity(final T z) {
1046 final T tmp = z.divide(EARTH_RADIUS).add(1);
1047 return tmp.multiply(tmp).reciprocal().multiply(Constants.G0_STANDARD_GRAVITY);
1048 }
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058 private static double semian08(final double doy, final double alt,
1059 final double f10B, final double s10B, final double xm10B) {
1060
1061 final double htz = alt / 1000.0;
1062
1063
1064 double fsmb = f10B - 0.70 * s10B - 0.04 * xm10B;
1065
1066
1067 final double fzz = FZM[0] + fsmb * (FZM[1] + htz * (FZM[2] + FZM[3] * htz + FZM[4] * fsmb));
1068
1069
1070 fsmb = f10B - 0.75 * s10B - 0.37 * xm10B;
1071
1072
1073 final double tau = MathUtils.TWO_PI * (doy - 1.0) / 365;
1074 final SinCos sc1P = FastMath.sinCos(tau);
1075 final SinCos sc2P = SinCos.sum(sc1P, sc1P);
1076 final double gtz = GTM[0] + GTM[1] * sc1P.sin() + GTM[2] * sc1P.cos() + GTM[3] * sc2P.sin() + GTM[4] * sc2P.cos() +
1077 fsmb * (GTM[5] + GTM[6] * sc1P.sin() + GTM[7] * sc1P.cos() + GTM[8] * sc2P.sin() + GTM[9] * sc2P.cos());
1078
1079 return FastMath.max(1.0e-6, fzz) * gtz;
1080
1081 }
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092 private static <T extends CalculusFieldElement<T>> T semian08(final T doy, final T alt,
1093 final double f10B, final double s10B, final double xm10B) {
1094
1095 final T htz = alt.divide(1000.0);
1096
1097
1098 double fsmb = f10B - 0.70 * s10B - 0.04 * xm10B;
1099
1100
1101 final T fzz = htz.multiply(FZM[3]).add(FZM[2] + FZM[4] * fsmb).multiply(htz).add(FZM[1]).multiply(fsmb).add(FZM[0]);
1102
1103
1104 fsmb = f10B - 0.75 * s10B - 0.37 * xm10B;
1105
1106
1107 final T tau = doy.subtract(1).divide(365).multiply(MathUtils.TWO_PI);
1108 final FieldSinCos<T> sc1P = FastMath.sinCos(tau);
1109 final FieldSinCos<T> sc2P = FieldSinCos.sum(sc1P, sc1P);
1110 final T gtz = sc2P.cos().multiply(GTM[9]).
1111 add(sc2P.sin().multiply(GTM[8])).
1112 add(sc1P.cos().multiply(GTM[7])).
1113 add(sc1P.sin().multiply(GTM[6])).
1114 add(GTM[5]).multiply(fsmb).
1115 add( sc2P.cos().multiply(GTM[4]).
1116 add(sc2P.sin().multiply(GTM[3])).
1117 add(sc1P.cos().multiply(GTM[2])).
1118 add(sc1P.sin().multiply(GTM[1])).
1119 add(GTM[0]));
1120
1121 return fzz.getReal() > 1.0e-6 ? gtz.multiply(fzz) : gtz.multiply(1.0e-6);
1122
1123 }
1124
1125
1126
1127
1128
1129 private static double dayOfYear(final double dateMJD) {
1130 final double d1950 = dateMJD - 33281;
1131
1132 int iyday = (int) d1950;
1133 final double frac = d1950 - iyday;
1134 iyday = iyday + 364;
1135
1136 int itemp = iyday / 1461;
1137
1138 iyday = iyday - itemp * 1461;
1139 itemp = iyday / 365;
1140 if (itemp >= 3) {
1141 itemp = 3;
1142 }
1143 iyday = iyday - 365 * itemp + 1;
1144 return iyday + frac;
1145 }
1146
1147
1148
1149
1150
1151
1152 private static <T extends CalculusFieldElement<T>> T dayOfYear(final T dateMJD) {
1153 final T d1950 = dateMJD.subtract(33281);
1154
1155 int iyday = (int) d1950.getReal();
1156 final T frac = d1950.subtract(iyday);
1157 iyday = iyday + 364;
1158
1159 int itemp = iyday / 1461;
1160
1161 iyday = iyday - itemp * 1461;
1162 itemp = iyday / 365;
1163 if (itemp >= 3) {
1164 itemp = 3;
1165 }
1166 iyday = iyday - 365 * itemp + 1;
1167 return frac.add(iyday);
1168 }
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178 private <T extends CalculusFieldElement<T>> T min(final double d, final T f) {
1179 return (f.getReal() > d) ? f.getField().getZero().newInstance(d) : f;
1180 }
1181
1182
1183
1184
1185
1186
1187
1188 private <T extends CalculusFieldElement<T>> T max(final double d, final T f) {
1189 return (f.getReal() <= d) ? f.getField().getZero().newInstance(d) : f;
1190 }
1191
1192
1193
1194
1195
1196
1197
1198 public double getDensity(final AbsoluteDate date, final Vector3D position,
1199 final Frame frame) {
1200
1201 if (date.compareTo(inputParams.getMaxDate()) > 0 ||
1202 date.compareTo(inputParams.getMinDate()) < 0) {
1203 throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
1204 date, inputParams.getMinDate(), inputParams.getMaxDate());
1205 }
1206
1207
1208 final DateTimeComponents dt = date.getComponents(utc);
1209 final double dateMJD = dt.getDate().getMJD() +
1210 dt.getTime().getSecondsInLocalDay() / Constants.JULIAN_DAY;
1211
1212
1213 final GeodeticPoint inBody = earth.transform(position, frame, date);
1214
1215
1216 final Frame ecef = earth.getBodyFrame();
1217 final Vector3D sunPos = getSunPosition(date, ecef);
1218 final GeodeticPoint sunInBody = earth.transform(sunPos, ecef, date);
1219
1220 return getDensity(dateMJD,
1221 sunInBody.getLongitude(), sunInBody.getLatitude(),
1222 inBody.getLongitude(), inBody.getLatitude(), inBody.getAltitude(),
1223 inputParams.getF10(date), inputParams.getF10B(date),
1224 inputParams.getS10(date), inputParams.getS10B(date),
1225 inputParams.getXM10(date), inputParams.getXM10B(date),
1226 inputParams.getY10(date), inputParams.getY10B(date),
1227 inputParams.getDSTDTC(date));
1228
1229 }
1230
1231
1232 @Override
1233 public <T extends CalculusFieldElement<T>> T getDensity(final FieldAbsoluteDate<T> date,
1234 final FieldVector3D<T> position,
1235 final Frame frame) {
1236
1237
1238 final AbsoluteDate dateD = date.toAbsoluteDate();
1239 if (dateD.compareTo(inputParams.getMaxDate()) > 0 ||
1240 dateD.compareTo(inputParams.getMinDate()) < 0) {
1241 throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
1242 dateD, inputParams.getMinDate(), inputParams.getMaxDate());
1243 }
1244
1245
1246 final DateTimeComponents components = date.getComponents(utc);
1247 final T dateMJD = date
1248 .durationFrom(new FieldAbsoluteDate<>(date.getField(), components, utc))
1249 .add(components.getTime().getSecondsInLocalDay())
1250 .divide(Constants.JULIAN_DAY)
1251 .add(components.getDate().getMJD());
1252
1253
1254 final FieldGeodeticPoint<T> inBody = earth.transform(position, frame, date);
1255
1256
1257 final Frame ecef = earth.getBodyFrame();
1258 final FieldVector3D<T> sunPos = getSunPosition(date, ecef);
1259 final FieldGeodeticPoint<T> sunInBody = earth.transform(sunPos, ecef, date);
1260
1261 return getDensity(dateMJD,
1262 sunInBody.getLongitude(), sunInBody.getLatitude(),
1263 inBody.getLongitude(), inBody.getLatitude(), inBody.getAltitude(),
1264 inputParams.getF10(dateD), inputParams.getF10B(dateD),
1265 inputParams.getS10(dateD), inputParams.getS10B(dateD),
1266 inputParams.getXM10(dateD), inputParams.getXM10B(dateD),
1267 inputParams.getY10(dateD), inputParams.getY10B(dateD),
1268 inputParams.getDSTDTC(dateD));
1269
1270 }
1271
1272 }