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