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.CalculusFieldElement;
20 import org.hipparchus.exception.DummyLocalizable;
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.SinCos;
27 import org.orekit.annotation.DefaultDataContext;
28 import org.orekit.bodies.BodyShape;
29 import org.orekit.bodies.FieldGeodeticPoint;
30 import org.orekit.bodies.GeodeticPoint;
31 import org.orekit.data.DataContext;
32 import org.orekit.errors.OrekitException;
33 import org.orekit.errors.OrekitMessages;
34 import org.orekit.frames.Frame;
35 import org.orekit.time.AbsoluteDate;
36 import org.orekit.time.FieldAbsoluteDate;
37 import org.orekit.time.TimeScale;
38 import org.orekit.utils.PVCoordinatesProvider;
39
40 import java.io.BufferedReader;
41 import java.io.IOException;
42 import java.io.InputStream;
43 import java.io.InputStreamReader;
44 import java.nio.charset.StandardCharsets;
45 import java.util.Arrays;
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
79
80
81
82
83 public class DTM2000 implements Atmosphere {
84
85
86 private static final long serialVersionUID = 20170705L;
87
88
89
90
91 private static final int NLATM = 96;
92
93
94 private static final double[] ALEFA = {
95 0, -0.40, -0.38, 0, 0, 0, 0
96 };
97
98
99 private static final double[] MA = {
100 0, 1, 4, 16, 28, 32, 14
101 };
102
103
104 private static final double[] VMA = {
105 0, 1.6606e-24, 6.6423e-24, 26.569e-24, 46.4958e-24, 53.1381e-24, 23.2479e-24
106 };
107
108
109 private static final double RE = 6356.77;
110
111
112 private static final double ZLB0 = 120.0;
113
114
115 private static final double CPMG = .19081;
116
117
118 private static final double SPMG = .98163;
119
120
121 private static final double XLMG = -1.2392;
122
123
124 private static final double GSURF = 980.665;
125
126
127 private static final double RGAS = 831.4;
128
129
130 private static final double ROT = 0.017214206;
131
132
133 private static final double ROT2 = 0.034428412;
134
135
136 private static final String DTM2000 = "/assets/org/orekit/dtm_2000.txt";
137
138
139
140
141 private static double[] tt = null;
142 private static double[] h = null;
143 private static double[] he = null;
144 private static double[] o = null;
145 private static double[] az2 = null;
146 private static double[] o2 = null;
147 private static double[] az = null;
148 private static double[] t0 = null;
149 private static double[] tp = null;
150
151
152 private PVCoordinatesProvider sun;
153
154
155 private DTM2000InputParameters inputParams;
156
157
158 private BodyShape earth;
159
160
161 private final TimeScale utc;
162
163
164
165
166
167
168
169
170
171
172 @DefaultDataContext
173 public DTM2000(final DTM2000InputParameters parameters,
174 final PVCoordinatesProvider sun, final BodyShape earth) {
175 this(parameters, sun, earth, DataContext.getDefault().getTimeScales().getUTC());
176 }
177
178
179
180
181
182
183
184
185 public DTM2000(final DTM2000InputParameters parameters,
186 final PVCoordinatesProvider sun,
187 final BodyShape earth,
188 final TimeScale utc) {
189
190 synchronized (DTM2000.class) {
191
192 if (tt == null) {
193 readcoefficients();
194 }
195 }
196
197 this.earth = earth;
198 this.sun = sun;
199 this.inputParams = parameters;
200
201 this.utc = utc;
202 }
203
204
205 public Frame getFrame() {
206 return earth.getBodyFrame();
207 }
208
209
210
211
212
213
214
215
216
217
218
219
220
221 public double getDensity(final int day,
222 final double alti, final double lon, final double lat,
223 final double hl, final double f, final double fbar,
224 final double akp3, final double akp24) {
225 final double threshold = 120000;
226 if (alti < threshold) {
227 throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD,
228 alti, threshold);
229 }
230 final Computation result = new Computation(day, alti / 1000, lon, lat, hl,
231 new double[] {
232 0, f, 0
233 }, new double[] {
234 0, fbar, 0
235 }, new double[] {
236 0, akp3, 0, akp24, 0
237 });
238 return result.ro * 1000;
239 }
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255 public <T extends CalculusFieldElement<T>> T getDensity(final int day,
256 final T alti, final T lon, final T lat,
257 final T hl, final double f, final double fbar,
258 final double akp3, final double akp24) {
259 final double threshold = 120000;
260 if (alti.getReal() < threshold) {
261 throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD,
262 alti, threshold);
263 }
264 final FieldComputation<T> result = new FieldComputation<>(day, alti.divide(1000), lon, lat, hl,
265 new double[] {
266 0, f, 0
267 }, new double[] {
268 0, fbar, 0
269 }, new double[] {
270 0, akp3, 0, akp24, 0
271 });
272 return result.ro.multiply(1000);
273 }
274
275
276
277 private static void readcoefficients() {
278
279 final int size = NLATM + 1;
280 tt = new double[size];
281 h = new double[size];
282 he = new double[size];
283 o = new double[size];
284 az2 = new double[size];
285 o2 = new double[size];
286 az = new double[size];
287 t0 = new double[size];
288 tp = new double[size];
289
290 try (InputStream in = checkNull(DTM2000.class.getResourceAsStream(DTM2000));
291 BufferedReader r = new BufferedReader(new InputStreamReader(in, StandardCharsets.UTF_8))) {
292
293 r.readLine();
294 r.readLine();
295 for (String line = r.readLine(); line != null; line = r.readLine()) {
296 final int num = Integer.parseInt(line.substring(0, 4).replace(' ', '0'));
297 line = line.substring(4);
298 tt[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
299 line = line.substring(13 + 9);
300 h[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
301 line = line.substring(13 + 9);
302 he[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
303 line = line.substring(13 + 9);
304 o[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
305 line = line.substring(13 + 9);
306 az2[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
307 line = line.substring(13 + 9);
308 o2[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
309 line = line.substring(13 + 9);
310 az[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
311 line = line.substring(13 + 9);
312 t0[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
313 line = line.substring(13 + 9);
314 tp[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
315 }
316 } catch (IOException ioe) {
317 throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
318 }
319 }
320
321
322
323
324
325
326
327 public double getDensity(final AbsoluteDate date, final Vector3D position,
328 final Frame frame) {
329
330
331 if (date.compareTo(inputParams.getMaxDate()) > 0 ||
332 date.compareTo(inputParams.getMinDate()) < 0) {
333 throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
334 date, inputParams.getMinDate(), inputParams.getMaxDate());
335 }
336
337
338 final int day = date.getComponents(utc).getDate().getDayOfYear();
339
340 final Frame ecef = earth.getBodyFrame();
341 final Vector3D pEcef = frame.getStaticTransformTo(ecef, date)
342 .transformPosition(position);
343
344 final GeodeticPoint inBody = earth.transform(pEcef, ecef, date);
345 final double alti = inBody.getAltitude();
346 final double lon = inBody.getLongitude();
347 final double lat = inBody.getLatitude();
348
349
350 final Vector3D sunPos = sun.getPosition(date, ecef);
351 final double hl = FastMath.PI + FastMath.atan2(
352 sunPos.getX() * pEcef.getY() - sunPos.getY() * pEcef.getX(),
353 sunPos.getX() * pEcef.getX() + sunPos.getY() * pEcef.getY());
354
355
356 return getDensity(day, alti, lon, lat, hl, inputParams.getInstantFlux(date),
357 inputParams.getMeanFlux(date), inputParams.getThreeHourlyKP(date),
358 inputParams.get24HoursKp(date));
359
360 }
361
362
363 @Override
364 public <T extends CalculusFieldElement<T>> T
365 getDensity(final FieldAbsoluteDate<T> date, final FieldVector3D<T> position,
366 final Frame frame) {
367
368 final AbsoluteDate dateD = date.toAbsoluteDate();
369 if (dateD.compareTo(inputParams.getMaxDate()) > 0 ||
370 dateD.compareTo(inputParams.getMinDate()) < 0) {
371 throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
372 dateD, inputParams.getMinDate(), inputParams.getMaxDate());
373 }
374
375
376 final int day = date.getComponents(utc).getDate().getDayOfYear();
377
378 final Frame ecef = earth.getBodyFrame();
379 final FieldVector3D<T> pEcef = frame.getStaticTransformTo(ecef, date).transformPosition(position);
380
381 final FieldGeodeticPoint<T> inBody = earth.transform(pEcef, ecef, date);
382 final T alti = inBody.getAltitude();
383 final T lon = inBody.getLongitude();
384 final T lat = inBody.getLatitude();
385
386
387 final Vector3D sunPos = sun.getPosition(dateD, ecef);
388 final T y = pEcef.getY().multiply(sunPos.getX()).subtract(pEcef.getX().multiply(sunPos.getY()));
389 final T x = pEcef.getX().multiply(sunPos.getX()).add(pEcef.getY().multiply(sunPos.getY()));
390 final T hl = y.atan2(x).add(y.getPi());
391
392
393 return getDensity(day, alti, lon, lat, hl, inputParams.getInstantFlux(dateD),
394 inputParams.getMeanFlux(dateD), inputParams.getThreeHourlyKP(dateD),
395 inputParams.get24HoursKp(dateD));
396 }
397
398
399
400
401
402
403
404
405 private static InputStream checkNull(final InputStream stream) {
406 if (stream == null) {
407 throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_RESOURCE, DTM2000);
408 }
409 return stream;
410 }
411
412
413 private static class Computation {
414
415
416 private final int day;
417
418
419 private final double[] f;
420
421
422 private final double[] fbar;
423
424
425
426
427
428
429
430
431
432 private final double[] akp;
433
434
435 private final double clfl;
436
437
438 private final double slfl;
439
440
441 private final double ro;
442
443
444
445
446 private final double p10;
447 private final double p20;
448 private final double p30;
449 private final double p40;
450 private final double p50;
451 private final double p60;
452 private final double p11;
453 private final double p21;
454 private final double p31;
455 private final double p41;
456 private final double p51;
457 private final double p22;
458 private final double p32;
459 private final double p42;
460 private final double p52;
461 private final double p62;
462 private final double p33;
463 private final double p10mg;
464 private final double p20mg;
465 private final double p40mg;
466
467
468 private final double hl0;
469 private final double ch;
470 private final double sh;
471 private final double c2h;
472 private final double s2h;
473 private final double c3h;
474 private final double s3h;
475
476
477
478
479
480
481
482
483
484
485
486 Computation(final int day,
487 final double altiKM, final double lon, final double lat,
488 final double hl, final double[] f, final double[] fbar,
489 final double[] akp) {
490
491 this.day = day;
492 this.f = f;
493 this.fbar = fbar;
494 this.akp = akp;
495
496
497 final SinCos scLat = FastMath.sinCos(lat);
498 final SinCos scLon = FastMath.sinCos(lon);
499
500
501 final double c = scLat.sin();
502 final double c2 = c * c;
503 final double c4 = c2 * c2;
504 final double s = scLat.cos();
505 final double s2 = s * s;
506 p10 = c;
507 p20 = 1.5 * c2 - 0.5;
508 p30 = c * (2.5 * c2 - 1.5);
509 p40 = 4.375 * c4 - 3.75 * c2 + 0.375;
510 p50 = c * (7.875 * c4 - 8.75 * c2 + 1.875);
511 p60 = (5.5 * c * p50 - 2.5 * p40) / 3.0;
512 p11 = s;
513 p21 = 3.0 * c * s;
514 p31 = s * (7.5 * c2 - 1.5);
515 p41 = c * s * (17.5 * c2 - 7.5);
516 p51 = s * (39.375 * c4 - 26.25 * c2 + 1.875);
517 p22 = 3.0 * s2;
518 p32 = 15.0 * c * s2;
519 p42 = s2 * (52.5 * c2 - 7.5);
520 p52 = 3.0 * c * p42 - 2.0 * p32;
521 p62 = 2.75 * c * p52 - 1.75 * p42;
522 p33 = 15.0 * s * s2;
523
524
525 final double clmlmg = FastMath.cos(lon - XLMG);
526 final double cmg = s * CPMG * clmlmg + c * SPMG;
527 final double cmg2 = cmg * cmg;
528 final double cmg4 = cmg2 * cmg2;
529 p10mg = cmg;
530 p20mg = 1.5 * cmg2 - 0.5;
531 p40mg = 4.375 * cmg4 - 3.75 * cmg2 + 0.375;
532
533 clfl = scLon.cos();
534 slfl = scLon.sin();
535
536
537 hl0 = hl;
538 final SinCos scHlo = FastMath.sinCos(hl0);
539 ch = scHlo.cos();
540 sh = scHlo.sin();
541 c2h = ch * ch - sh * sh;
542 s2h = 2.0 * ch * sh;
543 c3h = c2h * ch - s2h * sh;
544 s3h = s2h * ch + c2h * sh;
545
546 final double zlb = ZLB0;
547
548 final double[] dtt = new double[tt.length];
549 final double[] dh = new double[tt.length];
550 final double[] dhe = new double[tt.length];
551 final double[] dox = new double[tt.length];
552 final double[] daz2 = new double[tt.length];
553 final double[] do2 = new double[tt.length];
554 final double[] daz = new double[tt.length];
555 final double[] dt0 = new double[tt.length];
556 final double[] dtp = new double[tt.length];
557
558 Arrays.fill(dtt, Double.NaN);
559 Arrays.fill(dh, Double.NaN);
560 Arrays.fill(dhe, Double.NaN);
561 Arrays.fill(dox, Double.NaN);
562 Arrays.fill(daz2, Double.NaN);
563 Arrays.fill(do2, Double.NaN);
564 Arrays.fill(daz, Double.NaN);
565 Arrays.fill(dt0, Double.NaN);
566 Arrays.fill(dtp, Double.NaN);
567
568
569 int kleq = 1;
570 final double gdelt = gFunction(tt, dtt, 1, kleq);
571 dtt[1] = 1.0 + gdelt;
572 final double tinf = tt[1] * dtt[1];
573
574 kleq = 0;
575
576 if (day < 59 || day > 284) {
577 kleq = -1;
578 }
579 if (day > 99 && day < 244) {
580 kleq = 1;
581 }
582
583 final double gdelt0 = gFunction(t0, dt0, 0, kleq);
584 dt0[1] = (t0[1] + gdelt0) / t0[1];
585 final double t120 = t0[1] + gdelt0;
586 final double gdeltp = gFunction(tp, dtp, 0, kleq);
587 dtp[1] = (tp[1] + gdeltp) / tp[1];
588 final double tp120 = tp[1] + gdeltp;
589
590
591 final double sigma = tp120 / (tinf - t120);
592 final double dzeta = (RE + zlb) / (RE + altiKM);
593 final double zeta = (altiKM - zlb) * dzeta;
594 final double sigzeta = sigma * zeta;
595 final double expsz = FastMath.exp(-sigzeta);
596 final double tz = tinf - (tinf - t120) * expsz;
597
598 final double[] dbase = new double[7];
599
600 kleq = 1;
601
602 final double gdelh = gFunction(h, dh, 0, kleq);
603 dh[1] = FastMath.exp(gdelh);
604 dbase[1] = h[1] * dh[1];
605
606 final double gdelhe = gFunction(he, dhe, 0, kleq);
607 dhe[1] = FastMath.exp(gdelhe);
608 dbase[2] = he[1] * dhe[1];
609
610 final double gdelo = gFunction(o, dox, 1, kleq);
611 dox[1] = FastMath.exp(gdelo);
612 dbase[3] = o[1] * dox[1];
613
614 final double gdelaz2 = gFunction(az2, daz2, 1, kleq);
615 daz2[1] = FastMath.exp(gdelaz2);
616 dbase[4] = az2[1] * daz2[1];
617
618 final double gdelo2 = gFunction(o2, do2, 1, kleq);
619 do2[1] = FastMath.exp(gdelo2);
620 dbase[5] = o2[1] * do2[1];
621
622 final double gdelaz = gFunction(az, daz, 1, kleq);
623 daz[1] = FastMath.exp(gdelaz);
624 dbase[6] = az[1] * daz[1];
625
626 final double zlbre = 1.0 + zlb / RE;
627 final double glb = (GSURF / (zlbre * zlbre)) / (sigma * RGAS * tinf);
628 final double t120tz = t120 / tz;
629
630
631
632
633
634
635
636
637 double tmpro = 0;
638 for (int i = 1; i <= 6; i++) {
639 final double gamma = MA[i] * glb;
640 final double upapg = 1.0 + ALEFA[i] + gamma;
641 final double fzI = FastMath.exp(FastMath.log(t120tz) * upapg - sigzeta * gamma);
642
643 final double ccI = dbase[i] * fzI;
644
645 tmpro += ccI * VMA[i];
646 }
647 this.ro = tmpro;
648
649 }
650
651
652
653
654
655
656
657
658 private double gFunction(final double[] a, final double[] da,
659 final int ff0, final int kle_eq) {
660
661 final double[] fmfb = new double[3];
662 final double[] fbm150 = new double[3];
663
664
665 da[2] = p20;
666 da[3] = p40;
667 da[74] = p10;
668 double a74 = a[74];
669 double a77 = a[77];
670 double a78 = a[78];
671 if (kle_eq == -1) {
672
673 a74 = -a74;
674 a77 = -a77;
675 a78 = -a78;
676 }
677 if (kle_eq == 0 ) {
678
679 a74 = semestrialCorrection(a74);
680 a77 = semestrialCorrection(a77);
681 a78 = semestrialCorrection(a78);
682 }
683 da[77] = p30;
684 da[78] = p50;
685 da[79] = p60;
686
687
688 fmfb[1] = f[1] - fbar[1];
689 fmfb[2] = f[2] - fbar[2];
690 fbm150[1] = fbar[1] - 150.0;
691 fbm150[2] = fbar[2];
692 da[4] = fmfb[1];
693 da[6] = fbm150[1];
694 da[4] = da[4] + a[70] * fmfb[2];
695 da[6] = da[6] + a[71] * fbm150[2];
696 da[70] = fmfb[2] * (a[4] + 2.0 * a[5] * da[4] + a[82] * p10 +
697 a[83] * p20 + a[84] * p30);
698 da[71] = fbm150[2] * (a[6] + 2.0 * a[69] * da[6] + a[85] * p10 +
699 a[86] * p20 + a[87] * p30);
700 da[5] = da[4] * da[4];
701 da[69] = da[6] * da[6];
702 da[82] = da[4] * p10;
703 da[83] = da[4] * p20;
704 da[84] = da[4] * p30;
705 da[85] = da[6] * p20;
706 da[86] = da[6] * p30;
707 da[87] = da[6] * p40;
708
709
710 final int ikp = 62;
711 final int ikpm = 67;
712 final double c2fi = 1.0 - p10mg * p10mg;
713 final double dkp = akp[1] + (a[ikp] + c2fi * a[ikp + 1]) * akp[2];
714 double dakp = a[7] + a[8] * p20mg + a[68] * p40mg +
715 2.0 * dkp * (a[60] + a[61] * p20mg +
716 a[75] * 2.0 * dkp * dkp);
717 da[ikp] = dakp * akp[2];
718 da[ikp + 1] = da[ikp] * c2fi;
719 final double dkpm = akp[3] + a[ikpm] * akp[4];
720 final double dakpm = a[64] + a[65] * p20mg + a[72] * p40mg +
721 2.0 * dkpm * (a[66] + a[73] * p20mg +
722 a[76] * 2.0 * dkpm * dkpm);
723 da[ikpm] = dakpm * akp[4];
724 da[7] = dkp;
725 da[8] = p20mg * dkp;
726 da[68] = p40mg * dkp;
727 da[60] = dkp * dkp;
728 da[61] = p20mg * da[60];
729 da[75] = da[60] * da[60];
730 da[64] = dkpm;
731 da[65] = p20mg * dkpm;
732 da[72] = p40mg * dkpm;
733 da[66] = dkpm * dkpm;
734 da[73] = p20mg * da[66];
735 da[76] = da[66] * da[66];
736
737
738 double f0 = a[4] * da[4] + a[5] * da[5] + a[6] * da[6] +
739 a[69] * da[69] + a[82] * da[82] + a[83] * da[83] +
740 a[84] * da[84] + a[85] * da[85] + a[86] * da[86] +
741 a[87] * da[87];
742 final double f1f = 1.0 + f0 * ff0;
743
744 f0 = f0 + a[2] * da[2] + a[3] * da[3] + a74 * da[74] +
745 a77 * da[77] + a[7] * da[7] + a[8] * da[8] +
746 a[60] * da[60] + a[61] * da[61] + a[68] * da[68] +
747 a[64] * da[64] + a[65] * da[65] + a[66] * da[66] +
748 a[72] * da[72] + a[73] * da[73] + a[75] * da[75] +
749 a[76] * da[76] + a78 * da[78] + a[79] * da[79];
750
751 da[9] = FastMath.cos(ROT * (day - a[11]));
752 da[10] = p20 * da[9];
753
754 da[12] = FastMath.cos(ROT2 * (day - a[14]));
755 da[13] = p20 * da[12];
756
757 final double coste = FastMath.cos(ROT * (day - a[18]));
758 da[15] = p10 * coste;
759 da[16] = p30 * coste;
760 da[17] = p50 * coste;
761
762 final double cos2te = FastMath.cos(ROT2 * (day - a[20]));
763 da[19] = p10 * cos2te;
764 da[39] = p30 * cos2te;
765 da[59] = p50 * cos2te;
766
767 da[21] = p11 * ch;
768 da[22] = p31 * ch;
769 da[23] = p51 * ch;
770 da[24] = da[21] * coste;
771 da[25] = p21 * ch * coste;
772 da[26] = p11 * sh;
773 da[27] = p31 * sh;
774 da[28] = p51 * sh;
775 da[29] = da[26] * coste;
776 da[30] = p21 * sh * coste;
777
778 da[31] = p22 * c2h;
779 da[37] = p42 * c2h;
780 da[32] = p32 * c2h * coste;
781 da[33] = p22 * s2h;
782 da[38] = p42 * s2h;
783 da[34] = p32 * s2h * coste;
784 da[88] = p32 * c2h;
785 da[89] = p32 * s2h;
786 da[90] = p52 * c2h;
787 da[91] = p52 * s2h;
788 double a88 = a[88];
789 double a89 = a[89];
790 double a90 = a[90];
791 double a91 = a[91];
792 if (kle_eq == -1) {
793 a88 = -a88;
794 a89 = -a89;
795 a90 = -a90;
796 a91 = -a91;
797 }
798 if (kle_eq == 0) {
799 a88 = semestrialCorrection(a88);
800 a89 = semestrialCorrection(a89);
801 a90 = semestrialCorrection(a90);
802 a91 = semestrialCorrection(a91);
803 }
804 da[92] = p62 * c2h;
805 da[93] = p62 * s2h;
806
807 da[35] = p33 * c3h;
808 da[36] = p33 * s3h;
809
810 double fp = a[9] * da[9] + a[10] * da[10] + a[12] * da[12] + a[13] * da[13] +
811 a[15] * da[15] + a[16] * da[16] + a[17] * da[17] + a[19] * da[19] +
812 a[21] * da[21] + a[22] * da[22] + a[23] * da[23] + a[24] * da[24] +
813 a[25] * da[25] + a[26] * da[26] + a[27] * da[27] + a[28] * da[28] +
814 a[29] * da[29] + a[30] * da[30] + a[31] * da[31] + a[32] * da[32] +
815 a[33] * da[33] + a[34] * da[34] + a[35] * da[35] + a[36] * da[36] +
816 a[37] * da[37] + a[38] * da[38] + a[39] * da[39] + a[59] * da[59] +
817 a88 * da[88] + a89 * da[89] + a90 * da[90] + a91 * da[91] +
818 a[92] * da[92] + a[93] * da[93];
819
820 da[40] = p10 * coste * dkp;
821 da[41] = p30 * coste * dkp;
822 da[42] = p50 * coste * dkp;
823 da[43] = p11 * ch * dkp;
824 da[44] = p31 * ch * dkp;
825 da[45] = p51 * ch * dkp;
826 da[46] = p11 * sh * dkp;
827 da[47] = p31 * sh * dkp;
828 da[48] = p51 * sh * dkp;
829
830
831 fp += a[40] * da[40] + a[41] * da[41] + a[42] * da[42] + a[43] * da[43] +
832 a[44] * da[44] + a[45] * da[45] + a[46] * da[46] + a[47] * da[47] +
833 a[48] * da[48];
834
835 dakp = (a[40] * p10 + a[41] * p30 + a[42] * p50) * coste +
836 (a[43] * p11 + a[44] * p31 + a[45] * p51) * ch +
837 (a[46] * p11 + a[47] * p31 + a[48] * p51) * sh;
838 da[ikp] += dakp * akp[2];
839 da[ikp + 1] = da[ikp] + dakp * c2fi * akp[2];
840
841 da[49] = p11 * clfl;
842 da[50] = p21 * clfl;
843 da[51] = p31 * clfl;
844 da[52] = p41 * clfl;
845 da[53] = p51 * clfl;
846 da[54] = p11 * slfl;
847 da[55] = p21 * slfl;
848 da[56] = p31 * slfl;
849 da[57] = p41 * slfl;
850 da[58] = p51 * slfl;
851
852
853 fp += a[49] * da[49] + a[50] * da[50] + a[51] * da[51] + a[52] * da[52] +
854 a[53] * da[53] + a[54] * da[54] + a[55] * da[55] + a[56] * da[56] +
855 a[57] * da[57] + a[58] * da[58];
856
857
858 return f0 + fp * f1f;
859
860 }
861
862
863
864
865
866
867 private double semestrialCorrection(final double param) {
868 final int debeq_pr = 59;
869 final int debeq_au = 244;
870 final double result;
871 if (day >= 100) {
872 final double xmult = (day - debeq_au) / 40.0;
873 result = param - 2.0 * param * xmult;
874 } else {
875 final double xmult = (day - debeq_pr) / 40.0;
876 result = 2.0 * param * xmult - param;
877 }
878 return result;
879 }
880
881
882 }
883
884
885
886
887 private static class FieldComputation<T extends CalculusFieldElement<T>> {
888
889
890 private int day;
891
892
893 private double[] f = new double[3];
894
895
896 private double[] fbar = new double[3];
897
898
899
900
901
902
903
904
905
906 private double[] akp = new double[5];
907
908
909 private final T clfl;
910
911
912 private final T slfl;
913
914
915 private final T ro;
916
917
918
919
920 private final T p10;
921 private final T p20;
922 private final T p30;
923 private final T p40;
924 private final T p50;
925 private final T p60;
926 private final T p11;
927 private final T p21;
928 private final T p31;
929 private final T p41;
930 private final T p51;
931 private final T p22;
932 private final T p32;
933 private final T p42;
934 private final T p52;
935 private final T p62;
936 private final T p33;
937 private final T p10mg;
938 private final T p20mg;
939 private final T p40mg;
940
941
942 private final T hl0;
943 private final T ch;
944 private final T sh;
945 private final T c2h;
946 private final T s2h;
947 private final T c3h;
948 private final T s3h;
949
950
951
952
953
954
955
956
957
958
959
960 FieldComputation(final int day,
961 final T altiKM, final T lon, final T lat,
962 final T hl, final double[] f, final double[] far,
963 final double[] akp) {
964
965 this.day = day;
966 this.f = f;
967 this.fbar = far;
968 this.akp = akp;
969
970
971 final FieldSinCos<T> scLat = FastMath.sinCos(lat);
972 final FieldSinCos<T> scLon = FastMath.sinCos(lon);
973
974
975 final T c = scLat.sin();
976 final T c2 = c.square();
977 final T c4 = c2.square();
978 final T s = scLat.cos();
979 final T s2 = s.multiply(s);
980 p10 = c;
981 p20 = c2.multiply(1.5).subtract(0.5);
982 p30 = c.multiply(c2.multiply(2.5).subtract(1.5));
983 p40 = c4.multiply(4.375).subtract(c2.multiply(3.75)).add(0.375);
984 p50 = c.multiply(c4.multiply(7.875).subtract(c2.multiply(8.75)).add(1.875));
985 p60 = (c.multiply(5.5).multiply(p50).subtract(p40.multiply(2.5))).divide(3.0);
986 p11 = s;
987 p21 = c.multiply(3.0).multiply(s);
988 p31 = s.multiply(c2.multiply(7.5).subtract(1.5));
989 p41 = c.multiply(s).multiply(c2.multiply(17.5).subtract(7.5));
990 p51 = s.multiply(c4.multiply(39.375).subtract(c2.multiply(26.25)).add(1.875));
991 p22 = s2.multiply(3.0);
992 p32 = c.multiply(15.0).multiply(s2);
993 p42 = s2.multiply(c2.multiply(52.5).subtract(7.5));
994 p52 = c.multiply(3.0).multiply(p42).subtract(p32.multiply(2.0));
995 p62 = c.multiply(2.75).multiply(p52).subtract(p42.multiply(1.75));
996 p33 = s.multiply(15.0).multiply(s2);
997
998
999 final T clmlmg = lon.subtract(XLMG).cos();
1000 final T cmg = s.multiply(CPMG).multiply(clmlmg).add(c.multiply(SPMG));
1001 final T cmg2 = cmg.multiply(cmg);
1002 final T cmg4 = cmg2.multiply(cmg2);
1003 p10mg = cmg;
1004 p20mg = cmg2.multiply(1.5).subtract(0.5);
1005 p40mg = cmg4.multiply(4.375).subtract(cmg2.multiply(3.75)).add(0.375);
1006
1007 clfl = scLon.cos();
1008 slfl = scLon.sin();
1009
1010
1011 hl0 = hl;
1012 final FieldSinCos<T> scHlo = FastMath.sinCos(hl0);
1013 ch = scHlo.cos();
1014 sh = scHlo.sin();
1015 c2h = ch.multiply(ch).subtract(sh.multiply(sh));
1016 s2h = ch.multiply(sh).multiply(2);
1017 c3h = c2h.multiply(ch).subtract(s2h.multiply(sh));
1018 s3h = s2h.multiply(ch).add(c2h.multiply(sh));
1019
1020 final double zlb = ZLB0;
1021
1022 final T[] dtt = MathArrays.buildArray(altiKM.getField(), tt.length);
1023 final T[] dh = MathArrays.buildArray(altiKM.getField(), tt.length);
1024 final T[] dhe = MathArrays.buildArray(altiKM.getField(), tt.length);
1025 final T[] dox = MathArrays.buildArray(altiKM.getField(), tt.length);
1026 final T[] daz2 = MathArrays.buildArray(altiKM.getField(), tt.length);
1027 final T[] do2 = MathArrays.buildArray(altiKM.getField(), tt.length);
1028 final T[] daz = MathArrays.buildArray(altiKM.getField(), tt.length);
1029 final T[] dt0 = MathArrays.buildArray(altiKM.getField(), tt.length);
1030 final T[] dtp = MathArrays.buildArray(altiKM.getField(), tt.length);
1031
1032
1033 int kleq = 1;
1034 final T gdelt = gFunction(tt, dtt, 1, kleq);
1035 dtt[1] = gdelt.add(1);
1036 final T tinf = dtt[1].multiply(tt[1]);
1037
1038 kleq = 0;
1039
1040 if (day < 59 || day > 284) {
1041 kleq = -1;
1042 }
1043 if (day > 99 && day < 244) {
1044 kleq = 1;
1045 }
1046
1047 final T gdelt0 = gFunction(t0, dt0, 0, kleq);
1048 dt0[1] = gdelt0.add(t0[1]).divide(t0[1]);
1049 final T t120 = gdelt0.add(t0[1]);
1050 final T gdeltp = gFunction(tp, dtp, 0, kleq);
1051 dtp[1] = gdeltp.add(tp[1]).divide(tp[1]);
1052 final T tp120 = gdeltp.add(tp[1]);
1053
1054
1055 final T sigma = tp120.divide(tinf.subtract(t120));
1056 final T dzeta = altiKM.add(RE).reciprocal().multiply(zlb + RE);
1057 final T zeta = altiKM.subtract(zlb).multiply(dzeta);
1058 final T sigzeta = sigma.multiply(zeta);
1059 final T expsz = sigzeta.negate().exp();
1060 final T tz = tinf.subtract(tinf.subtract(t120).multiply(expsz));
1061
1062 final T[] dbase = MathArrays.buildArray(altiKM.getField(), 7);
1063
1064 kleq = 1;
1065
1066 final T gdelh = gFunction(h, dh, 0, kleq);
1067 dh[1] = gdelh.exp();
1068 dbase[1] = dh[1].multiply(h[1]);
1069
1070 final T gdelhe = gFunction(he, dhe, 0, kleq);
1071 dhe[1] = gdelhe.exp();
1072 dbase[2] = dhe[1].multiply(he[1]);
1073
1074 final T gdelo = gFunction(o, dox, 1, kleq);
1075 dox[1] = gdelo.exp();
1076 dbase[3] = dox[1].multiply(o[1]);
1077
1078 final T gdelaz2 = gFunction(az2, daz2, 1, kleq);
1079 daz2[1] = gdelaz2.exp();
1080 dbase[4] = daz2[1].multiply(az2[1]);
1081
1082 final T gdelo2 = gFunction(o2, do2, 1, kleq);
1083 do2[1] = gdelo2.exp();
1084 dbase[5] = do2[1].multiply(o2[1]);
1085
1086 final T gdelaz = gFunction(az, daz, 1, kleq);
1087 daz[1] = gdelaz.exp();
1088 dbase[6] = daz[1].multiply(az[1]);
1089
1090 final double zlbre = 1.0 + zlb / RE;
1091 final T glb = sigma.multiply(RGAS).multiply(tinf).reciprocal().multiply(GSURF / (zlbre * zlbre));
1092 final T t120tz = t120.divide(tz);
1093
1094
1095
1096
1097
1098
1099
1100
1101 T tmpro = altiKM.getField().getZero();
1102 for (int i = 1; i <= 6; i++) {
1103 final T gamma = glb.multiply(MA[i]);
1104 final T upapg = gamma.add(1.0 + ALEFA[i]);
1105 final T fzI = (t120tz.log().multiply(upapg).subtract(sigzeta.multiply(gamma))).exp();
1106
1107 final T ccI = dbase[i].multiply(fzI);
1108
1109 tmpro = tmpro.add(ccI.multiply(VMA[i]));
1110 }
1111 this.ro = tmpro;
1112
1113 }
1114
1115
1116
1117
1118
1119
1120
1121
1122 private T gFunction(final double[] a, final T[] da,
1123 final int ff0, final int kle_eq) {
1124
1125 final T zero = da[0].getField().getZero();
1126 final double[] fmfb = new double[3];
1127 final double[] fbm150 = new double[3];
1128
1129
1130 da[2] = p20;
1131 da[3] = p40;
1132 da[74] = p10;
1133 double a74 = a[74];
1134 double a77 = a[77];
1135 double a78 = a[78];
1136 if (kle_eq == -1) {
1137
1138 a74 = -a74;
1139 a77 = -a77;
1140 a78 = -a78;
1141 }
1142 if (kle_eq == 0 ) {
1143
1144 a74 = semestrialCorrection(a74);
1145 a77 = semestrialCorrection(a77);
1146 a78 = semestrialCorrection(a78);
1147 }
1148 da[77] = p30;
1149 da[78] = p50;
1150 da[79] = p60;
1151
1152
1153 fmfb[1] = f[1] - fbar[1];
1154 fmfb[2] = f[2] - fbar[2];
1155 fbm150[1] = fbar[1] - 150.0;
1156 fbm150[2] = fbar[2];
1157 da[4] = zero.newInstance(fmfb[1]);
1158 da[6] = zero.newInstance(fbm150[1]);
1159 da[4] = da[4].add(a[70] * fmfb[2]);
1160 da[6] = da[6].add(a[71] * fbm150[2]);
1161 da[70] = da[4].multiply(a[ 5]).multiply(2).
1162 add(p10.multiply(a[82])).
1163 add(p20.multiply(a[83])).
1164 add(p30.multiply(a[84])).
1165 add(a[4]).
1166 multiply(fmfb[2]);
1167 da[71] = da[6].multiply(a[69]).multiply(2).
1168 add(p10.multiply(a[85])).
1169 add(p20.multiply(a[86])).
1170 add(p30.multiply(a[87])).
1171 add(a[6]).
1172 multiply(fbm150[2]);
1173 da[5] = da[4].multiply(da[4]);
1174 da[69] = da[6].multiply(da[6]);
1175 da[82] = da[4].multiply(p10);
1176 da[83] = da[4].multiply(p20);
1177 da[84] = da[4].multiply(p30);
1178 da[85] = da[6].multiply(p20);
1179 da[86] = da[6].multiply(p30);
1180 da[87] = da[6].multiply(p40);
1181
1182
1183 final int ikp = 62;
1184 final int ikpm = 67;
1185 final T c2fi = p10mg.multiply(p10mg).negate().add(1);
1186 final T dkp = c2fi.multiply(a[ikp + 1]).add(a[ikp]).multiply(akp[2]).add(akp[1]);
1187 T dakp = p20mg.multiply(a[8]).add(p40mg.multiply(a[68])).
1188 add(p20mg.multiply(a[61]).add(dkp.square().multiply(2 * a[75]).add(a[60])).multiply(dkp.multiply(2))).
1189 add(a[7]);
1190 da[ikp] = dakp.multiply(akp[2]);
1191 da[ikp + 1] = da[ikp].multiply(c2fi);
1192 final double dkpm = akp[3] + a[ikpm] * akp[4];
1193 final T dakpm = p20mg.multiply(a[65]).add(p40mg.multiply(a[72])).
1194 add(p20mg.multiply(a[73]).add(a[66] + a[76] * 2.0 * dkpm * dkpm).multiply( 2.0 * dkpm)).
1195 add(a[64]);
1196 da[ikpm] = dakpm.multiply(akp[4]);
1197 da[7] = dkp;
1198 da[8] = p20mg.multiply(dkp);
1199 da[68] = p40mg.multiply(dkp);
1200 da[60] = dkp.square();
1201 da[61] = p20mg.multiply(da[60]);
1202 da[75] = da[60].square();
1203 da[64] = zero.newInstance(dkpm);
1204 da[65] = p20mg.multiply(dkpm);
1205 da[72] = p40mg.multiply(dkpm);
1206 da[66] = zero.newInstance(dkpm * dkpm);
1207 da[73] = p20mg.multiply(da[66]);
1208 da[76] = da[66].square();
1209
1210
1211 T f0 = da[4].multiply(a[4]).
1212 add(da[5].multiply(a[5])).
1213 add(da[6].multiply(a[6])).
1214 add(da[69].multiply(a[69])).
1215 add(da[82].multiply(a[82])).
1216 add(da[83].multiply(a[83])).
1217 add(da[84].multiply(a[84])).
1218 add(da[85].multiply(a[85])).
1219 add(da[86].multiply(a[86])).
1220 add(da[87].multiply(a[87]));
1221 final T f1f = f0.multiply(ff0).add(1);
1222
1223 f0 = f0.
1224 add(da[2].multiply(a[2])).
1225 add(da[3].multiply(a[3])).
1226 add(da[7].multiply(a[7])).
1227 add(da[8].multiply(a[8])).
1228 add(da[60].multiply(a[60])).
1229 add(da[61].multiply(a[61])).
1230 add(da[68].multiply(a[68])).
1231 add(da[64].multiply(a[64])).
1232 add(da[65].multiply(a[65])).
1233 add(da[66].multiply(a[66])).
1234 add(da[72].multiply(a[72])).
1235 add(da[73].multiply(a[73])).
1236 add(da[74].multiply(a74)).
1237 add(da[75].multiply(a[75])).
1238 add(da[76].multiply(a[76])).
1239 add(da[77].multiply(a77)).
1240 add(da[78].multiply(a78)).
1241 add(da[79].multiply(a[79]));
1242
1243 da[9] = zero.newInstance(FastMath.cos(ROT * (day - a[11])));
1244 da[10] = p20.multiply(da[9]);
1245
1246 da[12] = zero.newInstance(FastMath.cos(ROT2 * (day - a[14])));
1247 da[13] = p20.multiply(da[12]);
1248
1249 final double coste = FastMath.cos(ROT * (day - a[18]));
1250 da[15] = p10.multiply(coste);
1251 da[16] = p30.multiply(coste);
1252 da[17] = p50.multiply(coste);
1253
1254 final double cos2te = FastMath.cos(ROT2 * (day - a[20]));
1255 da[19] = p10.multiply(cos2te);
1256 da[39] = p30.multiply(cos2te);
1257 da[59] = p50.multiply(cos2te);
1258
1259 da[21] = p11.multiply(ch);
1260 da[22] = p31.multiply(ch);
1261 da[23] = p51.multiply(ch);
1262 da[24] = da[21].multiply(coste);
1263 da[25] = p21.multiply(ch).multiply(coste);
1264 da[26] = p11.multiply(sh);
1265 da[27] = p31.multiply(sh);
1266 da[28] = p51.multiply(sh);
1267 da[29] = da[26].multiply(coste);
1268 da[30] = p21.multiply(sh).multiply(coste);
1269
1270 da[31] = p22.multiply(c2h);
1271 da[37] = p42.multiply(c2h);
1272 da[32] = p32.multiply(c2h).multiply(coste);
1273 da[33] = p22.multiply(s2h);
1274 da[38] = p42.multiply(s2h);
1275 da[34] = p32.multiply(s2h).multiply(coste);
1276 da[88] = p32.multiply(c2h);
1277 da[89] = p32.multiply(s2h);
1278 da[90] = p52.multiply(c2h);
1279 da[91] = p52.multiply(s2h);
1280 double a88 = a[88];
1281 double a89 = a[89];
1282 double a90 = a[90];
1283 double a91 = a[91];
1284 if (kle_eq == -1) {
1285 a88 = -a88;
1286 a89 = -a89;
1287 a90 = -a90;
1288 a91 = -a91;
1289 }
1290 if (kle_eq == 0) {
1291 a88 = semestrialCorrection(a88);
1292 a89 = semestrialCorrection(a89);
1293 a90 = semestrialCorrection(a90);
1294 a91 = semestrialCorrection(a91);
1295 }
1296 da[92] = p62.multiply(c2h);
1297 da[93] = p62.multiply(s2h);
1298
1299 da[35] = p33.multiply(c3h);
1300 da[36] = p33.multiply(s3h);
1301
1302 T fp = da[ 9].multiply(a[ 9]) .add(da[10].multiply(a[10])).add(da[12].multiply(a[12])).add(da[13].multiply(a[13])).
1303 add(da[15].multiply(a[15])).add(da[16].multiply(a[16])).add(da[17].multiply(a[17])).add(da[19].multiply(a[19])).
1304 add(da[21].multiply(a[21])).add(da[22].multiply(a[22])).add(da[23].multiply(a[23])).add(da[24].multiply(a[24])).
1305 add(da[25].multiply(a[25])).add(da[26].multiply(a[26])).add(da[27].multiply(a[27])).add(da[28].multiply(a[28])).
1306 add(da[29].multiply(a[29])).add(da[30].multiply(a[30])).add(da[31].multiply(a[31])).add(da[32].multiply(a[32])).
1307 add(da[33].multiply(a[33])).add(da[34].multiply(a[34])).add(da[35].multiply(a[35])).add(da[36].multiply(a[36])).
1308 add(da[37].multiply(a[37])).add(da[38].multiply(a[38])).add(da[39].multiply(a[39])).add(da[59].multiply(a[59])).
1309 add(da[88].multiply(a88)) .add(da[89].multiply(a89) ).add(da[90].multiply(a90) ).add(da[91].multiply(a91) ).
1310 add(da[92].multiply(a[92])).add(da[93].multiply(a[93]));
1311
1312 da[40] = p10.multiply(coste).multiply(dkp);
1313 da[41] = p30.multiply(coste).multiply(dkp);
1314 da[42] = p50.multiply(coste).multiply(dkp);
1315 da[43] = p11.multiply(ch).multiply(dkp);
1316 da[44] = p31.multiply(ch).multiply(dkp);
1317 da[45] = p51.multiply(ch).multiply(dkp);
1318 da[46] = p11.multiply(sh).multiply(dkp);
1319 da[47] = p31.multiply(sh).multiply(dkp);
1320 da[48] = p51.multiply(sh).multiply(dkp);
1321
1322
1323 fp = fp.
1324 add(da[40].multiply(a[40])).
1325 add(da[41].multiply(a[41])).
1326 add(da[42].multiply(a[42])).
1327 add(da[43].multiply(a[43])).
1328 add(da[44].multiply(a[44])).
1329 add(da[45].multiply(a[45])).
1330 add(da[46].multiply(a[46])).
1331 add(da[47].multiply(a[47])).
1332 add(da[48].multiply(a[48]));
1333
1334 dakp = p10.multiply(a[40]).add(p30.multiply(a[41])).add(p50.multiply(a[42])).multiply(coste).
1335 add(p11.multiply(a[40]).add(p31.multiply(a[44])).add(p51.multiply(a[45])).multiply(ch)).
1336 add(p11.multiply(a[46]).add(p31.multiply(a[47])).add(p51.multiply(a[48])).multiply(sh));
1337 da[ikp] = da[ikp].add(dakp.multiply(akp[2]));
1338 da[ikp + 1] = da[ikp].add(dakp.multiply(c2fi).multiply(akp[2]));
1339
1340 da[49] = p11.multiply(clfl);
1341 da[50] = p21.multiply(clfl);
1342 da[51] = p31.multiply(clfl);
1343 da[52] = p41.multiply(clfl);
1344 da[53] = p51.multiply(clfl);
1345 da[54] = p11.multiply(slfl);
1346 da[55] = p21.multiply(slfl);
1347 da[56] = p31.multiply(slfl);
1348 da[57] = p41.multiply(slfl);
1349 da[58] = p51.multiply(slfl);
1350
1351
1352 fp = fp.
1353 add(da[49].multiply(a[49])).
1354 add(da[50].multiply(a[50])).
1355 add(da[51].multiply(a[51])).
1356 add(da[52].multiply(a[52])).
1357 add(da[53].multiply(a[53])).
1358 add(da[54].multiply(a[54])).
1359 add(da[55].multiply(a[55])).
1360 add(da[56].multiply(a[56])).
1361 add(da[57].multiply(a[57])).
1362 add(da[58].multiply(a[58]));
1363
1364
1365 return f0.add(fp.multiply(f1f));
1366
1367 }
1368
1369
1370
1371
1372
1373
1374 private double semestrialCorrection(final double param) {
1375 final int debeq_pr = 59;
1376 final int debeq_au = 244;
1377 final double result;
1378 if (day >= 100) {
1379 final double xmult = (day - debeq_au) / 40.0;
1380 result = param - 2.0 * param * xmult;
1381 } else {
1382 final double xmult = (day - debeq_pr) / 40.0;
1383 result = 2.0 * param * xmult - param;
1384 }
1385 return result;
1386 }
1387
1388 }
1389
1390 }