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