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.Field;
21 import org.hipparchus.exception.LocalizedCoreFormats;
22 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
23 import org.hipparchus.geometry.euclidean.threed.Vector3D;
24 import org.hipparchus.util.FastMath;
25 import org.hipparchus.util.FieldSinCos;
26 import org.hipparchus.util.MathArrays;
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.TimeComponents;
40 import org.orekit.time.TimeScale;
41 import org.orekit.utils.IERSConventions;
42 import org.orekit.utils.ExtendedPositionProvider;
43
44 import java.util.Arrays;
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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128 public class NRLMSISE00 extends AbstractSunInfluencedAtmosphere {
129
130
131
132 private static final int HELIUM = 0;
133
134
135 private static final int ATOMIC_OXYGEN = 1;
136
137
138 private static final int MOLECULAR_NITROGEN = 2;
139
140
141 private static final int MOLECULAR_OXYGEN = 3;
142
143
144 private static final int ARGON = 4;
145
146
147 private static final int TOTAL_MASS = 5;
148
149
150 private static final int HYDROGEN = 6;
151
152
153 private static final int ATOMIC_NITROGEN = 7;
154
155
156 private static final int ANOMALOUS_OXYGEN = 8;
157
158
159 private static final int EXOSPHERIC = 0;
160
161
162 private static final int ALTITUDE = 1;
163
164
165
166
167 private static final double DEG_TO_RAD = 1.74533e-2;
168
169
170 private static final double DAY_TO_RAD = 1.72142e-2;
171
172
173 private static final double HOUR_TO_RAD = 0.2618;
174
175
176 private static final double SEC_TO_RAD = 7.2722e-5;
177
178
179
180
181 private static final double LAT_REF = 45.;
182
183
184 private static final double G_REF = 980.616;
185
186
187
188
189 private static final double AMU = 1.66e-27;
190
191
192 private static final double R_GAS = 831.4;
193
194
195 private static final double H_MASS = 1.;
196
197
198 private static final double HE_MASS = 4.;
199
200
201 private static final double N_MASS = 14.;
202
203
204 private static final double N2_MASS = 2. * N_MASS;
205
206
207 private static final double O_MASS = 16.;
208
209
210 private static final double O2_MASS = 2. * O_MASS;
211
212
213 private static final double AR_MASS = 40.;
214
215
216
217
218 private static final double FLUX_REF = 150.;
219
220
221 private static final double[] ZN1 = {123.435, 110.0, 100.0, 90.0, 72.5};
222
223
224 private static final double[] ZN2 = {72.5, 55.0, 45.0, 32.5};
225
226
227 private static final double[] ZN3 = {32.5, 20.0, 15.0, 10.0, 0.0};
228
229
230 private static final double ZMIX = 62.5;
231
232
233 private static final double[] PT = {
234 9.86573e-01, 1.62228e-02, 1.55270e-02, -1.04323e-01, -3.75801e-03,
235 -1.18538e-03, -1.24043e-01, 4.56820e-03, 8.76018e-03, -1.36235e-01,
236 -3.52427e-02, 8.84181e-03, -5.92127e-03, -8.61650e+00, 0.00000e+00,
237 1.28492e-02, 0.00000e+00, 1.30096e+02, 1.04567e-02, 1.65686e-03,
238 -5.53887e-06, 2.97810e-03, 0.00000e+00, 5.13122e-03, 8.66784e-02,
239 1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, -7.27026e-06,
240 0.00000e+00, 6.74494e+00, 4.93933e-03, 2.21656e-03, 2.50802e-03,
241 0.00000e+00, 0.00000e+00, -2.08841e-02, -1.79873e+00, 1.45103e-03,
242 2.81769e-04, -1.44703e-03, -5.16394e-05, 8.47001e-02, 1.70147e-01,
243 5.72562e-03, 5.07493e-05, 4.36148e-03, 1.17863e-04, 4.74364e-03,
244 6.61278e-03, 4.34292e-05, 1.44373e-03, 2.41470e-05, 2.84426e-03,
245 8.56560e-04, 2.04028e-03, 0.00000e+00, -3.15994e+03, -2.46423e-03,
246 1.13843e-03, 4.20512e-04, 0.00000e+00, -9.77214e+01, 6.77794e-03,
247 5.27499e-03, 1.14936e-03, 0.00000e+00, -6.61311e-03, -1.84255e-02,
248 -1.96259e-02, 2.98618e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
249 6.44574e+02, 8.84668e-04, 5.05066e-04, 0.00000e+00, 4.02881e+03,
250 -1.89503e-03, 0.00000e+00, 0.00000e+00, 8.21407e-04, 2.06780e-03,
251 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
252 -1.20410e-02, -3.63963e-03, 9.92070e-05, -1.15284e-04, -6.33059e-05,
253 -6.05545e-01, 8.34218e-03, -9.13036e+01, 3.71042e-04, 0.00000e+00,
254 4.19000e-04, 2.70928e-03, 3.31507e-03, -4.44508e-03, -4.96334e-03,
255 -1.60449e-03, 3.95119e-03, 2.48924e-03, 5.09815e-04, 4.05302e-03,
256 2.24076e-03, 0.00000e+00, 6.84256e-03, 4.66354e-04, 0.00000e+00,
257 -3.68328e-04, 0.00000e+00, 0.00000e+00, -1.46870e+02, 0.00000e+00,
258 0.00000e+00, 1.09501e-03, 4.65156e-04, 5.62583e-04, 3.21596e+00,
259 6.43168e-04, 3.14860e-03, 3.40738e-03, 1.78481e-03, 9.62532e-04,
260 5.58171e-04, 3.43731e+00, -2.33195e-01, 5.10289e-04, 0.00000e+00,
261 0.00000e+00, -9.25347e+04, 0.00000e+00, -1.99639e-03, 0.00000e+00,
262 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
263 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
264 };
265
266
267 private static final double[][] PD = {
268
269 {
270 1.09979e+00, -4.88060e-02, -1.97501e-01, -9.10280e-02, -6.96558e-03,
271 2.42136e-02, 3.91333e-01, -7.20068e-03, -3.22718e-02, 1.41508e+00,
272 1.68194e-01, 1.85282e-02, 1.09384e-01, -7.24282e+00, 0.00000e+00,
273 2.96377e-01, -4.97210e-02, 1.04114e+02, -8.61108e-02, -7.29177e-04,
274 1.48998e-06, 1.08629e-03, 0.00000e+00, 0.00000e+00, 8.31090e-02,
275 1.12818e-01, -5.75005e-02, -1.29919e-02, -1.78849e-02, -2.86343e-06,
276 0.00000e+00, -1.51187e+02, -6.65902e-03, 0.00000e+00, -2.02069e-03,
277 0.00000e+00, 0.00000e+00, 4.32264e-02, -2.80444e+01, -3.26789e-03,
278 2.47461e-03, 0.00000e+00, 0.00000e+00, 9.82100e-02, 1.22714e-01,
279 -3.96450e-02, 0.00000e+00, -2.76489e-03, 0.00000e+00, 1.87723e-03,
280 -8.09813e-03, 4.34428e-05, -7.70932e-03, 0.00000e+00, -2.28894e-03,
281 -5.69070e-03, -5.22193e-03, 6.00692e-03, -7.80434e+03, -3.48336e-03,
282 -6.38362e-03, -1.82190e-03, 0.00000e+00, -7.58976e+01, -2.17875e-02,
283 -1.72524e-02, -9.06287e-03, 0.00000e+00, 2.44725e-02, 8.66040e-02,
284 1.05712e-01, 3.02543e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
285 -6.01364e+03, -5.64668e-03, -2.54157e-03, 0.00000e+00, 3.15611e+02,
286 -5.69158e-03, 0.00000e+00, 0.00000e+00, -4.47216e-03, -4.49523e-03,
287 4.64428e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
288 4.51236e-02, 2.46520e-02, 6.17794e-03, 0.00000e+00, 0.00000e+00,
289 -3.62944e-01, -4.80022e-02, -7.57230e+01, -1.99656e-03, 0.00000e+00,
290 -5.18780e-03, -1.73990e-02, -9.03485e-03, 7.48465e-03, 1.53267e-02,
291 1.06296e-02, 1.18655e-02, 2.55569e-03, 1.69020e-03, 3.51936e-02,
292 -1.81242e-02, 0.00000e+00, -1.00529e-01, -5.10574e-03, 0.00000e+00,
293 2.10228e-03, 0.00000e+00, 0.00000e+00, -1.73255e+02, 5.07833e-01,
294 -2.41408e-01, 8.75414e-03, 2.77527e-03, -8.90353e-05, -5.25148e+00,
295 -5.83899e-03, -2.09122e-02, -9.63530e-03, 9.77164e-03, 4.07051e-03,
296 2.53555e-04, -5.52875e+00, -3.55993e-01, -2.49231e-03, 0.00000e+00,
297 0.00000e+00, 2.86026e+01, 0.00000e+00, 3.42722e-04, 0.00000e+00,
298 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
299 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
300 },
301
302 {
303 1.02315e+00, -1.59710e-01, -1.06630e-01, -1.77074e-02, -4.42726e-03,
304 3.44803e-02, 4.45613e-02, -3.33751e-02, -5.73598e-02, 3.50360e-01,
305 6.33053e-02, 2.16221e-02, 5.42577e-02, -5.74193e+00, 0.00000e+00,
306 1.90891e-01, -1.39194e-02, 1.01102e+02, 8.16363e-02, 1.33717e-04,
307 6.54403e-06, 3.10295e-03, 0.00000e+00, 0.00000e+00, 5.38205e-02,
308 1.23910e-01, -1.39831e-02, 0.00000e+00, 0.00000e+00, -3.95915e-06,
309 0.00000e+00, -7.14651e-01, -5.01027e-03, 0.00000e+00, -3.24756e-03,
310 0.00000e+00, 0.00000e+00, 4.42173e-02, -1.31598e+01, -3.15626e-03,
311 1.24574e-03, -1.47626e-03, -1.55461e-03, 6.40682e-02, 1.34898e-01,
312 -2.42415e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 6.13666e-04,
313 -5.40373e-03, 2.61635e-05, -3.33012e-03, 0.00000e+00, -3.08101e-03,
314 -2.42679e-03, -3.36086e-03, 0.00000e+00, -1.18979e+03, -5.04738e-02,
315 -2.61547e-03, -1.03132e-03, 1.91583e-04, -8.38132e+01, -1.40517e-02,
316 -1.14167e-02, -4.08012e-03, 1.73522e-04, -1.39644e-02, -6.64128e-02,
317 -6.85152e-02, -1.34414e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
318 6.07916e+02, -4.12220e-03, -2.20996e-03, 0.00000e+00, 1.70277e+03,
319 -4.63015e-03, 0.00000e+00, 0.00000e+00, -2.25360e-03, -2.96204e-03,
320 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
321 3.92786e-02, 1.31186e-02, -1.78086e-03, 0.00000e+00, 0.00000e+00,
322 -3.90083e-01, -2.84741e-02, -7.78400e+01, -1.02601e-03, 0.00000e+00,
323 -7.26485e-04, -5.42181e-03, -5.59305e-03, 1.22825e-02, 1.23868e-02,
324 6.68835e-03, -1.03303e-02, -9.51903e-03, 2.70021e-04, -2.57084e-02,
325 -1.32430e-02, 0.00000e+00, -3.81000e-02, -3.16810e-03, 0.00000e+00,
326 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
327 0.00000e+00, -9.05762e-04, -2.14590e-03, -1.17824e-03, 3.66732e+00,
328 -3.79729e-04, -6.13966e-03, -5.09082e-03, -1.96332e-03, -3.08280e-03,
329 -9.75222e-04, 4.03315e+00, -2.52710e-01, 0.00000e+00, 0.00000e+00,
330 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
331 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
332 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
333 },
334
335 {
336 1.16112e+00, 0.00000e+00, 0.00000e+00, 3.33725e-02, 0.00000e+00,
337 3.48637e-02, -5.44368e-03, 0.00000e+00, -6.73940e-02, 1.74754e-01,
338 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.74712e+02, 0.00000e+00,
339 1.26733e-01, 0.00000e+00, 1.03154e+02, 5.52075e-02, 0.00000e+00,
340 0.00000e+00, 8.13525e-04, 0.00000e+00, 0.00000e+00, 8.66784e-02,
341 1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
342 0.00000e+00, -2.50482e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
343 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.48894e-03,
344 6.16053e-04, -5.79716e-04, 2.95482e-03, 8.47001e-02, 1.70147e-01,
345 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
346 0.00000e+00, 2.47425e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
347 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
348 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
349 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
350 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
351 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
352 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
353 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
354 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
355 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
356 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
357 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
358 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
359 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
360 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
361 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
362 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
363 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
364 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
365 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
366 },
367
368 {
369 9.44846e-01, 0.00000e+00, 0.00000e+00, -3.08617e-02, 0.00000e+00,
370 -2.44019e-02, 6.48607e-03, 0.00000e+00, 3.08181e-02, 4.59392e-02,
371 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.74712e+02, 0.00000e+00,
372 2.13260e-02, 0.00000e+00, -3.56958e+02, 0.00000e+00, 1.82278e-04,
373 0.00000e+00, 3.07472e-04, 0.00000e+00, 0.00000e+00, 8.66784e-02,
374 1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
375 0.00000e+00, 0.00000e+00, 3.83054e-03, 0.00000e+00, 0.00000e+00,
376 -1.93065e-03, -1.45090e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
377 0.00000e+00, -1.23493e-03, 1.36736e-03, 8.47001e-02, 1.70147e-01,
378 3.71469e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
379 5.10250e-03, 2.47425e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
380 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
381 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
382 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
383 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
384 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
385 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
386 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
387 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
388 0.00000e+00, 3.68756e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
389 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
390 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
391 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
392 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
393 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
394 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
395 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
396 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
397 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
398 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
399 },
400
401 {
402 1.35580e+00, 1.44816e-01, 0.00000e+00, 6.07767e-02, 0.00000e+00,
403 2.94777e-02, 7.46900e-02, 0.00000e+00, -9.23822e-02, 8.57342e-02,
404 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.38636e+01, 0.00000e+00,
405 7.71653e-02, 0.00000e+00, 8.18751e+01, 1.87736e-02, 0.00000e+00,
406 0.00000e+00, 1.49667e-02, 0.00000e+00, 0.00000e+00, 8.66784e-02,
407 1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
408 0.00000e+00, -3.67874e+02, 5.48158e-03, 0.00000e+00, 0.00000e+00,
409 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
410 0.00000e+00, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,
411 1.22631e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
412 8.17187e-03, 3.71617e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
413 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
414 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.10826e-03,
415 -3.13640e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
416 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
417 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
418 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
419 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
420 -7.35742e-02, -5.00266e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
421 0.00000e+00, 1.94965e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
422 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
423 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
424 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
425 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
426 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
427 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
428 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
429 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
430 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
431 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
432 },
433
434 {
435 1.04761e+00, 2.00165e-01, 2.37697e-01, 3.68552e-02, 0.00000e+00,
436 3.57202e-02, -2.14075e-01, 0.00000e+00, -1.08018e-01, -3.73981e-01,
437 0.00000e+00, 3.10022e-02, -1.16305e-03, -2.07596e+01, 0.00000e+00,
438 8.64502e-02, 0.00000e+00, 9.74908e+01, 5.16707e-02, 0.00000e+00,
439 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 8.66784e-02,
440 1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
441 0.00000e+00, 3.46193e+02, 1.34297e-02, 0.00000e+00, 0.00000e+00,
442 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -3.48509e-03,
443 -1.54689e-04, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,
444 1.47753e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
445 1.89320e-02, 3.68181e-05, 1.32570e-02, 0.00000e+00, 0.00000e+00,
446 3.59719e-03, 7.44328e-03, -1.00023e-03, -6.50528e+03, 0.00000e+00,
447 1.03485e-02, -1.00983e-03, -4.06916e-03, -6.60864e+01, -1.71533e-02,
448 1.10605e-02, 1.20300e-02, -5.20034e-03, 0.00000e+00, 0.00000e+00,
449 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
450 -2.62769e+03, 7.13755e-03, 4.17999e-03, 0.00000e+00, 1.25910e+04,
451 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.23595e-03, 4.60217e-03,
452 5.71794e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
453 -3.18353e-02, -2.35526e-02, -1.36189e-02, 0.00000e+00, 0.00000e+00,
454 0.00000e+00, 2.03522e-02, -6.67837e+01, -1.09724e-03, 0.00000e+00,
455 -1.38821e-02, 1.60468e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
456 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.51574e-02,
457 -5.44470e-04, 0.00000e+00, 7.28224e-02, 6.59413e-02, 0.00000e+00,
458 -5.15692e-03, 0.00000e+00, 0.00000e+00, -3.70367e+03, 0.00000e+00,
459 0.00000e+00, 1.36131e-02, 5.38153e-03, 0.00000e+00, 4.76285e+00,
460 -1.75677e-02, 2.26301e-02, 0.00000e+00, 1.76631e-02, 4.77162e-03,
461 0.00000e+00, 5.39354e+00, 0.00000e+00, -7.51710e-03, 0.00000e+00,
462 0.00000e+00, -8.82736e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
463 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
464 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
465 },
466
467 {
468 1.26376e+00, -2.14304e-01, -1.49984e-01, 2.30404e-01, 2.98237e-02,
469 2.68673e-02, 2.96228e-01, 2.21900e-02, -2.07655e-02, 4.52506e-01,
470 1.20105e-01, 3.24420e-02, 4.24816e-02, -9.14313e+00, 0.00000e+00,
471 2.47178e-02, -2.88229e-02, 8.12805e+01, 5.10380e-02, -5.80611e-03,
472 2.51236e-05, -1.24083e-02, 0.00000e+00, 0.00000e+00, 8.66784e-02,
473 1.58727e-01, -3.48190e-02, 0.00000e+00, 0.00000e+00, 2.89885e-05,
474 0.00000e+00, 1.53595e+02, -1.68604e-02, 0.00000e+00, 1.01015e-02,
475 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.84552e-04,
476 -1.22181e-03, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,
477 -1.04927e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, -5.91313e-03,
478 -2.30501e-02, 3.14758e-05, 0.00000e+00, 0.00000e+00, 1.26956e-02,
479 8.35489e-03, 3.10513e-04, 0.00000e+00, 3.42119e+03, -2.45017e-03,
480 -4.27154e-04, 5.45152e-04, 1.89896e-03, 2.89121e+01, -6.49973e-03,
481 -1.93855e-02, -1.48492e-02, 0.00000e+00, -5.10576e-02, 7.87306e-02,
482 9.51981e-02, -1.49422e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
483 2.65503e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
484 0.00000e+00, 0.00000e+00, 0.00000e+00, 6.37110e-03, 3.24789e-04,
485 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
486 6.14274e-02, 1.00376e-02, -8.41083e-04, 0.00000e+00, 0.00000e+00,
487 0.00000e+00, -1.27099e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
488 -3.94077e-03, -1.28601e-02, -7.97616e-03, 0.00000e+00, 0.00000e+00,
489 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
490 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
491 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
492 0.00000e+00, -6.71465e-03, -1.69799e-03, 1.93772e-03, 3.81140e+00,
493 -7.79290e-03, -1.82589e-02, -1.25860e-02, -1.04311e-02, -3.02465e-03,
494 2.43063e-03, 3.63237e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
495 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
496 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
497 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
498 },
499
500 {
501 7.09557e+01, -3.26740e-01, 0.00000e+00, -5.16829e-01, -1.71664e-03,
502 9.09310e-02, -6.71500e-01, -1.47771e-01, -9.27471e-02, -2.30862e-01,
503 -1.56410e-01, 1.34455e-02, -1.19717e-01, 2.52151e+00, 0.00000e+00,
504 -2.41582e-01, 5.92939e-02, 4.39756e+00, 9.15280e-02, 4.41292e-03,
505 0.00000e+00, 8.66807e-03, 0.00000e+00, 0.00000e+00, 8.66784e-02,
506 1.58727e-01, 9.74701e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
507 0.00000e+00, 6.70217e+01, -1.31660e-03, 0.00000e+00, -1.65317e-02,
508 0.00000e+00, 0.00000e+00, 8.50247e-02, 2.77428e+01, 4.98658e-03,
509 6.15115e-03, 9.50156e-03, -2.12723e-02, 8.47001e-02, 1.70147e-01,
510 -2.38645e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.37380e-03,
511 -8.41918e-03, 2.80145e-05, 7.12383e-03, 0.00000e+00, -1.66209e-02,
512 1.03533e-04, -1.68898e-02, 0.00000e+00, 3.64526e+03, 0.00000e+00,
513 6.54077e-03, 3.69130e-04, 9.94419e-04, 8.42803e+01, -1.16124e-02,
514 -7.74414e-03, -1.68844e-03, 1.42809e-03, -1.92955e-03, 1.17225e-01,
515 -2.41512e-02, 1.50521e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
516 1.60261e+03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
517 0.00000e+00, 0.00000e+00, 0.00000e+00, -3.54403e-04, -1.87270e-02,
518 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
519 2.76439e-02, 6.43207e-03, -3.54300e-02, 0.00000e+00, 0.00000e+00,
520 0.00000e+00, -2.80221e-02, 8.11228e+01, -6.75255e-04, 0.00000e+00,
521 -1.05162e-02, -3.48292e-03, -6.97321e-03, 0.00000e+00, 0.00000e+00,
522 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
523 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
524 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
525 0.00000e+00, -1.45546e-03, -1.31970e-02, -3.57751e-03, -1.09021e+00,
526 -1.50181e-02, -7.12841e-03, -6.64590e-03, -3.52610e-03, -1.87773e-02,
527 -2.22432e-03, -3.93895e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
528 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
529 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
530 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
531 },
532
533 {
534 6.04050e-02, 1.57034e+00, 2.99387e-02, 0.00000e+00, 0.00000e+00,
535 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -1.51018e+00,
536 0.00000e+00, 0.00000e+00, 0.00000e+00, -8.61650e+00, 1.26454e-02,
537 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
538 0.00000e+00, 5.50878e-03, 0.00000e+00, 0.00000e+00, 8.66784e-02,
539 1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
540 0.00000e+00, 0.00000e+00, 6.23881e-02, 0.00000e+00, 0.00000e+00,
541 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
542 0.00000e+00, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,
543 -9.45934e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
544 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
545 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
546 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
547 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
548 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
549 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
550 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
551 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
552 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
553 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
554 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
555 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
556 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
557 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
558 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
559 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
560 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
561 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
562 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
563 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
564 }
565 };
566
567
568 private static final double[] PS = {
569 9.56827e-01, 6.20637e-02, 3.18433e-02, 0.00000e+00, 0.00000e+00,
570 3.94900e-02, 0.00000e+00, 0.00000e+00, -9.24882e-03, -7.94023e-03,
571 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.74712e+02, 0.00000e+00,
572 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
573 0.00000e+00, 2.74677e-03, 0.00000e+00, 1.54951e-02, 8.66784e-02,
574 1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
575 0.00000e+00, 0.00000e+00, 0.00000e+00, -6.99007e-04, 0.00000e+00,
576 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
577 0.00000e+00, 1.24362e-02, -5.28756e-03, 8.47001e-02, 1.70147e-01,
578 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
579 0.00000e+00, 2.47425e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
580 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
581 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
582 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
583 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
584 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
585 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
586 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
587 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
588 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
589 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
590 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
591 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
592 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
593 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
594 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
595 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
596 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
597 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
598 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
599 };
600
601
602 private static final double[][] PDL = {
603 {
604 1.09930e+00, 3.90631e+00, 3.07165e+00, 9.86161e-01, 1.63536e+01,
605 4.63830e+00, 1.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
606 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
607 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
608 0.00000e+00, 0.00000e+00, 1.28840e+00, 3.10302e-02, 1.18339e-01
609 },
610 {
611 1.00000e+00, 7.00000e-01, 1.15020e+00, 3.44689e+00, 1.28840e+00,
612 1.00000e+00, 1.08738e+00, 1.22947e+00, 1.10016e+00, 7.34129e-01,
613 1.15241e+00, 2.22784e+00, 7.95046e-01, 4.01612e+00, 4.47749e+00,
614 1.23435e+02, -7.60535e-02, 1.68986e-06, 7.44294e-01, 1.03604e+00,
615 1.72783e+02, 1.15020e+00, 3.44689e+00, -7.46230e-01, 9.49154e-01
616 }
617 };
618
619
620 private static final double[] PTM = {
621 1.04130e+03, 3.86000e+02, 1.95000e+02, 1.66728e+01, 2.13000e+02,
622 1.20000e+02, 2.40000e+02, 1.87000e+02, -2.00000e+00, 0.00000e+00
623 };
624
625
626 private static final double[][] PDM = {
627 {
628 2.45600e+07, 6.71072e-06, 1.00000e+02, 0.00000e+00, 1.10000e+02,
629 1.00000e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
630 },
631 {
632 8.59400E+10, 1.00000e+00, 1.05000e+02, -8.00000e+00, 1.10000e+02,
633 1.00000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 0.00000e+00
634 },
635 {
636 2.81000E+11, 0.00000e+00, 1.05000e+02, 2.80000e+01, 2.89500e+01,
637 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
638 },
639 {
640 3.30000E+10, 2.68270e-01, 1.05000e+02, 1.00000e+00, 1.10000e+02,
641 1.00000e+01, 1.10000e+02, -1.00000e+01, 0.00000e+00, 0.00000e+00
642 },
643 {
644 1.33000e+09, 1.19615e-02, 1.05000e+02, 0.00000e+00, 1.10000e+02,
645 1.00000e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
646 },
647 {
648 1.76100e+05, 1.00000e+00, 9.50000e+01, -8.00000e+00, 1.10000e+02,
649 1.00000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 0.00000e+00,
650 },
651 {
652 1.00000e+07, 1.00000e+00, 1.05000e+02, -8.00000e+00, 1.10000e+02,
653 1.00000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 0.00000e+00
654 },
655 {
656 1.00000e+06, 1.00000e+00, 1.05000e+02, -8.00000e+00, 5.50000e+02,
657 7.60000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 4.00000e+03
658 }
659 };
660
661
662 private static final double[][] PTL = {
663
664 {
665 1.00858e+00, 4.56011e-02, -2.22972e-02, -5.44388e-02, 5.23136e-04,
666 -1.88849e-02, 5.23707e-02, -9.43646e-03, 6.31707e-03, -7.80460e-02,
667 -4.88430e-02, 0.00000e+00, 0.00000e+00, -7.60250e+00, 0.00000e+00,
668 -1.44635e-02, -1.76843e-02, -1.21517e+02, 2.85647e-02, 0.00000e+00,
669 0.00000e+00, 6.31792e-04, 0.00000e+00, 5.77197e-03, 8.66784e-02,
670 1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
671 0.00000e+00, -8.90272e+03, 3.30611e-03, 3.02172e-03, 0.00000e+00,
672 -2.13673e-03, -3.20910e-04, 0.00000e+00, 0.00000e+00, 2.76034e-03,
673 2.82487e-03, -2.97592e-04, -4.21534e-03, 8.47001e-02, 1.70147e-01,
674 8.96456e-03, 0.00000e+00, -1.08596e-02, 0.00000e+00, 0.00000e+00,
675 5.57917e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
676 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
677 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
678 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
679 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
680 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
681 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
682 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
683 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
684 0.00000e+00, 9.65405e-03, 0.00000e+00, 0.00000e+00, 2.00000e+00
685 },
686
687 {
688 9.39664e-01, 8.56514e-02, -6.79989e-03, 2.65929e-02, -4.74283e-03,
689 1.21855e-02, -2.14905e-02, 6.49651e-03, -2.05477e-02, -4.24952e-02,
690 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.19148e+01, 0.00000e+00,
691 1.18777e-02, -7.28230e-02, -8.15965e+01, 1.73887e-02, 0.00000e+00,
692 0.00000e+00, 0.00000e+00, -1.44691e-02, 2.80259e-04, 8.66784e-02,
693 1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
694 0.00000e+00, 2.16584e+02, 3.18713e-03, 7.37479e-03, 0.00000e+00,
695 -2.55018e-03, -3.92806e-03, 0.00000e+00, 0.00000e+00, -2.89757e-03,
696 -1.33549e-03, 1.02661e-03, 3.53775e-04, 8.47001e-02, 1.70147e-01,
697 -9.17497e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
698 3.56082e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
699 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
700 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
701 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
702 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
703 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
704 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
705 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
706 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
707 0.00000e+00, -1.00902e-02, 0.00000e+00, 0.00000e+00, 2.00000e+00
708 },
709
710 {
711 9.85982e-01, -4.55435e-02, 1.21106e-02, 2.04127e-02, -2.40836e-03,
712 1.11383e-02, -4.51926e-02, 1.35074e-02, -6.54139e-03, 1.15275e-01,
713 1.28247e-01, 0.00000e+00, 0.00000e+00, -5.30705e+00, 0.00000e+00,
714 -3.79332e-02, -6.24741e-02, 7.71062e-01, 2.96315e-02, 0.00000e+00,
715 0.00000e+00, 0.00000e+00, 6.81051e-03, -4.34767e-03, 8.66784e-02,
716 1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
717 0.00000e+00, 1.07003e+01, -2.76907e-03, 4.32474e-04, 0.00000e+00,
718 1.31497e-03, -6.47517e-04, 0.00000e+00, -2.20621e+01, -1.10804e-03,
719 -8.09338e-04, 4.18184e-04, 4.29650e-03, 8.47001e-02, 1.70147e-01,
720 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
721 -4.04337e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
722 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
723 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -9.52550e-04,
724 8.56253e-04, 4.33114e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
725 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.21223e-03,
726 2.38694e-04, 9.15245e-04, 1.28385e-03, 8.67668e-04, -5.61425e-06,
727 1.04445e+00, 3.41112e+01, 0.00000e+00, -8.40704e-01, -2.39639e+02,
728 7.06668e-01, -2.05873e+01, -3.63696e-01, 2.39245e+01, 0.00000e+00,
729 -1.06657e-03, -7.67292e-04, 1.54534e-04, 0.00000e+00, 0.00000e+00,
730 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
731 },
732
733 {
734 1.00320e+00, 3.83501e-02, -2.38983e-03, 2.83950e-03, 4.20956e-03,
735 5.86619e-04, 2.19054e-02, -1.00946e-02, -3.50259e-03, 4.17392e-02,
736 -8.44404e-03, 0.00000e+00, 0.00000e+00, 4.96949e+00, 0.00000e+00,
737 -7.06478e-03, -1.46494e-02, 3.13258e+01, -1.86493e-03, 0.00000e+00,
738 -1.67499e-02, 0.00000e+00, 0.00000e+00, 5.12686e-04, 8.66784e-02,
739 1.58727e-01, -4.64167e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
740 4.37353e-03, -1.99069e+02, 0.00000e+00, -5.34884e-03, 0.00000e+00,
741 1.62458e-03, 2.93016e-03, 2.67926e-03, 5.90449e+02, 0.00000e+00,
742 0.00000e+00, -1.17266e-03, -3.58890e-04, 8.47001e-02, 1.70147e-01,
743 0.00000e+00, 0.00000e+00, 1.38673e-02, 0.00000e+00, 0.00000e+00,
744 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
745 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
746 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.60571e-03,
747 6.28078e-04, 5.05469e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
748 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -1.57829e-03,
749 -4.00855e-04, 5.04077e-05, -1.39001e-03, -2.33406e-03, -4.81197e-04,
750 1.46758e+00, 6.20332e+00, 0.00000e+00, 3.66476e-01, -6.19760e+01,
751 3.09198e-01, -1.98999e+01, 0.00000e+00, -3.29933e+02, 0.00000e+00,
752 -1.10080e-03, -9.39310e-05, 1.39638e-04, 0.00000e+00, 0.00000e+00,
753 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
754 }
755 };
756
757
758 private static final double[][] PMA = {
759
760 {
761 9.81637e-01, -1.41317e-03, 3.87323e-02, 0.00000e+00, 0.00000e+00,
762 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -3.58707e-02,
763 -8.63658e-03, 0.00000e+00, 0.00000e+00, -2.02226e+00, 0.00000e+00,
764 -8.69424e-03, -1.91397e-02, 8.76779e+01, 4.52188e-03, 0.00000e+00,
765 2.23760e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
766 0.00000e+00, -7.07572e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
767 -4.11210e-03, 3.50060e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
768 0.00000e+00, 0.00000e+00, -8.36657e-03, 1.61347e+01, 0.00000e+00,
769 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
770 0.00000e+00, 0.00000e+00, -1.45130e-02, 0.00000e+00, 0.00000e+00,
771 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
772 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
773 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.24152e-03,
774 6.43365e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
775 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.33255e-03,
776 2.42657e-03, 1.60666e-03, -1.85728e-03, -1.46874e-03, -4.79163e-06,
777 1.22464e+00, 3.53510e+01, 0.00000e+00, 4.49223e-01, -4.77466e+01,
778 4.70681e-01, 8.41861e+00, -2.88198e-01, 1.67854e+02, 0.00000e+00,
779 7.11493e-04, 6.05601e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
780 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
781 },
782
783 {
784 1.00422e+00, -7.11212e-03, 5.24480e-03, 0.00000e+00, 0.00000e+00,
785 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -5.28914e-02,
786 -2.41301e-02, 0.00000e+00, 0.00000e+00, -2.12219e+01, -1.03830e-02,
787 -3.28077e-03, 1.65727e-02, 1.68564e+00, -6.68154e-03, 0.00000e+00,
788 1.45155e-02, 0.00000e+00, 8.42365e-03, 0.00000e+00, 0.00000e+00,
789 0.00000e+00, -4.34645e-03, 0.00000e+00, 0.00000e+00, 2.16780e-02,
790 0.00000e+00, -1.38459e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
791 0.00000e+00, 0.00000e+00, 7.04573e-03, -4.73204e+01, 0.00000e+00,
792 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
793 0.00000e+00, 0.00000e+00, 1.08767e-02, 0.00000e+00, 0.00000e+00,
794 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
795 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -8.08279e-03,
796 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.21769e-04,
797 -2.27387e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
798 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.26769e-03,
799 3.16901e-03, 4.60316e-04, -1.01431e-04, 1.02131e-03, 9.96601e-04,
800 1.25707e+00, 2.50114e+01, 0.00000e+00, 4.24472e-01, -2.77655e+01,
801 3.44625e-01, 2.75412e+01, 0.00000e+00, 7.94251e+02, 0.00000e+00,
802 2.45835e-03, 1.38871e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
803 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
804 },
805
806 {
807 1.01890e+00, -2.46603e-02, 1.00078e-02, 0.00000e+00, 0.00000e+00,
808 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -6.70977e-02,
809 -4.02286e-02, 0.00000e+00, 0.00000e+00, -2.29466e+01, -7.47019e-03,
810 2.26580e-03, 2.63931e-02, 3.72625e+01, -6.39041e-03, 0.00000e+00,
811 9.58383e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
812 0.00000e+00, -1.85291e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
813 0.00000e+00, 1.39717e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
814 0.00000e+00, 0.00000e+00, 9.19771e-03, -3.69121e+02, 0.00000e+00,
815 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
816 0.00000e+00, 0.00000e+00, -1.57067e-02, 0.00000e+00, 0.00000e+00,
817 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
818 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -7.07265e-03,
819 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.92953e-03,
820 -2.77739e-03, -4.40092e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
821 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.47280e-03,
822 2.95035e-04, -1.81246e-03, 2.81945e-03, 4.27296e-03, 9.78863e-04,
823 1.40545e+00, -6.19173e+00, 0.00000e+00, 0.00000e+00, -7.93632e+01,
824 4.44643e-01, -4.03085e+02, 0.00000e+00, 1.15603e+01, 0.00000e+00,
825 2.25068e-03, 8.48557e-04, -2.98493e-04, 0.00000e+00, 0.00000e+00,
826 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
827 },
828
829 {
830 9.75801e-01, 3.80680e-02, -3.05198e-02, 0.00000e+00, 0.00000e+00,
831 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.85575e-02,
832 5.04057e-02, 0.00000e+00, 0.00000e+00, -1.76046e+02, 1.44594e-02,
833 -1.48297e-03, -3.68560e-03, 3.02185e+01, -3.23338e-03, 0.00000e+00,
834 1.53569e-02, 0.00000e+00, -1.15558e-02, 0.00000e+00, 0.00000e+00,
835 0.00000e+00, 4.89620e-03, 0.00000e+00, 0.00000e+00, -1.00616e-02,
836 -8.21324e-03, -1.57757e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
837 0.00000e+00, 0.00000e+00, 6.63564e-03, 4.58410e+01, 0.00000e+00,
838 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
839 0.00000e+00, 0.00000e+00, -2.51280e-02, 0.00000e+00, 0.00000e+00,
840 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
841 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 9.91215e-03,
842 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -8.73148e-04,
843 -1.29648e-03, -7.32026e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
844 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -4.68110e-03,
845 -4.66003e-03, -1.31567e-03, -7.39390e-04, 6.32499e-04, -4.65588e-04,
846 -1.29785e+00, -1.57139e+02, 0.00000e+00, 2.58350e-01, -3.69453e+01,
847 4.10672e-01, 9.78196e+00, -1.52064e-01, -3.85084e+03, 0.00000e+00,
848 -8.52706e-04, -1.40945e-03, -7.26786e-04, 0.00000e+00, 0.00000e+00,
849 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
850 },
851
852 {
853 9.60722e-01, 7.03757e-02, -3.00266e-02, 0.00000e+00, 0.00000e+00,
854 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.22671e-02,
855 4.10423e-02, 0.00000e+00, 0.00000e+00, -1.63070e+02, 1.06073e-02,
856 5.40747e-04, 7.79481e-03, 1.44908e+02, 1.51484e-04, 0.00000e+00,
857 1.97547e-02, 0.00000e+00, -1.41844e-02, 0.00000e+00, 0.00000e+00,
858 0.00000e+00, 5.77884e-03, 0.00000e+00, 0.00000e+00, 9.74319e-03,
859 0.00000e+00, -2.88015e+03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
860 0.00000e+00, 0.00000e+00, -4.44902e-03, -2.92760e+01, 0.00000e+00,
861 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
862 0.00000e+00, 0.00000e+00, 2.34419e-02, 0.00000e+00, 0.00000e+00,
863 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
864 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.36685e-03,
865 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -4.65325e-04,
866 -5.50628e-04, 3.31465e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
867 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.06179e-03,
868 -3.08575e-03, -7.93589e-04, -1.08629e-04, 5.95511e-04, -9.05050e-04,
869 1.18997e+00, 4.15924e+01, 0.00000e+00, -4.72064e-01, -9.47150e+02,
870 3.98723e-01, 1.98304e+01, 0.00000e+00, 3.73219e+03, 0.00000e+00,
871 -1.50040e-03, -1.14933e-03, -1.56769e-04, 0.00000e+00, 0.00000e+00,
872 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
873 },
874
875 {
876 1.03123e+00, -7.05124e-02, 8.71615e-03, 0.00000e+00, 0.00000e+00,
877 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -3.82621e-02,
878 -9.80975e-03, 0.00000e+00, 0.00000e+00, 2.89286e+01, 9.57341e-03,
879 0.00000e+00, 0.00000e+00, 8.66153e+01, 7.91938e-04, 0.00000e+00,
880 0.00000e+00, 0.00000e+00, 4.68917e-03, 0.00000e+00, 0.00000e+00,
881 0.00000e+00, 7.86638e-03, 0.00000e+00, 0.00000e+00, 9.90827e-03,
882 0.00000e+00, 6.55573e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
883 0.00000e+00, 0.00000e+00, 0.00000e+00, -4.00200e+01, 0.00000e+00,
884 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
885 0.00000e+00, 0.00000e+00, 7.07457e-03, 0.00000e+00, 0.00000e+00,
886 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
887 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.72268e-03,
888 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.04970e-04,
889 1.21560e-03, -8.05579e-06, 0.00000e+00, 0.00000e+00, 0.00000e+00,
890 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.49941e-03,
891 -4.57256e-04, -1.59311e-04, 2.96481e-04, -1.77318e-03, -6.37918e-04,
892 1.02395e+00, 1.28172e+01, 0.00000e+00, 1.49903e-01, -2.63818e+01,
893 0.00000e+00, 4.70628e+01, -2.22139e-01, 4.82292e-02, 0.00000e+00,
894 -8.67075e-04, -5.86479e-04, 5.32462e-04, 0.00000e+00, 0.00000e+00,
895 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
896 },
897
898 {
899 1.00828e+00, -9.10404e-02, -2.26549e-02, 0.00000e+00, 0.00000e+00,
900 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.32420e-02,
901 -9.08925e-03, 0.00000e+00, 0.00000e+00, 3.36105e+01, 0.00000e+00,
902 0.00000e+00, 0.00000e+00, -1.24957e+01, -5.87939e-03, 0.00000e+00,
903 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
904 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
905 0.00000e+00, 2.79765e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
906 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.01237e+03, 0.00000e+00,
907 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
908 0.00000e+00, 0.00000e+00, -1.75553e-02, 0.00000e+00, 0.00000e+00,
909 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
910 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
911 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.29699e-03,
912 1.26659e-03, 2.68402e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
913 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.17894e-03,
914 1.48746e-03, 1.06478e-04, 1.34743e-04, -2.20939e-03, -6.23523e-04,
915 6.36539e-01, 1.13621e+01, 0.00000e+00, -3.93777e-01, 2.38687e+03,
916 0.00000e+00, 6.61865e+02, -1.21434e-01, 9.27608e+00, 0.00000e+00,
917 1.68478e-04, 1.24892e-03, 1.71345e-03, 0.00000e+00, 0.00000e+00,
918 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
919 },
920
921 {
922 1.57293e+00, -6.78400e-01, 6.47500e-01, 0.00000e+00, 0.00000e+00,
923 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -7.62974e-02,
924 -3.60423e-01, 0.00000e+00, 0.00000e+00, 1.28358e+02, 0.00000e+00,
925 0.00000e+00, 0.00000e+00, 4.68038e+01, 0.00000e+00, 0.00000e+00,
926 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
927 0.00000e+00, -1.67898e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
928 0.00000e+00, 2.90994e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
929 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.15706e+01, 0.00000e+00,
930 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
931 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
932 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
933 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
934 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
935 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
936 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
937 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
938 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
939 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
940 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
941 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
942 },
943
944 {
945 8.60028e-01, 3.77052e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
946 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -1.17570e+00,
947 0.00000e+00, 0.00000e+00, 0.00000e+00, 7.77757e-03, 0.00000e+00,
948 0.00000e+00, 0.00000e+00, 1.01024e+02, 0.00000e+00, 0.00000e+00,
949 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
950 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
951 0.00000e+00, 6.54251e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
952 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
953 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
954 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
955 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
956 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
957 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
958 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
959 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -1.56959e-02,
960 1.91001e-02, 3.15971e-02, 1.00982e-02, -6.71565e-03, 2.57693e-03,
961 1.38692e+00, 2.82132e-01, 0.00000e+00, 0.00000e+00, 3.81511e+02,
962 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
963 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
964 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
965 },
966
967 {
968 1.06029e+00, -5.25231e-02, 3.73034e-01, 0.00000e+00, 0.00000e+00,
969 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.31072e-02,
970 -3.88409e-01, 0.00000e+00, 0.00000e+00, -1.65295e+02, -2.13801e-01,
971 -4.38916e-02, -3.22716e-01, -8.82393e+01, 1.18458e-01, 0.00000e+00,
972 -4.35863e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
973 0.00000e+00, -1.19782e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
974 0.00000e+00, 2.62229e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
975 0.00000e+00, 0.00000e+00, 0.00000e+00, -5.37443e+01, 0.00000e+00,
976 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
977 0.00000e+00, 0.00000e+00, -4.55788e-01, 0.00000e+00, 0.00000e+00,
978 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
979 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
980 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.84009e-02,
981 3.96733e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
982 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.05494e-02,
983 7.39617e-02, 1.92200e-02, -8.46151e-03, -1.34244e-02, 1.96338e-02,
984 1.50421e+00, 1.88368e+01, 0.00000e+00, 0.00000e+00, -5.13114e+01,
985 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
986 5.11923e-02, 3.61225e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
987 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
988 }
989 };
990
991
992 private static final double[] PAVGM = {
993 2.61000e+02, 2.64000e+02, 2.29000e+02, 2.17000e+02, 2.17000e+02,
994 2.23000e+02, 2.86760e+02, -2.93940e+00, 2.50000e+00, 0.00000e+00
995 };
996
997
998 private static final double MIN_TEMP = 50.;
999
1000
1001
1002
1003 private final NRLMSISE00InputParameters inputParams;
1004
1005
1006 private final BodyShape earth;
1007
1008
1009 private final int[] sw;
1010
1011
1012 private final int[] swc;
1013
1014
1015 private final TimeScale ut;
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035 @DefaultDataContext
1036 public NRLMSISE00(final NRLMSISE00InputParameters parameters,
1037 final ExtendedPositionProvider sun,
1038 final BodyShape earth) {
1039 this(parameters, sun, earth,
1040 DataContext.getDefault().getTimeScales()
1041 .getUT1(IERSConventions.IERS_2010, true));
1042 }
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061 public NRLMSISE00(final NRLMSISE00InputParameters parameters,
1062 final ExtendedPositionProvider sun,
1063 final BodyShape earth,
1064 final TimeScale ut) {
1065 this(parameters, sun, earth, allOnes(), allOnes(), ut);
1066 }
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084 private NRLMSISE00(final NRLMSISE00InputParameters parameters,
1085 final ExtendedPositionProvider sun,
1086 final BodyShape earth,
1087 final int[] sw,
1088 final int[] swc,
1089 final TimeScale ut) {
1090 super(sun);
1091 this.inputParams = parameters;
1092 this.earth = earth;
1093 this.sw = sw;
1094 this.swc = swc;
1095 this.ut = ut;
1096 }
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107 public NRLMSISE00 withSwitch(final int number, final int value) {
1108 if (number < 1 || number > 23) {
1109 throw new OrekitException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, number, 1, 23);
1110 }
1111
1112 final int[] newSw = sw.clone();
1113 final int[] newSwc = swc.clone();
1114 if (number != 9) {
1115 newSw[number] = (value == 1) ? 1 : 0;
1116 newSwc[number] = (value > 0) ? 1 : 0;
1117 } else {
1118 if (value == -1 || value == 1) {
1119 newSw[number] = value;
1120 } else {
1121 newSw[number] = 0;
1122 }
1123 newSwc[number] = newSw[number];
1124 }
1125
1126 return new NRLMSISE00(inputParams, getSun(), earth, newSwc, newSwc, ut);
1127
1128 }
1129
1130
1131
1132
1133 private static int[] allOnes() {
1134 final int[] array = new int[24];
1135 Arrays.fill(array, 1);
1136 return array;
1137 }
1138
1139
1140 @Override
1141 public Frame getFrame() {
1142 return earth.getBodyFrame();
1143 }
1144
1145
1146 @Override
1147 public double getDensity(final AbsoluteDate date,
1148 final Vector3D position,
1149 final Frame frame) {
1150
1151
1152 if (!date.isBetweenOrEqualTo(inputParams.getMinDate(), inputParams.getMaxDate())) {
1153 throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
1154 date, inputParams.getMinDate(), inputParams.getMaxDate());
1155 }
1156
1157
1158 final DateTimeComponents dtc = date.getComponents(ut);
1159 final int doy = dtc.getDate().getDayOfYear();
1160 final double sec = dtc.getTime().getSecondsInLocalDay();
1161
1162
1163 final GeodeticPoint inBody = earth.transform(position, frame, date);
1164 final double alt = inBody.getAltitude() / 1000.;
1165 final double lon = FastMath.toDegrees(inBody.getLongitude());
1166 final double lat = FastMath.toDegrees(inBody.getLatitude());
1167
1168
1169 final double lst = localSolarTime(date, position, frame);
1170
1171
1172 final Output out = new Output(doy, sec, lat, lon, lst, inputParams.getAverageFlux(date),
1173 inputParams.getDailyFlux(date), inputParams.getAp(date));
1174 out.gtd7d(alt);
1175
1176
1177 return out.getDensity(TOTAL_MASS);
1178
1179 }
1180
1181
1182 @Override
1183 public <T extends CalculusFieldElement<T>> T getDensity(final FieldAbsoluteDate<T> date,
1184 final FieldVector3D<T> position,
1185 final Frame frame) {
1186
1187 final AbsoluteDate dateD = date.toAbsoluteDate();
1188 if (!dateD.isBetweenOrEqualTo(inputParams.getMinDate(), inputParams.getMaxDate())) {
1189 throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
1190 dateD, inputParams.getMinDate(), inputParams.getMaxDate());
1191 }
1192
1193
1194 final DateTimeComponents dtc = dateD.getComponents(ut);
1195 final int doy = dtc.getDate().getDayOfYear();
1196 final T sec = date.durationFrom(new AbsoluteDate(dtc.getDate(), TimeComponents.H00, ut));
1197
1198
1199 final FieldGeodeticPoint<T> inBody = earth.transform(position, frame, date);
1200 final T alt = inBody.getAltitude().divide(1000.);
1201 final T lon = FastMath.toDegrees(inBody.getLongitude());
1202 final T lat = FastMath.toDegrees(inBody.getLatitude());
1203
1204
1205 final T lst = localSolarTime(date, position, frame);
1206
1207
1208 final FieldOutput<T> out = new FieldOutput<>(doy, sec, lat, lon, lst,
1209 inputParams.getAverageFlux(dateD),
1210 inputParams.getDailyFlux(dateD), inputParams.getAp(dateD));
1211 out.gtd7d(alt);
1212
1213
1214 return out.getDensity(TOTAL_MASS);
1215
1216 }
1217
1218
1219
1220
1221
1222
1223
1224 private double localSolarTime(final AbsoluteDate date,
1225 final Vector3D position,
1226 final Frame frame) {
1227 final Vector3D sunPos = getSunPosition(date, frame);
1228 final double lst = FastMath.PI + FastMath.atan2(
1229 sunPos.getX() * position.getY() - sunPos.getY() * position.getX(),
1230 sunPos.getX() * position.getX() + sunPos.getY() * position.getY());
1231 return lst * 12. / FastMath.PI;
1232 }
1233
1234
1235
1236
1237
1238
1239
1240
1241 private <T extends CalculusFieldElement<T>> T localSolarTime(final FieldAbsoluteDate<T> date,
1242 final FieldVector3D<T> position,
1243 final Frame frame) {
1244 final FieldVector3D<T> sunPos = getSunPosition(date, frame);
1245 final T y = position.getY().multiply(sunPos.getX()).subtract(position.getX().multiply(sunPos.getY()));
1246 final T x = position.getX().multiply(sunPos.getX()).add(position.getY().multiply(sunPos.getY()));
1247 final T hl = y.atan2(x).add(y.getPi());
1248
1249 return hl.divide(y.getPi()).multiply(12.);
1250
1251 }
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287 private class Output {
1288
1289
1290 private final int doy;
1291
1292
1293 private final double sec;
1294
1295
1296 private final double lat;
1297
1298
1299 private final double lon;
1300
1301
1302 private final double hl;
1303
1304
1305 private final double f107a;
1306
1307
1308 private final double f107;
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320 private final double[] ap;
1321
1322
1323 private final double glat;
1324
1325
1326 private final double rlat;
1327
1328
1329 private double dm28;
1330
1331
1332 private final double[][] plg;
1333
1334
1335 private final double ctloc;
1336
1337 private final double stloc;
1338
1339 private final double c2tloc;
1340
1341 private final double s2tloc;
1342
1343 private final double c3tloc;
1344
1345 private final double s3tloc;
1346
1347
1348 private double apdf;
1349
1350
1351 private double apt;
1352
1353
1354 private final double[] meso_tn1;
1355
1356
1357 private final double[] meso_tn2;
1358
1359
1360 private final double[] meso_tn3;
1361
1362
1363 private final double[] meso_tgn1;
1364
1365
1366 private final double[] meso_tgn2;
1367
1368
1369 private final double[] meso_tgn3;
1370
1371
1372 private final double[] densities;
1373
1374
1375 private final double[] temperatures;
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396 Output(final int doy, final double sec,
1397 final double lat, final double lon, final double hl,
1398 final double f107a, final double f107, final double[] ap) {
1399
1400 this.doy = doy;
1401 this.sec = sec;
1402 this.lat = lat;
1403 this.lon = lon;
1404 this.hl = hl;
1405 this.f107a = f107a;
1406 this.f107 = f107;
1407 this.ap = ap.clone();
1408
1409 this.plg = new double[4][8];
1410
1411 this.meso_tn1 = new double[ZN1.length];
1412 this.meso_tn2 = new double[ZN2.length];
1413 this.meso_tn3 = new double[ZN3.length];
1414 this.meso_tgn1 = new double[2];
1415 this.meso_tgn2 = new double[2];
1416 this.meso_tgn3 = new double[2];
1417
1418 densities = new double[9];
1419 temperatures = new double[2];
1420
1421
1422 final double xlat = (sw[2] == 0) ? LAT_REF : lat;
1423 final double c2 = FastMath.cos(2 * DEG_TO_RAD * xlat);
1424 glat = G_REF * (1. - .0026373 * c2);
1425 rlat = 2. * glat / (3.085462e-6 + 2.27e-9 * c2) * 1.e-5;
1426
1427
1428 final double latr = DEG_TO_RAD * lat;
1429
1430
1431 final SinCos scLatr = FastMath.sinCos(latr);
1432 final double c = scLatr.sin();
1433 final double s = scLatr.cos();
1434
1435 plg[0][1] = c;
1436 plg[0][2] = ( 3.0 * c * plg[0][1] - 1.0) / 2.0;
1437 plg[0][3] = ( 5.0 * c * plg[0][2] - 2.0 * plg[0][1]) / 3.0;
1438 plg[0][4] = ( 7.0 * c * plg[0][3] - 3.0 * plg[0][2]) / 4.0;
1439 plg[0][5] = ( 9.0 * c * plg[0][4] - 4.0 * plg[0][3]) / 5.0;
1440 plg[0][6] = (11.0 * c * plg[0][5] - 5.0 * plg[0][4]) / 6.0;
1441
1442 plg[1][1] = s;
1443 plg[1][2] = 3.0 * c * plg[1][1];
1444 plg[1][3] = ( 5.0 * c * plg[1][2] - 3.0 * plg[1][1]) / 2.0;
1445 plg[1][4] = ( 7.0 * c * plg[1][3] - 4.0 * plg[1][2]) / 3.0;
1446 plg[1][5] = ( 9.0 * c * plg[1][4] - 5.0 * plg[1][3]) / 4.0;
1447 plg[1][6] = (11.0 * c * plg[1][5] - 6.0 * plg[1][4]) / 5.0;
1448
1449 plg[2][2] = 3.0 * s * plg[1][1];
1450 plg[2][3] = 5.0 * c * plg[2][2];
1451 plg[2][4] = ( 7.0 * c * plg[2][3] - 5.0 * plg[2][2]) / 2.0;
1452 plg[2][5] = ( 9.0 * c * plg[2][4] - 6.0 * plg[2][3]) / 3.0;
1453 plg[2][6] = (11.0 * c * plg[2][5] - 7.0 * plg[2][4]) / 4.0;
1454 plg[2][7] = (13.0 * c * plg[2][6] - 8.0 * plg[2][5]) / 5.0;
1455
1456 plg[3][3] = 5.0 * s * plg[2][2];
1457 plg[3][4] = 7.0 * c * plg[3][3];
1458 plg[3][5] = ( 9.0 * c * plg[3][4] - 7.0 * plg[3][3]) / 2.0;
1459 plg[3][6] = (11.0 * c * plg[3][5] - 8.0 * plg[3][4]) / 3.0;
1460
1461
1462 if (!(sw[7] == 0 && sw[8] == 0 && sw[14] == 0)) {
1463 final double tloc = HOUR_TO_RAD * hl;
1464 final SinCos sc = FastMath.sinCos(tloc);
1465 final SinCos sc2 = SinCos.sum(sc, sc);
1466 final SinCos sc3 = SinCos.sum(sc, sc2);
1467 stloc = sc.sin();
1468 ctloc = sc.cos();
1469 s2tloc = sc2.sin();
1470 c2tloc = sc2.cos();
1471 s3tloc = sc3.sin();
1472 c3tloc = sc3.cos();
1473 } else {
1474 stloc = 0;
1475 ctloc = 0;
1476 s2tloc = 0;
1477 c2tloc = 0;
1478 s3tloc = 0;
1479 c3tloc = 0;
1480 }
1481
1482 }
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508 void gts7(final double alt) {
1509
1510
1511 final double[] alpha = {-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0};
1512
1513 final double[] altl = {200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0};
1514
1515 final double xmm = PDM[2][4];
1516
1517
1518 double tinf = PTM[0] * PT[0];
1519
1520 if (alt > ZN1[0]) {
1521 tinf *= 1.0 + sw[16] * globe7(PT);
1522 }
1523 setTemperature(EXOSPHERIC, tinf);
1524
1525
1526 double g0 = PTM[3] * PS[0];
1527 if (alt > ZN1[4]) {
1528 g0 *= 1.0 + sw[19] * globe7(PS);
1529 }
1530
1531
1532 double tlb = PTM[1] * PD[3][0];
1533 tlb *= 1.0 + sw[17] * globe7(PD[3]);
1534
1535
1536 final double s = g0 / (tinf - tlb);
1537
1538
1539 meso_tn1[1] = PTM[6] * PTL[0][0];
1540 meso_tn1[2] = PTM[2] * PTL[1][0];
1541 meso_tn1[3] = PTM[7] * PTL[2][0];
1542 meso_tn1[4] = PTM[4] * PTL[3][0];
1543 meso_tgn1[1] = PTM[8] * PMA[8][0];
1544 if (alt < 300.0) {
1545 final double r = PTM[4] * PTL[3][0];
1546 meso_tn1[1] /= 1.0 - sw[18] * glob7s(PTL[0]);
1547 meso_tn1[2] /= 1.0 - sw[18] * glob7s(PTL[1]);
1548 meso_tn1[3] /= 1.0 - sw[18] * glob7s(PTL[2]);
1549 meso_tn1[4] /= 1.0 - sw[18] * sw[20] * glob7s(PTL[3]);
1550 meso_tgn1[1] *= 1.0 + sw[18] * sw[20] * glob7s(PMA[8]);
1551 meso_tgn1[1] *= meso_tn1[4] * meso_tn1[4] / (r * r);
1552 }
1553
1554
1555 setTemperature(ALTITUDE, densu(alt, 1.0, tinf, tlb, 0.0, 0.0, PTM[5], s));
1556
1557
1558
1559 final double g28 = sw[21] * globe7(PD[2]);
1560
1561 final double db28 = PDM[2][0] * FastMath.exp(g28) * PD[2][0];
1562
1563 double diffusiveDensity = densu(alt, db28, tinf, tlb, N2_MASS, alpha[2], PTM[5], s);
1564 setDensity(MOLECULAR_NITROGEN, diffusiveDensity);
1565
1566 final double zhf = PDL[1][24] * (1.0 + sw[5] * PDL[0][24] *
1567 FastMath.sin(DEG_TO_RAD * lat) *
1568 FastMath.cos(DAY_TO_RAD * (doy - PT[13])));
1569
1570 final double zh28 = PDM[2][2] * zhf;
1571 final double zhm28 = PDM[2][3] * PDL[1][5];
1572
1573 final double b28 = densu(zh28, db28, tinf, tlb, N2_MASS - xmm, alpha[2] - 1.0, PTM[5], s);
1574 if (sw[15] != 0 && alt <= altl[2]) {
1575
1576 dm28 = densu(alt, b28, tinf, tlb, xmm, alpha[2], PTM[5], s);
1577
1578 setDensity(MOLECULAR_NITROGEN, dnet(diffusiveDensity, dm28, zhm28, xmm, N2_MASS));
1579 }
1580
1581
1582
1583 final double g4 = sw[21] * globe7(PD[0]);
1584
1585 final double db04 = PDM[0][0] * FastMath.exp(g4) * PD[0][0];
1586
1587 diffusiveDensity = densu(alt, db04, tinf, tlb, HE_MASS, alpha[0], PTM[5], s);
1588 setDensity(HELIUM, diffusiveDensity);
1589 if (sw[15] != 0 && alt < altl[0]) {
1590
1591 final double zh04 = PDM[0][2];
1592
1593 final double b04 = densu(zh04, db04, tinf, tlb, HE_MASS - xmm, alpha[0] - 1., PTM[5], s);
1594
1595 final double dm04 = densu(alt, b04, tinf, tlb, xmm, 0., PTM[5], s);
1596 final double zhm04 = zhm28;
1597
1598 diffusiveDensity = dnet(diffusiveDensity, dm04, zhm04, xmm, HE_MASS);
1599
1600 final double rl = FastMath.log(b28 * PDM[0][1] / b04);
1601 final double zc04 = PDM[0][4] * PDL[1][0];
1602 final double hc04 = PDM[0][5] * PDL[1][1];
1603
1604 setDensity(HELIUM, diffusiveDensity * ccor(alt, rl, hc04, zc04));
1605 }
1606
1607
1608
1609 final double g16 = sw[21] * globe7(PD[1]);
1610
1611 final double db16 = PDM[1][0] * FastMath.exp(g16) * PD[1][0];
1612
1613 diffusiveDensity = densu(alt, db16, tinf, tlb, O_MASS, alpha[1], PTM[5], s);
1614 setDensity(ATOMIC_OXYGEN, diffusiveDensity);
1615 if (sw[15] != 0 && alt < altl[1]) {
1616
1617 final double zh16 = PDM[1][2];
1618
1619 final double b16 = densu(zh16, db16, tinf, tlb, O_MASS - xmm, alpha[1] - 1.0, PTM[5], s);
1620
1621 final double dm16 = densu(alt, b16, tinf, tlb, xmm, 0., PTM[5], s);
1622 final double zhm16 = zhm28;
1623
1624 diffusiveDensity = dnet(diffusiveDensity, dm16, zhm16, xmm, O_MASS);
1625 final double rl = PDM[1][1] * PDL[1][16] * (1.0 + sw[1] * PDL[0][23] * (f107a - FLUX_REF));
1626 final double hc16 = PDM[1][5] * PDL[1][3];
1627 final double zc16 = PDM[1][4] * PDL[1][2];
1628 final double hc216 = PDM[1][5] * PDL[1][4];
1629 diffusiveDensity *= ccor2(alt, rl, hc16, zc16, hc216);
1630
1631 final double hcc16 = PDM[1][7] * PDL[1][13];
1632 final double zcc16 = PDM[1][6] * PDL[1][12];
1633 final double rc16 = PDM[1][3] * PDL[1][14];
1634
1635 setDensity(ATOMIC_OXYGEN, diffusiveDensity * ccor(alt, rc16, hcc16, zcc16));
1636 }
1637
1638
1639
1640 final double g32 = sw[21] * globe7(PD[4]);
1641
1642 final double db32 = PDM[3][0] * FastMath.exp(g32) * PD[4][0];
1643
1644 diffusiveDensity = densu(alt, db32, tinf, tlb, O2_MASS, alpha[3], PTM[5], s);
1645 setDensity(MOLECULAR_OXYGEN, diffusiveDensity);
1646 if (sw[15] != 0) {
1647 if (alt <= altl[3]) {
1648
1649 final double zh32 = PDM[3][2];
1650
1651 final double b32 = densu(zh32, db32, tinf, tlb, O2_MASS - xmm, alpha[3] - 1., PTM[5], s);
1652
1653 final double dm32 = densu(alt, b32, tinf, tlb, xmm, 0., PTM[5], s);
1654 final double zhm32 = zhm28;
1655
1656 diffusiveDensity = dnet(diffusiveDensity, dm32, zhm32, xmm, O2_MASS);
1657
1658 final double rl = FastMath.log(b28 * PDM[3][1] / b32);
1659 final double hc32 = PDM[3][5] * PDL[1][7];
1660 final double zc32 = PDM[3][4] * PDL[1][6];
1661 diffusiveDensity *= ccor(alt, rl, hc32, zc32);
1662 }
1663
1664 final double hcc32 = PDM[3][7] * PDL[1][22];
1665 final double hcc232 = PDM[3][7] * PDL[0][22];
1666 final double zcc32 = PDM[3][6] * PDL[1][21];
1667 final double rc32 = PDM[3][3] * PDL[1][23] * (1. + sw[1] * PDL[0][23] * (f107a - FLUX_REF));
1668
1669 setDensity(MOLECULAR_OXYGEN, diffusiveDensity * ccor2(alt, rc32, hcc32, zcc32, hcc232));
1670 }
1671
1672
1673
1674 final double g40 = sw[21] * globe7(PD[5]);
1675
1676 final double db40 = PDM[4][0] * FastMath.exp(g40) * PD[5][0];
1677
1678 diffusiveDensity = densu(alt, db40, tinf, tlb, AR_MASS, alpha[4], PTM[5], s);
1679 setDensity(ARGON, diffusiveDensity);
1680 if (sw[15] != 0 && alt <= altl[4]) {
1681
1682 final double zh40 = PDM[4][2];
1683
1684 final double b40 = densu(zh40, db40, tinf, tlb, AR_MASS - xmm, alpha[4] - 1., PTM[5], s);
1685
1686 final double dm40 = densu(alt, b40, tinf, tlb, xmm, 0., PTM[5], s);
1687 final double zhm40 = zhm28;
1688
1689 diffusiveDensity = dnet(diffusiveDensity, dm40, zhm40, xmm, AR_MASS);
1690
1691 final double rl = FastMath.log(b28 * PDM[4][1] / b40);
1692 final double hc40 = PDM[4][5] * PDL[1][9];
1693 final double zc40 = PDM[4][4] * PDL[1][8];
1694
1695 setDensity(ARGON, diffusiveDensity * ccor(alt, rl, hc40, zc40));
1696 }
1697
1698
1699
1700 final double g1 = sw[21] * globe7(PD[6]);
1701
1702 final double db01 = PDM[5][0] * FastMath.exp(g1) * PD[6][0];
1703
1704 diffusiveDensity = densu(alt, db01, tinf, tlb, H_MASS, alpha[6], PTM[5], s);
1705 setDensity(HYDROGEN, diffusiveDensity);
1706 if (sw[15] != 0 && alt <= altl[6]) {
1707
1708 final double zh01 = PDM[5][2];
1709
1710 final double b01 = densu(zh01, db01, tinf, tlb, H_MASS - xmm, alpha[6] - 1., PTM[5], s);
1711
1712 final double dm01 = densu(alt, b01, tinf, tlb, xmm, 0., PTM[5], s);
1713 final double zhm01 = zhm28;
1714
1715 diffusiveDensity = dnet(diffusiveDensity, dm01, zhm01, xmm, H_MASS);
1716
1717 final double rl = FastMath.log(b28 * PDM[5][1] * FastMath.sqrt(PDL[1][17] * PDL[1][17]) / b01);
1718 final double hc01 = PDM[5][5] * PDL[1][11];
1719 final double zc01 = PDM[5][4] * PDL[1][10];
1720 diffusiveDensity *= ccor(alt, rl, hc01, zc01);
1721
1722 final double hcc01 = PDM[5][7] * PDL[1][19];
1723 final double zcc01 = PDM[5][6] * PDL[1][18];
1724 final double rc01 = PDM[5][3] * PDL[1][20];
1725
1726 setDensity(HYDROGEN, diffusiveDensity * ccor(alt, rc01, hcc01, zcc01));
1727 }
1728
1729
1730
1731 final double g14 = sw[21] * globe7(PD[7]);
1732
1733 final double db14 = PDM[6][0] * FastMath.exp(g14) * PD[7][0];
1734
1735 diffusiveDensity = densu(alt, db14, tinf, tlb, N_MASS, alpha[7], PTM[5], s);
1736 setDensity(ATOMIC_NITROGEN, diffusiveDensity);
1737 if (sw[15] != 0 && alt <= altl[7]) {
1738
1739 final double zh14 = PDM[6][2];
1740
1741 final double b14 = densu(zh14, db14, tinf, tlb, N_MASS - xmm, alpha[7] - 1., PTM[5], s);
1742
1743 final double dm14 = densu(alt, b14, tinf, tlb, xmm, 0., PTM[5], s);
1744 final double zhm14 = zhm28;
1745
1746 diffusiveDensity = dnet(diffusiveDensity, dm14, zhm14, xmm, N_MASS);
1747
1748 final double rl = FastMath.log(b28 * PDM[6][1] * PDL[0][2] / b14);
1749 final double hc14 = PDM[6][5] * PDL[0][1];
1750 final double zc14 = PDM[6][4] * PDL[0][0];
1751 diffusiveDensity *= ccor(alt, rl, hc14, zc14);
1752
1753 final double hcc14 = PDM[6][7] * PDL[0][4];
1754 final double zcc14 = PDM[6][6] * PDL[0][3];
1755 final double rc14 = PDM[6][3] * PDL[0][5];
1756
1757 setDensity(ATOMIC_NITROGEN, diffusiveDensity * ccor(alt, rc14, hcc14, zcc14));
1758 }
1759
1760
1761 final double g16h = sw[21] * globe7(PD[8]);
1762 final double db16h = PDM[7][0] * FastMath.exp(g16h) * PD[8][0];
1763 final double tho = PDM[7][9] * PDL[0][6];
1764 diffusiveDensity = densu(alt, db16h, tho, tho, O_MASS, alpha[8], PTM[5], s);
1765 final double zsht = PDM[7][5];
1766 final double zmho = PDM[7][4];
1767 final double zsho = scalh(zmho, O_MASS, tho);
1768 diffusiveDensity *= FastMath.exp(-zsht / zsho * (FastMath.exp((zmho - alt ) / zsht) - 1.));
1769 setDensity(ANOMALOUS_OXYGEN, diffusiveDensity);
1770
1771
1772 for (int i = 0; i < 9; i++) {
1773 setDensity(i, getDensity(i) * 1.0e+06);
1774 }
1775
1776
1777 final double tmd = AMU * (HE_MASS * getDensity(HELIUM) +
1778 O_MASS * getDensity(ATOMIC_OXYGEN) +
1779 N2_MASS * getDensity(MOLECULAR_NITROGEN) +
1780 O2_MASS * getDensity(MOLECULAR_OXYGEN) +
1781 AR_MASS * getDensity(ARGON) +
1782 H_MASS * getDensity(HYDROGEN) +
1783 N_MASS * getDensity(ATOMIC_NITROGEN));
1784 setDensity(TOTAL_MASS, tmd);
1785
1786 }
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809 void gtd7(final double alt) {
1810
1811
1812 final double altt = (alt > ZN2[0]) ? alt : ZN2[0];
1813 gts7(altt);
1814 if (alt >= ZN2[0]) {
1815 return;
1816 }
1817
1818
1819
1820
1821 final double r = PMA[2][0] * PAVGM[2];
1822 meso_tgn2[0] = meso_tgn1[1];
1823 meso_tn2[0] = meso_tn1[4];
1824 meso_tn2[1] = PMA[0][0] * PAVGM[0] / (1.0 - sw[20] * glob7s(PMA[0]));
1825 meso_tn2[2] = PMA[1][0] * PAVGM[1] / (1.0 - sw[20] * glob7s(PMA[1]));
1826 meso_tn2[3] = PMA[2][0] * PAVGM[2] / (1.0 - sw[20] * sw[22] * glob7s(PMA[2]));
1827 meso_tgn2[1] = PMA[9][0] * PAVGM[8] * (1.0 + sw[20] * sw[22] * glob7s(PMA[9])) *
1828 meso_tn2[3] * meso_tn2[3] / (r * r);
1829 meso_tn3[0] = meso_tn2[3];
1830
1831
1832
1833
1834 if (alt <= ZN3[0]) {
1835 final double q = PMA[6][0] * PAVGM[6];
1836 meso_tgn3[0] = meso_tgn2[1];
1837 meso_tn3[1] = PMA[3][0] * PAVGM[3] / (1.0 - sw[22] * glob7s(PMA[3]));
1838 meso_tn3[2] = PMA[4][0] * PAVGM[4] / (1.0 - sw[22] * glob7s(PMA[4]));
1839 meso_tn3[3] = PMA[5][0] * PAVGM[5] / (1.0 - sw[22] * glob7s(PMA[5]));
1840 meso_tn3[4] = PMA[6][0] * PAVGM[6] / (1.0 - sw[22] * glob7s(PMA[6]));
1841 meso_tgn3[1] = PMA[7][0] * PAVGM[7] * (1.0 + sw[22] * glob7s(PMA[7])) *
1842 meso_tn3[4] * meso_tn3[4] / (q * q);
1843
1844 }
1845
1846
1847 final double dmc = (alt > ZMIX) ? 1.0 - (ZN2[0] - alt) / (ZN2[0] - ZMIX) : 0.;
1848 final double dz28 = getDensity(MOLECULAR_NITROGEN);
1849
1850
1851 final double dm28m = dm28 * 1.0e+06;
1852 double dmr = dz28 / dm28m - 1.0;
1853 double dst = densm(alt, dm28m, PDM[2][4]) * (1.0 + dmr * dmc);
1854 setDensity(MOLECULAR_NITROGEN, dst);
1855
1856
1857 dmr = getDensity(HELIUM) / (dz28 * PDM[0][1]) - 1.0;
1858 dst = getDensity(MOLECULAR_NITROGEN) * PDM[0][1] * (1.0 + dmr * dmc);
1859 setDensity(HELIUM, dst);
1860
1861
1862 setDensity(ATOMIC_OXYGEN, 0.);
1863 setDensity(ANOMALOUS_OXYGEN, 0.);
1864
1865
1866 dmr = getDensity(MOLECULAR_OXYGEN) / (dz28 * PDM[3][1]) - 1.0;
1867 dst = getDensity(MOLECULAR_NITROGEN) * PDM[3][1] * (1.0 + dmr * dmc);
1868 setDensity(MOLECULAR_OXYGEN, dst);
1869
1870
1871 dmr = getDensity(ARGON) / (dz28 * PDM[4][1]) - 1.0;
1872 dst = getDensity(MOLECULAR_NITROGEN) * PDM[4][1] * (1.0 + dmr * dmc);
1873 setDensity(ARGON, dst);
1874
1875
1876 setDensity(HYDROGEN, 0.);
1877
1878
1879 setDensity(ATOMIC_NITROGEN, 0.);
1880
1881
1882 final double tmd = AMU * (HE_MASS * getDensity(HELIUM) +
1883 O_MASS * getDensity(ATOMIC_OXYGEN) +
1884 N2_MASS * getDensity(MOLECULAR_NITROGEN) +
1885 O2_MASS * getDensity(MOLECULAR_OXYGEN) +
1886 AR_MASS * getDensity(ARGON) +
1887 H_MASS * getDensity(HYDROGEN) +
1888 N_MASS * getDensity(ATOMIC_NITROGEN));
1889 setDensity(TOTAL_MASS, tmd);
1890
1891
1892 setTemperature(ALTITUDE, densm(alt, 1.0, 0));
1893
1894 }
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918 void gtd7d(final double alt) {
1919
1920
1921 gtd7(alt);
1922
1923
1924 final double dTot = getDensity(TOTAL_MASS) + AMU * O_MASS * getDensity(ANOMALOUS_OXYGEN);
1925 setDensity(TOTAL_MASS, dTot);
1926
1927 }
1928
1929
1930
1931
1932
1933
1934
1935
1936
1937
1938
1939
1940
1941
1942
1943
1944 void setDensity(final int index, final double d) {
1945 densities[index] = d;
1946 }
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956 void setTemperature(final int index, final double t) {
1957 temperatures[index] = t;
1958 }
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975 public double getDensity(final int index) {
1976 return densities[index];
1977 }
1978
1979
1980
1981
1982
1983 private double globe7(final double[] p) {
1984
1985 final double[] t = new double[14];
1986 final double cd32 = FastMath.cos(DAY_TO_RAD * (doy - p[31]));
1987 final double cd18 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[17]));
1988 final double cd14 = FastMath.cos(DAY_TO_RAD * (doy - p[13]));
1989 final double cd39 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[38]));
1990
1991
1992 final double df = f107 - f107a;
1993 final double dfa = f107a - FLUX_REF;
1994 t[0] = p[19] * df * (1.0 + p[59] * dfa) + p[20] * df * df + p[21] * dfa + p[29] * dfa * dfa;
1995
1996 final double f1 = 1.0 + (p[47] * dfa + p[19] * df + p[20] * df * df) * swc[1];
1997 final double f2 = 1.0 + (p[49] * dfa + p[19] * df + p[20] * df * df) * swc[1];
1998
1999
2000 t[1] = (p[1] * plg[0][2] + p[2] * plg[0][4] + p[22] * plg[0][6]) +
2001 (p[14] * plg[0][2]) * dfa * swc[1] + p[26] * plg[0][1];
2002
2003
2004 t[2] = p[18] * cd32;
2005
2006
2007 t[3] = (p[15] + p[16] * plg[0][2]) * cd18;
2008
2009
2010 t[4] = f1 * (p[9] * plg[0][1] + p[10] * plg[0][3]) * cd14;
2011
2012
2013 t[5] = p[37] * plg[0][1] * cd39;
2014
2015
2016 if (sw[7] != 0) {
2017 final double t71 = (p[11] * plg[1][2]) * cd14 * swc[5];
2018 final double t72 = (p[12] * plg[1][2]) * cd14 * swc[5];
2019 t[6] = f2 * ((p[3] * plg[1][1] + p[4] * plg[1][3] + p[27] * plg[1][5] + t71) * ctloc +
2020 (p[6] * plg[1][1] + p[7] * plg[1][3] + p[28] * plg[1][5] + t72) * stloc);
2021 }
2022
2023
2024 if (sw[8] != 0) {
2025 final double t81 = (p[23] * plg[2][3] + p[35] * plg[2][5]) * cd14 * swc[5];
2026 final double t82 = (p[33] * plg[2][3] + p[36] * plg[2][5]) * cd14 * swc[5];
2027 t[7] = f2 * ((p[5] * plg[2][2] + p[41] * plg[2][4] + t81) * c2tloc +
2028 (p[8] * plg[2][2] + p[42] * plg[2][4] + t82) * s2tloc);
2029 }
2030
2031
2032 if (sw[14] != 0) {
2033 t[13] = f2 * ((p[39] * plg[3][3] + (p[93] * plg[3][4] + p[46] * plg[3][6]) * cd14 * swc[5]) * s3tloc +
2034 (p[40] * plg[3][3] + (p[94] * plg[3][4] + p[48] * plg[3][6]) * cd14 * swc[5]) * c3tloc);
2035 }
2036
2037
2038 if (sw[9] == -1) {
2039 if (p[51] != 0) {
2040 final double exp1 = FastMath.exp(-10800.0 * FastMath.abs(p[51]) /
2041 (1.0 + p[138] * (LAT_REF - FastMath.abs(lat))));
2042 final double p24 = FastMath.max(p[24], 1.0e-4);
2043 apt = sg0(FastMath.min(exp1, 0.99999), p24, p[25]);
2044 t[8] = apt * (p[50] + p[96] * plg[0][2] + p[54] * plg[0][4] +
2045 (p[125] * plg[0][1] + p[126] * plg[0][3] + p[127] * plg[0][5]) * cd14 * swc[5] +
2046 (p[128] * plg[1][1] + p[129] * plg[1][3] + p[130] * plg[1][5]) * swc[7] *
2047 FastMath.cos(HOUR_TO_RAD * (hl - p[131])));
2048 }
2049 } else {
2050 final double apd = ap[0] - 4.0;
2051 final double p44 = (p[43] < 0.) ? 1.0E-5 : p[43];
2052 final double p45 = p[44];
2053 apdf = apd + (p45 - 1.0) * (apd + (FastMath.exp(-p44 * apd) - 1.0) / p44);
2054 if (sw[9] != 0) {
2055 t[8] = apdf * (p[32] + p[45] * plg[0][2] + p[34] * plg[0][4] +
2056 (p[100] * plg[0][1] + p[101] * plg[0][3] + p[102] * plg[0][5]) * cd14 * swc[5] +
2057 (p[121] * plg[1][1] + p[122] * plg[1][3] + p[123] * plg[1][5]) * swc[7] *
2058 FastMath.cos(HOUR_TO_RAD * (hl - p[124])));
2059 }
2060 }
2061
2062 if (sw[10] != 0) {
2063 final double lonr = DEG_TO_RAD * lon;
2064 final SinCos scLonr = FastMath.sinCos(lonr);
2065
2066 if (sw[11] != 0) {
2067 t[10] = (1.0 + p[80] * dfa * swc[1]) *
2068 ((p[64] * plg[1][2] + p[65] * plg[1][4] + p[66] * plg[1][6] +
2069 p[103] * plg[1][1] + p[104] * plg[1][3] + p[105] * plg[1][5] +
2070 (p[109] * plg[1][1] + p[110] * plg[1][3] + p[111] * plg[1][5]) * swc[5] * cd14) *
2071 scLonr.cos() +
2072 (p[90] * plg[1][2] + p[91] * plg[1][4] + p[92] * plg[1][6] +
2073 p[106] * plg[1][1] + p[107] * plg[1][3] + p[108] * plg[1][5] +
2074 (p[112] * plg[1][1] + p[113] * plg[1][3] + p[114] * plg[1][5]) * swc[5] * cd14) *
2075 scLonr.sin());
2076 }
2077
2078
2079 if (sw[12] != 0) {
2080 t[11] = (1.0 + p[95] * plg[0][1]) * (1.0 + p[81] * dfa * swc[1]) *
2081 (1.0 + p[119] * plg[0][1] * swc[5] * cd14) *
2082 (p[68] * plg[0][1] + p[69] * plg[0][3] + p[70] * plg[0][5]) *
2083 FastMath.cos(SEC_TO_RAD * (sec - p[71]));
2084 t[11] += swc[11] * (1.0 + p[137] * dfa * swc[1]) *
2085 (p[76] * plg[2][3] + p[77] * plg[2][5] + p[78] * plg[2][7]) *
2086 FastMath.cos(SEC_TO_RAD * (sec - p[79]) + 2.0 * lonr);
2087 }
2088
2089
2090 if (sw[13] != 0) {
2091 if (sw[9] == -1) {
2092 if (p[51] != 0.) {
2093 t[12] = apt * swc[11] * (1. + p[132] * plg[0][1]) *
2094 (p[52] * plg[1][2] + p[98] * plg[1][4] + p[67] * plg[1][6]) *
2095 FastMath.cos(DEG_TO_RAD * (lon - p[97])) +
2096 apt * swc[11] * swc[5] * cd14 *
2097 (p[133] * plg[1][1] + p[134] * plg[1][3] + p[135] * plg[1][5]) *
2098 FastMath.cos(DEG_TO_RAD * (lon - p[136])) +
2099 apt * swc[12] *
2100 (p[55] * plg[0][1] + p[56] * plg[0][3] + p[57] * plg[0][5]) *
2101 FastMath.cos(SEC_TO_RAD * (sec - p[58]));
2102 }
2103 } else {
2104 t[12] = apdf * swc[11] * (1.0 + p[120] * plg[0][1]) *
2105 ((p[60] * plg[1][2] + p[61] * plg[1][4] + p[62] * plg[1][6]) *
2106 FastMath.cos(DEG_TO_RAD * (lon - p[63]))) +
2107 apdf * swc[11] * swc[5] * cd14 *
2108 (p[115] * plg[1][1] + p[116] * plg[1][3] + p[117] * plg[1][5]) *
2109 FastMath.cos(DEG_TO_RAD * (lon - p[118])) +
2110 apdf * swc[12] *
2111 (p[83] * plg[0][1] + p[84] * plg[0][3] + p[85] * plg[0][5]) *
2112 FastMath.cos(SEC_TO_RAD * (sec - p[75]));
2113 }
2114 }
2115 }
2116
2117
2118 double tinf = p[30];
2119 for (int i = 0; i < 14; i++) {
2120 tinf += FastMath.abs(sw[i + 1]) * t[i];
2121 }
2122
2123
2124 return tinf;
2125
2126 }
2127
2128
2129
2130
2131
2132 private double glob7s(final double[] p) {
2133
2134 final double[] t = new double[14];
2135 final double cd32 = FastMath.cos(DAY_TO_RAD * (doy - p[31]));
2136 final double cd18 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[17]));
2137 final double cd14 = FastMath.cos(DAY_TO_RAD * (doy - p[13]));
2138 final double cd39 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[38]));
2139
2140
2141 t[0] = p[21] * (f107a - FLUX_REF);
2142
2143
2144 t[1] = p[1] * plg[0][2] + p[2] * plg[0][4] + p[22] * plg[0][6] +
2145 p[26] * plg[0][1] + p[14] * plg[0][3] + p[59] * plg[0][5];
2146
2147
2148 t[2] = (p[18] + p[47] * plg[0][2] + p[29] * plg[0][4]) * cd32;
2149
2150
2151 t[3] = (p[15] + p[16] * plg[0][2] + p[30] * plg[0][4]) * cd18;
2152
2153
2154 t[4] = (p[9] * plg[0][1] + p[10] * plg[0][3] + p[20] * plg[0][5]) * cd14;
2155
2156
2157 t[5] = (p[37] * plg[0][1]) * cd39;
2158
2159
2160 if (sw[7] != 0) {
2161 final double t71 = p[11] * plg[1][2] * cd14 * swc[5];
2162 final double t72 = p[12] * plg[1][2] * cd14 * swc[5];
2163 t[6] = (p[3] * plg[1][1] + p[4] * plg[1][3] + t71) * ctloc +
2164 (p[6] * plg[1][1] + p[7] * plg[1][3] + t72) * stloc;
2165 }
2166
2167
2168 if (sw[8] != 0) {
2169 final double t81 = (p[23] * plg[2][3] + p[35] * plg[2][5]) * cd14 * swc[5];
2170 final double t82 = (p[33] * plg[2][3] + p[36] * plg[2][5]) * cd14 * swc[5];
2171 t[7] = (p[5] * plg[2][2] + p[41] * plg[2][4] + t81) * c2tloc +
2172 (p[8] * plg[2][2] + p[42] * plg[2][4] + t82) * s2tloc;
2173 }
2174
2175
2176 if (sw[14] != 0) {
2177 t[13] = p[39] * plg[3][3] * s3tloc + p[40] * plg[3][3] * c3tloc;
2178 }
2179
2180
2181 if (sw[9] == 1) {
2182 t[8] = apdf * (p[32] + p[45] * plg[0][2] * swc[2]);
2183 } else if (sw[9] == -1) {
2184 t[8] = apt * (p[50] + p[96] * plg[0][2] * swc[2]);
2185 }
2186
2187
2188 if (!(sw[10] == 0 || sw[11] == 0)) {
2189 final double lonr = DEG_TO_RAD * lon;
2190 final SinCos scLonr = FastMath.sinCos(lonr);
2191 t[10] = (1.0 + plg[0][1] * (p[80] * swc[5] * FastMath.cos(DAY_TO_RAD * (doy - p[81])) +
2192 p[85] * swc[6] * FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[86]))) +
2193 p[83] * swc[3] * FastMath.cos(DAY_TO_RAD * (doy - p[84])) +
2194 p[87] * swc[4] * FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[88]))) *
2195 ((p[64] * plg[1][2] + p[65] * plg[1][4] + p[66] * plg[1][6] +
2196 p[74] * plg[1][1] + p[75] * plg[1][3] + p[76] * plg[1][5]) * scLonr.cos() +
2197 (p[90] * plg[1][2] + p[91] * plg[1][4] + p[92] * plg[1][6] +
2198 p[77] * plg[1][1] + p[78] * plg[1][3] + p[79] * plg[1][5]) * scLonr.sin());
2199 }
2200
2201
2202 double gl = 0;
2203 for (int i = 0; i < 14; i++) {
2204 gl += FastMath.abs(sw[i + 1]) * t[i];
2205 }
2206
2207
2208 return gl;
2209 }
2210
2211
2212
2213
2214
2215
2216
2217 private double sg0(final double ex, final double p24, final double p25) {
2218 final double g01 = g0(ap[1], p24, p25);
2219 final double g02 = g0(ap[2], p24, p25);
2220 final double g03 = g0(ap[3], p24, p25);
2221 final double g04 = g0(ap[4], p24, p25);
2222 final double g05 = g0(ap[5], p24, p25);
2223 final double g06 = g0(ap[6], p24, p25);
2224 final double ex2 = ex * ex;
2225 final double ex3 = ex * ex2;
2226 final double ex4 = ex2 * ex2;
2227 final double ex8 = ex4 * ex4;
2228 final double ex12 = ex4 * ex8;
2229 final double g234 = g02 * ex + g03 * ex2 + g04 * ex3;
2230 final double g56 = g05 * ex4 + g06 * ex12;
2231 final double ex19 = ex3 * ex4 * ex12;
2232 final double omex = 1.0 - ex;
2233 final double sumex = 1.0 + (1.0 - ex19) / omex * FastMath.sqrt(ex);
2234 return (g01 + (g234 + g56 * (1.0 - ex8) / omex)) / sumex;
2235 }
2236
2237
2238
2239
2240
2241
2242
2243 private double g0(final double apI, final double p24, final double p25) {
2244 final double am4 = apI - 4.0;
2245 return am4 + (p25 - 1.0) * (am4 + (FastMath.exp(-p24 * am4) - 1.0) / p24);
2246 }
2247
2248
2249
2250
2251
2252
2253
2254
2255 private double ccor(final double alt, final double r, final double h1, final double zh) {
2256 final double e = (alt - zh) / h1;
2257 if (e > 70.) {
2258 return 1.;
2259 } else if (e < -70.) {
2260 return FastMath.exp(r);
2261 } else {
2262 return FastMath.exp(r / (1.0 + FastMath.exp(e)));
2263 }
2264 }
2265
2266
2267
2268
2269
2270
2271
2272
2273
2274
2275 private double ccor2(final double alt, final double r,
2276 final double h1, final double zh, final double h2) {
2277 final double e1 = (alt - zh) / h1;
2278 final double e2 = (alt - zh) / h2;
2279 if (e1 > 70. || e2 > 70.) {
2280 return 1.;
2281 } else if (e1 < -70. && e2 < -70.) {
2282 return FastMath.exp(r);
2283 } else {
2284 final double ex1 = FastMath.exp(e1);
2285 final double ex2 = FastMath.exp(e2);
2286 return FastMath.exp(r / (1.0 + 0.5 * (ex1 + ex2)));
2287 }
2288 }
2289
2290
2291
2292
2293
2294
2295
2296 private double scalh(final double alt, final double xm, final double temp) {
2297
2298 final double denom = 1.0 + alt / rlat;
2299 final double galt = glat / (denom * denom);
2300 return R_GAS * temp / (galt * xm);
2301 }
2302
2303
2304
2305
2306
2307
2308
2309
2310
2311 private double dnet(final double dd, final double dm,
2312 final double zhm, final double xmm, final double xm) {
2313 if (!(dm > 0 && dd > 0)) {
2314 double ddd = dd;
2315 if (dd == 0 && dm == 0) {
2316 ddd = 1;
2317 }
2318 if (dm == 0) {
2319 return ddd;
2320 }
2321 if (dd == 0) {
2322 return dm;
2323 }
2324 }
2325
2326 final double a = zhm / (xmm - xm);
2327 final double ylog = a * FastMath.log(dm / dd);
2328 if (ylog < -10.) {
2329 return dd;
2330 } else if (ylog > 10.) {
2331 return dm;
2332 } else {
2333 return dd * FastMath.pow(1.0 + FastMath.exp(ylog), 1.0 / a);
2334 }
2335 }
2336
2337
2338
2339
2340
2341
2342
2343
2344
2345 private double splini(final double[] xa, final double[] ya, final double[] y2a, final double x) {
2346 final int n = xa.length;
2347 double yi = 0;
2348 int klo = 0;
2349 int khi = 1;
2350 while (x > xa[klo] && khi < n) {
2351 double xx = x;
2352 if (khi < n - 1) {
2353 xx = (x < xa[khi]) ? x : xa[khi];
2354 }
2355 final double h = xa[khi] - xa[klo];
2356 final double a = (xa[khi] - xx) / h;
2357 final double b = (xx - xa[klo]) / h;
2358 final double a2 = a * a;
2359 final double b2 = b * b;
2360 yi += ((1.0 - a2) * ya[klo] / 2.0 + b2 * ya[khi] / 2.0 +
2361 ((-(1.0 + a2 * a2) / 4.0 + a2 / 2.0) * y2a[klo] +
2362 (b2 * b2 / 4.0 - b2 / 2.0) * y2a[khi]) * h * h / 6.0) * h;
2363 klo++;
2364 khi++;
2365 }
2366 return yi;
2367 }
2368
2369
2370
2371
2372
2373
2374
2375
2376
2377 private double splint(final double[] xa, final double[] ya, final double[] y2a, final double x) {
2378 final int n = xa.length;
2379 int klo = 0;
2380 int khi = n - 1;
2381 while (khi - klo > 1) {
2382 final int k = (khi + klo) >>> 1;
2383 if (xa[k] > x) {
2384 khi = k;
2385 } else {
2386 klo = k;
2387 }
2388 }
2389 final double h = xa[khi] - xa[klo];
2390 final double a = (xa[khi] - x) / h;
2391 final double b = (x - xa[klo]) / h;
2392 return a * ya[klo] + b * ya[khi] +
2393 ((a * a * a - a) * y2a[klo] + (b * b * b - b) * y2a[khi]) * h * h / 6.0;
2394 }
2395
2396
2397
2398
2399
2400
2401
2402
2403
2404 private double[] spline(final double[] x, final double[] y, final double yp1, final double ypn) {
2405 final int n = x.length;
2406 final double[] y2 = new double[n];
2407 final double[] u = new double[n];
2408
2409 if (yp1 < 1e+30) {
2410 y2[0] = -0.5;
2411 u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - yp1);
2412 }
2413 for (int i = 1; i < n - 1; i++) {
2414 final double sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
2415 final double p = sig * y2[i - 1] + 2.0;
2416 y2[i] = (sig - 1.0) / p;
2417 u[i] = (6.0 * ((y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1])) /
2418 (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
2419 }
2420
2421 double qn = 0;
2422 double un = 0;
2423 if (ypn < 1e+30) {
2424 qn = 0.5;
2425 un = (3.0 / (x[n - 1] - x[n - 2])) * (ypn - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
2426 }
2427
2428 y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
2429 for (int k = n - 2; k >= 0; k--) {
2430 y2[k] = y2[k] * y2[k + 1] + u[k];
2431 }
2432
2433 return y2;
2434 }
2435
2436
2437
2438
2439
2440
2441
2442 private double densm(final double alt, final double d0, final double xm) {
2443
2444 double densm = d0;
2445
2446
2447 int mn = ZN2.length;
2448 double z = (alt > ZN2[mn - 1]) ? alt : ZN2[mn - 1];
2449
2450 double z1 = ZN2[0];
2451 double z2 = ZN2[mn - 1];
2452 double t1 = meso_tn2[0];
2453 double t2 = meso_tn2[mn - 1];
2454 double zg = zeta(z, z1);
2455 double zgdif = zeta(z2, z1);
2456
2457
2458 double[] xs = new double[mn];
2459 double[] ys = new double[mn];
2460 for (int k = 0; k < mn; k++) {
2461 xs[k] = zeta(ZN2[k], z1) / zgdif;
2462 ys[k] = 1.0 / meso_tn2[k];
2463 }
2464 final double qSM = (rlat + z2) / (rlat + z1);
2465 double yd1 = -meso_tgn2[0] / (t1 * t1) * zgdif;
2466 double yd2 = -meso_tgn2[1] / (t2 * t2) * zgdif * qSM * qSM;
2467
2468
2469 double[] y2out = spline(xs, ys, yd1, yd2);
2470 double x = zg / zgdif;
2471 double y = splint(xs, ys, y2out, x);
2472
2473
2474 double tz = 1.0 / y;
2475
2476 if (xm != 0.0) {
2477
2478 final double glb = galt(z1);
2479 final double gamm = xm * glb * zgdif / R_GAS;
2480
2481
2482 final double yi = splini(xs, ys, y2out, x);
2483 final double expl = FastMath.min(MIN_TEMP, gamm * yi);
2484
2485
2486 densm *= (t1 / tz) * FastMath.exp(-expl);
2487 }
2488
2489 if (alt > ZN3[0]) {
2490 return (xm == 0.0) ? tz : densm;
2491 }
2492
2493
2494 z = alt;
2495 mn = ZN3.length;
2496 z1 = ZN3[0];
2497 z2 = ZN3[mn - 1];
2498 t1 = meso_tn3[0];
2499 t2 = meso_tn3[mn - 1];
2500 zg = zeta(z, z1);
2501 zgdif = zeta(z2, z1);
2502
2503
2504 xs = new double[mn];
2505 ys = new double[mn];
2506 for (int k = 0; k < mn; k++) {
2507 xs[k] = zeta(ZN3[k], z1) / zgdif;
2508 ys[k] = 1.0 / meso_tn3[k];
2509 }
2510 final double qTS = (rlat + z2) / (rlat + z1);
2511 yd1 = -meso_tgn3[0] / (t1 * t1) * zgdif;
2512 yd2 = -meso_tgn3[1] / (t2 * t2) * zgdif * qTS * qTS;
2513
2514
2515 y2out = spline(xs, ys, yd1, yd2);
2516 x = zg / zgdif;
2517 y = splint(xs, ys, y2out, x);
2518
2519
2520 tz = 1.0 / y;
2521
2522 if (xm != 0.0) {
2523
2524 final double glb = galt(z1);
2525 final double gamm = xm * glb * zgdif / R_GAS;
2526
2527
2528 final double yi = splini(xs, ys, y2out, x);
2529 final double expl = FastMath.min(MIN_TEMP, gamm * yi);
2530
2531
2532 densm *= (t1 / tz) * FastMath.exp(-expl);
2533 }
2534
2535 return (xm == 0.0) ? tz : densm;
2536 }
2537
2538
2539
2540
2541
2542
2543
2544
2545
2546
2547
2548
2549 private double densu(final double alt, final double dlb, final double tinf,
2550 final double tlb, final double xm, final double alpha,
2551 final double zlb, final double s2) {
2552
2553 double z = (alt > ZN1[0]) ? alt : ZN1[0];
2554
2555
2556 final double zg2 = zeta(z, zlb);
2557
2558
2559 final double tt = tinf - (tinf - tlb) * FastMath.exp(-s2 * zg2);
2560 final double ta = tt;
2561 double tz = tt;
2562
2563 final int mn = ZN1.length;
2564 final double[] xs = new double[mn];
2565 final double[] ys = new double[mn];
2566 double x = 0.;
2567 double[] y2out = new double[mn];
2568 double zgdif = 0.;
2569 if (alt < ZN1[0]) {
2570
2571
2572 final double p = (rlat + zlb) / (rlat + ZN1[0]);
2573 final double dta = (tinf - ta) * s2 * p * p;
2574 meso_tgn1[0] = dta;
2575 meso_tn1[0] = ta;
2576 z = (alt > ZN1[mn - 1]) ? alt : ZN1[mn - 1];
2577
2578 final double t1 = meso_tn1[0];
2579 final double t2 = meso_tn1[mn - 1];
2580
2581 final double zg = zeta(z, ZN1[0]);
2582 zgdif = zeta(ZN1[mn - 1], ZN1[0]);
2583
2584 for (int k = 0; k < mn; k++) {
2585 xs[k] = zeta(ZN1[k], ZN1[0]) / zgdif;
2586 ys[k] = 1.0 / meso_tn1[k];
2587 }
2588
2589 final double q = (rlat + ZN1[mn - 1]) / (rlat + ZN1[0]);
2590 final double yd1 = -meso_tgn1[0] / (t1 * t1) * zgdif;
2591 final double yd2 = -meso_tgn1[1] / (t2 * t2) * zgdif * q * q;
2592
2593 y2out = spline(xs, ys, yd1, yd2);
2594 x = zg / zgdif;
2595 final double y = splint(xs, ys, y2out, x);
2596
2597 tz = 1.0 / y;
2598 }
2599
2600 if (xm == 0) {
2601 return tz;
2602 }
2603
2604
2605 double glb = galt(zlb);
2606 double gamma = xm * glb / (R_GAS * s2 * tinf);
2607 double expl = (tt <= 0) ? MIN_TEMP : FastMath.min(MIN_TEMP, FastMath.exp(-s2 * gamma * zg2));
2608 double densu = dlb * expl * FastMath.pow(tlb / tt, 1.0 + alpha + gamma);
2609
2610
2611 if (!Double.isFinite(densu)) {
2612 if (expl < MIN_TEMP) {
2613 densu = dlb * FastMath.exp(FastMath.log(tlb / tt) * (1.0 + alpha + gamma) - s2 * gamma * zg2);
2614 } else {
2615 throw new OrekitException( OrekitMessages.INFINITE_NRLMSISE00_DENSITY);
2616 }
2617 }
2618
2619
2620 if (alt < ZN1[0]) {
2621 glb = galt(ZN1[0]);
2622 gamma = xm * glb * zgdif / R_GAS;
2623
2624 expl = (tz <= 0) ? MIN_TEMP : FastMath.min(MIN_TEMP, gamma * splini(xs, ys, y2out, x));
2625
2626 densu *= FastMath.pow(meso_tn1[0] / tz, 1.0 + alpha) * FastMath.exp(-expl);
2627 }
2628
2629
2630 return densu;
2631 }
2632
2633
2634
2635
2636
2637 private double galt(final double alt) {
2638 final double r = 1.0 + alt / rlat;
2639 return glat / (r * r);
2640 }
2641
2642
2643
2644
2645
2646
2647 private double zeta(final double zz, final double zl) {
2648 return (zz - zl) * (rlat + zl) / (rlat + zz);
2649 }
2650
2651 }
2652
2653
2654
2655
2656
2657
2658
2659
2660
2661
2662
2663
2664
2665
2666
2667
2668
2669
2670
2671
2672
2673
2674
2675
2676
2677
2678
2679
2680
2681
2682
2683
2684
2685
2686
2687
2688
2689 public class FieldOutput<T extends CalculusFieldElement<T>> {
2690
2691
2692 private final Field<T> field;
2693
2694
2695 private final T zero;
2696
2697
2698 private final int doy;
2699
2700
2701 private final T sec;
2702
2703
2704 private final T lat;
2705
2706
2707 private final T lon;
2708
2709
2710 private final T hl;
2711
2712
2713 private final double f107a;
2714
2715
2716 private final double f107;
2717
2718
2719
2720
2721
2722
2723
2724
2725
2726
2727
2728 private final double[] ap;
2729
2730
2731 private final T glat;
2732
2733
2734 private final T rlat;
2735
2736
2737 private T dm28;
2738
2739
2740 private final T[][] plg;
2741
2742
2743 private final T ctloc;
2744
2745 private final T stloc;
2746
2747 private final T c2tloc;
2748
2749 private final T s2tloc;
2750
2751 private final T c3tloc;
2752
2753 private final T s3tloc;
2754
2755
2756 private double apdf;
2757
2758
2759 private T apt;
2760
2761
2762 private final T[] meso_tn1;
2763
2764
2765 private final T[] meso_tn2;
2766
2767
2768 private final T[] meso_tn3;
2769
2770
2771 private final T[] meso_tgn1;
2772
2773
2774 private final T[] meso_tgn2;
2775
2776
2777 private final T[] meso_tgn3;
2778
2779
2780 private final T[] densities;
2781
2782
2783 private final T[] temperatures;
2784
2785
2786
2787
2788
2789
2790
2791
2792
2793
2794
2795
2796
2797
2798
2799
2800
2801
2802
2803
2804 FieldOutput(final int doy, final T sec,
2805 final T lat, final T lon, final T hl,
2806 final double f107a, final double f107, final double[] ap) {
2807
2808 this.field = sec.getField();
2809 this.zero = field.getZero();
2810
2811 this.doy = doy;
2812 this.sec = sec;
2813 this.lat = lat;
2814 this.lon = lon;
2815 this.hl = hl;
2816 this.f107a = f107a;
2817 this.f107 = f107;
2818 this.ap = ap.clone();
2819
2820 this.plg = MathArrays.buildArray(field, 4, 8);
2821
2822 this.meso_tn1 = MathArrays.buildArray(field, ZN1.length);
2823 this.meso_tn2 = MathArrays.buildArray(field, ZN2.length);
2824 this.meso_tn3 = MathArrays.buildArray(field, ZN3.length);
2825 this.meso_tgn1 = MathArrays.buildArray(field, 2);
2826 this.meso_tgn2 = MathArrays.buildArray(field, 2);
2827 this.meso_tgn3 = MathArrays.buildArray(field, 2);
2828
2829 densities = MathArrays.buildArray(field, 9);
2830 temperatures = MathArrays.buildArray(field, 2);
2831
2832
2833 final T xlat = (sw[2] == 0) ? zero.newInstance(LAT_REF) : lat;
2834 final T c2 = xlat.multiply(2 * DEG_TO_RAD).cos();
2835 glat = c2.multiply(-0.0026373).add(1).multiply(G_REF);
2836 rlat = glat.multiply(2).divide(c2.multiply(2.27e-9).add(3.085462e-6)).multiply(1.e-5);
2837
2838
2839 final T latr = lat.multiply(DEG_TO_RAD);
2840
2841
2842 final FieldSinCos<T> scLatr = FastMath.sinCos(latr);
2843 final T c = scLatr.sin();
2844 final T s = scLatr.cos();
2845
2846 plg[0][1] = c;
2847 plg[0][2] = c.multiply( 3.0).multiply(plg[0][1]).subtract(1.0).divide(2.0);
2848 plg[0][3] = c.multiply( 5.0).multiply(plg[0][2]).subtract(plg[0][1].multiply(2.0)).divide(3.0);
2849 plg[0][4] = c.multiply( 7.0).multiply(plg[0][3]).subtract(plg[0][2].multiply(3.0)).divide(4.0);
2850 plg[0][5] = c.multiply( 9.0).multiply(plg[0][4]).subtract(plg[0][3].multiply(4.0)).divide(5.0);
2851 plg[0][6] = c.multiply(11.0).multiply(plg[0][5]).subtract(plg[0][4].multiply(5.0)).divide(6.0);
2852
2853 plg[1][1] = s;
2854 plg[1][2] = c.multiply( 3.0).multiply(plg[1][1]);
2855 plg[1][3] = c.multiply( 5.0).multiply(plg[1][2]).subtract(plg[1][1].multiply(3.0)).divide(2.0);
2856 plg[1][4] = c.multiply( 7.0).multiply(plg[1][3]).subtract(plg[1][2].multiply(4.0)).divide(3.0);
2857 plg[1][5] = c.multiply( 9.0).multiply(plg[1][4]).subtract(plg[1][3].multiply(5.0)).divide(4.0);
2858 plg[1][6] = c.multiply(11.0).multiply(plg[1][5]).subtract(plg[1][4].multiply(6.0)).divide(5.0);
2859
2860 plg[2][2] = s.multiply( 3.0).multiply(plg[1][1]);
2861 plg[2][3] = c.multiply( 5.0).multiply(plg[2][2]);
2862 plg[2][4] = c.multiply( 7.0).multiply(plg[2][3]).subtract(plg[2][2].multiply(5.0)).divide(2.0);
2863 plg[2][5] = c.multiply( 9.0).multiply(plg[2][4]).subtract(plg[2][3].multiply(6.0)).divide(3.0);
2864 plg[2][6] = c.multiply(11.0).multiply(plg[2][5]).subtract(plg[2][4].multiply(7.0)).divide(4.0);
2865 plg[2][7] = c.multiply(13.0).multiply(plg[2][6]).subtract(plg[2][5].multiply(8.0)).divide(5.0);
2866
2867 plg[3][3] = s.multiply( 5.0).multiply(plg[2][2]);
2868 plg[3][4] = c.multiply( 7.0).multiply(plg[3][3]);
2869 plg[3][5] = c.multiply( 9.0).multiply(plg[3][4]).subtract(plg[3][3].multiply(7.0)).divide(2.0);
2870 plg[3][6] = c.multiply(11.0).multiply(plg[3][5]).subtract(plg[3][4].multiply(8.0)).divide(3.0);
2871
2872
2873 if (!(sw[7] == 0 && sw[8] == 0 && sw[14] == 0)) {
2874 final T tloc = hl.multiply(HOUR_TO_RAD);
2875 final FieldSinCos<T> sc = FastMath.sinCos(tloc);
2876 final FieldSinCos<T> sc2 = FieldSinCos.sum(sc, sc);
2877 final FieldSinCos<T> sc3 = FieldSinCos.sum(sc, sc2);
2878 stloc = sc.sin();
2879 ctloc = sc.cos();
2880 s2tloc = sc2.sin();
2881 c2tloc = sc2.cos();
2882 s3tloc = sc3.sin();
2883 c3tloc = sc3.cos();
2884 } else {
2885 stloc = zero;
2886 ctloc = zero;
2887 s2tloc = zero;
2888 c2tloc = zero;
2889 s3tloc = zero;
2890 c3tloc = zero;
2891 }
2892
2893 }
2894
2895
2896
2897
2898
2899
2900
2901
2902
2903
2904
2905
2906
2907
2908
2909
2910
2911
2912
2913
2914
2915
2916
2917
2918
2919 void gts7(final T alt) {
2920
2921
2922 final double[] alpha = {-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0};
2923
2924 final double[] altl = {200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0};
2925
2926 final double xmm = PDM[2][4];
2927
2928
2929 T tinf = zero.newInstance(PTM[0] * PT[0]);
2930
2931 if (alt.getReal() > ZN1[0]) {
2932 tinf = tinf.multiply(globe7(PT).multiply(sw[16]).add(1));
2933 }
2934 setTemperature(EXOSPHERIC, tinf);
2935
2936
2937 T g0 = zero.newInstance(PTM[3] * PS[0]);
2938 if (alt.getReal() > ZN1[4]) {
2939 g0 = g0.multiply(globe7(PS).multiply(sw[19]).add(1));
2940 }
2941
2942
2943 T tlb = zero.newInstance(PTM[1] * PD[3][0]);
2944 tlb = tlb.multiply(globe7(PD[3]).multiply(sw[17]).add(1));
2945
2946
2947 final T s = g0.divide(tinf.subtract(tlb));
2948
2949
2950 meso_tn1[1] = zero.newInstance(PTM[6] * PTL[0][0]);
2951 meso_tn1[2] = zero.newInstance(PTM[2] * PTL[1][0]);
2952 meso_tn1[3] = zero.newInstance(PTM[7] * PTL[2][0]);
2953 meso_tn1[4] = zero.newInstance(PTM[4] * PTL[3][0]);
2954 meso_tgn1[1] = zero.newInstance(PTM[8] * PMA[8][0]);
2955 if (alt.getReal() < 300.0) {
2956 final double r = PTM[4] * PTL[3][0];
2957 meso_tn1[1] = meso_tn1[1].divide(glob7s(PTL[0]).multiply(sw[18] ).negate().add(1));
2958 meso_tn1[2] = meso_tn1[2].divide(glob7s(PTL[1]).multiply(sw[18] ).negate().add(1));
2959 meso_tn1[3] = meso_tn1[3].divide(glob7s(PTL[2]).multiply(sw[18] ).negate().add(1));
2960 meso_tn1[4] = meso_tn1[4].divide(glob7s(PTL[3]).multiply(sw[18] * sw[20]).negate().add(1));
2961 meso_tgn1[1] = meso_tgn1[1].multiply(glob7s(PMA[8]).multiply(sw[18] * sw[20]).add(1));
2962 meso_tgn1[1] = meso_tgn1[1].multiply(meso_tn1[4].multiply(meso_tn1[4]).divide(r * r));
2963 }
2964
2965
2966 setTemperature(ALTITUDE, densu(alt, zero.newInstance(1.0), tinf, tlb, 0, 0, PTM[5], s));
2967
2968
2969
2970 final T g28 = globe7(PD[2]).multiply(sw[21]);
2971
2972 final T db28 = g28.exp().multiply(PDM[2][0] * PD[2][0]);
2973
2974 T diffusiveDensity = densu(alt, db28, tinf, tlb, N2_MASS, alpha[2], PTM[5], s);
2975 setDensity(MOLECULAR_NITROGEN, diffusiveDensity);
2976
2977 final T zhf = lat.multiply(DEG_TO_RAD).sin().
2978 multiply(sw[5] * PDL[0][24] * FastMath.cos(DAY_TO_RAD * (doy - PT[13]))).
2979 add(1).
2980 multiply(PDL[1][24]);
2981
2982 final T zh28 = zhf.multiply(PDM[2][2]);
2983 final double zhm28 = PDM[2][3] * PDL[1][5];
2984
2985 final T b28 = densu(zh28, db28, tinf, tlb, N2_MASS - xmm, alpha[2] - 1.0, PTM[5], s);
2986 if (sw[15] != 0 && alt.getReal() <= altl[2]) {
2987
2988 dm28 = densu(alt, b28, tinf, tlb, xmm, alpha[2], PTM[5], s);
2989
2990 setDensity(MOLECULAR_NITROGEN, dnet(diffusiveDensity, dm28, zhm28, xmm, N2_MASS));
2991 } else {
2992 dm28 = zero;
2993 }
2994
2995
2996
2997 final T g4 = globe7(PD[0]).multiply(sw[21]);
2998
2999 final T db04 = g4.exp().multiply(PDM[0][0] * PD[0][0]);
3000
3001 diffusiveDensity = densu(alt, db04, tinf, tlb, HE_MASS, alpha[0], PTM[5], s);
3002 setDensity(HELIUM, diffusiveDensity);
3003 if (sw[15] != 0 && alt.getReal() < altl[0]) {
3004
3005 final double zh04 = PDM[0][2];
3006
3007 final T b04 = densu(zero.newInstance(zh04), db04, tinf, tlb, HE_MASS - xmm, alpha[0] - 1., PTM[5], s);
3008
3009 final T dm04 = densu(alt, b04, tinf, tlb, xmm, 0., PTM[5], s);
3010 final double zhm04 = zhm28;
3011
3012 diffusiveDensity = dnet(diffusiveDensity, dm04, zhm04, xmm, HE_MASS);
3013
3014 final T rl = b28.multiply(PDM[0][1]).divide(b04).log();
3015 final double zc04 = PDM[0][4] * PDL[1][0];
3016 final double hc04 = PDM[0][5] * PDL[1][1];
3017
3018 setDensity(HELIUM, diffusiveDensity.multiply(ccor(alt, rl, hc04, zc04)));
3019 }
3020
3021
3022
3023 final T g16 = globe7(PD[1]).multiply(sw[21]);
3024
3025 final T db16 = g16.exp().multiply(PDM[1][0] * PD[1][0]);
3026
3027 diffusiveDensity = densu(alt, db16, tinf, tlb, O_MASS, alpha[1], PTM[5], s);
3028 setDensity(ATOMIC_OXYGEN, diffusiveDensity);
3029 if (sw[15] != 0 && alt.getReal() < altl[1]) {
3030
3031 final double zh16 = PDM[1][2];
3032
3033 final T b16 = densu(zero.newInstance(zh16), db16, tinf, tlb, O_MASS - xmm, alpha[1] - 1.0, PTM[5], s);
3034
3035 final T dm16 = densu(alt, b16, tinf, tlb, xmm, 0., PTM[5], s);
3036 final double zhm16 = zhm28;
3037
3038 diffusiveDensity = dnet(diffusiveDensity, dm16, zhm16, xmm, O_MASS);
3039 final double rl = PDM[1][1] * PDL[1][16] * (1.0 + sw[1] * PDL[0][23] * (f107a - FLUX_REF));
3040 final double hc16 = PDM[1][5] * PDL[1][3];
3041 final double zc16 = PDM[1][4] * PDL[1][2];
3042 final double hc216 = PDM[1][5] * PDL[1][4];
3043 diffusiveDensity = diffusiveDensity.multiply(ccor2(alt, rl, hc16, zc16, hc216));
3044
3045 final double hcc16 = PDM[1][7] * PDL[1][13];
3046 final double zcc16 = PDM[1][6] * PDL[1][12];
3047 final double rc16 = PDM[1][3] * PDL[1][14];
3048
3049 setDensity(ATOMIC_OXYGEN, diffusiveDensity.multiply(ccor(alt, zero.newInstance(rc16), hcc16, zcc16)));
3050 }
3051
3052
3053
3054 final T g32 = globe7(PD[4]).multiply(sw[21]);
3055
3056 final T db32 = g32.exp().multiply(PDM[3][0] * PD[4][0]);
3057
3058 diffusiveDensity = densu(alt, db32, tinf, tlb, O2_MASS, alpha[3], PTM[5], s);
3059 setDensity(MOLECULAR_OXYGEN, diffusiveDensity);
3060 if (sw[15] != 0) {
3061 if (alt.getReal() <= altl[3]) {
3062
3063 final double zh32 = PDM[3][2];
3064
3065 final T b32 = densu(zero.newInstance(zh32), db32, tinf, tlb, O2_MASS - xmm, alpha[3] - 1., PTM[5], s);
3066
3067 final T dm32 = densu(alt, b32, tinf, tlb, xmm, 0., PTM[5], s);
3068 final double zhm32 = zhm28;
3069
3070 diffusiveDensity = dnet(diffusiveDensity, dm32, zhm32, xmm, O2_MASS);
3071
3072 final T rl = b28.multiply(PDM[3][1]).divide(b32).log();
3073 final double hc32 = PDM[3][5] * PDL[1][7];
3074 final double zc32 = PDM[3][4] * PDL[1][6];
3075 diffusiveDensity = diffusiveDensity.multiply(ccor(alt, rl, hc32, zc32));
3076 }
3077
3078 final double hcc32 = PDM[3][7] * PDL[1][22];
3079 final double hcc232 = PDM[3][7] * PDL[0][22];
3080 final double zcc32 = PDM[3][6] * PDL[1][21];
3081 final double rc32 = PDM[3][3] * PDL[1][23] * (1. + sw[1] * PDL[0][23] * (f107a - FLUX_REF));
3082
3083 setDensity(MOLECULAR_OXYGEN, diffusiveDensity.multiply(ccor2(alt, rc32, hcc32, zcc32, hcc232)));
3084 }
3085
3086
3087
3088 final T g40 = globe7(PD[5]).multiply(sw[21]);
3089
3090 final T db40 = g40.exp().multiply(PDM[4][0] * PD[5][0]);
3091
3092 diffusiveDensity = densu(alt, db40, tinf, tlb, AR_MASS, alpha[4], PTM[5], s);
3093 setDensity(ARGON, diffusiveDensity);
3094 if (sw[15] != 0 && alt.getReal() <= altl[4]) {
3095
3096 final double zh40 = PDM[4][2];
3097
3098 final T b40 = densu(zero.newInstance(zh40), db40, tinf, tlb, AR_MASS - xmm, alpha[4] - 1., PTM[5], s);
3099
3100 final T dm40 = densu(alt, b40, tinf, tlb, xmm, 0., PTM[5], s);
3101 final double zhm40 = zhm28;
3102
3103 diffusiveDensity = dnet(diffusiveDensity, dm40, zhm40, xmm, AR_MASS);
3104
3105 final T rl = b28.multiply(PDM[4][1]).divide(b40).log();
3106 final double hc40 = PDM[4][5] * PDL[1][9];
3107 final double zc40 = PDM[4][4] * PDL[1][8];
3108
3109 setDensity(ARGON, diffusiveDensity.multiply(ccor(alt, rl, hc40, zc40)));
3110 }
3111
3112
3113
3114 final T g1 = globe7(PD[6]).multiply(sw[21]);
3115
3116 final T db01 = g1.exp().multiply(PDM[5][0] * PD[6][0]);
3117
3118 diffusiveDensity = densu(alt, db01, tinf, tlb, H_MASS, alpha[6], PTM[5], s);
3119 setDensity(HYDROGEN, diffusiveDensity);
3120 if (sw[15] != 0 && alt.getReal() <= altl[6]) {
3121
3122 final double zh01 = PDM[5][2];
3123
3124 final T b01 = densu(zero.newInstance(zh01), db01, tinf, tlb, H_MASS - xmm, alpha[6] - 1., PTM[5], s);
3125
3126 final T dm01 = densu(alt, b01, tinf, tlb, xmm, 0., PTM[5], s);
3127 final double zhm01 = zhm28;
3128
3129 diffusiveDensity = dnet(diffusiveDensity, dm01, zhm01, xmm, H_MASS);
3130
3131 final T rl = b28.multiply(PDM[5][1] * FastMath.sqrt(PDL[1][17] * PDL[1][17])).divide(b01).log();
3132 final double hc01 = PDM[5][5] * PDL[1][11];
3133 final double zc01 = PDM[5][4] * PDL[1][10];
3134 diffusiveDensity = diffusiveDensity.multiply(ccor(alt, rl, hc01, zc01));
3135
3136 final double hcc01 = PDM[5][7] * PDL[1][19];
3137 final double zcc01 = PDM[5][6] * PDL[1][18];
3138 final double rc01 = PDM[5][3] * PDL[1][20];
3139
3140 setDensity(HYDROGEN, diffusiveDensity.multiply(ccor(alt, zero.newInstance(rc01), hcc01, zcc01)));
3141 }
3142
3143
3144
3145 final T g14 = globe7(PD[7]).multiply(sw[21]);
3146
3147 final T db14 = g14.exp().multiply(PDM[6][0] * PD[7][0]);
3148
3149 diffusiveDensity = densu(alt, db14, tinf, tlb, N_MASS, alpha[7], PTM[5], s);
3150 setDensity(ATOMIC_NITROGEN, diffusiveDensity);
3151 if (sw[15] != 0 && alt.getReal() <= altl[7]) {
3152
3153 final double zh14 = PDM[6][2];
3154
3155 final T b14 = densu(zero.newInstance(zh14), db14, tinf, tlb, N_MASS - xmm, alpha[7] - 1., PTM[5], s);
3156
3157 final T dm14 = densu(alt, b14, tinf, tlb, xmm, 0., PTM[5], s);
3158 final double zhm14 = zhm28;
3159
3160 diffusiveDensity = dnet(diffusiveDensity, dm14, zhm14, xmm, N_MASS);
3161
3162 final T rl = b28.multiply(PDM[6][1] * PDL[0][2]).divide(b14).log();
3163 final double hc14 = PDM[6][5] * PDL[0][1];
3164 final double zc14 = PDM[6][4] * PDL[0][0];
3165 diffusiveDensity = diffusiveDensity.multiply(ccor(alt, rl, hc14, zc14));
3166
3167 final double hcc14 = PDM[6][7] * PDL[0][4];
3168 final double zcc14 = PDM[6][6] * PDL[0][3];
3169 final double rc14 = PDM[6][3] * PDL[0][5];
3170
3171 setDensity(ATOMIC_NITROGEN, diffusiveDensity.multiply(ccor(alt, zero.newInstance(rc14), hcc14, zcc14)));
3172 }
3173
3174
3175 final T g16h = globe7(PD[8]).multiply(sw[21]);
3176 final T db16h = g16h.exp().multiply(PDM[7][0] * PD[8][0]);
3177 final double tho = PDM[7][9] * PDL[0][6];
3178 diffusiveDensity = densu(alt, db16h, zero.newInstance(tho), zero.newInstance(tho), O_MASS, alpha[8], PTM[5], s);
3179 final double zsht = PDM[7][5];
3180 final double zmho = PDM[7][4];
3181 final T zsho = scalh(zmho, O_MASS, tho);
3182 diffusiveDensity = diffusiveDensity.multiply(alt.negate().add(zmho).divide(zsht).exp().subtract(1).multiply(-zsht).divide(zsho).exp());
3183 setDensity(ANOMALOUS_OXYGEN, diffusiveDensity);
3184
3185
3186 for (int i = 0; i < 9; i++) {
3187 setDensity(i, getDensity(i).multiply(1.0e+06));
3188 }
3189
3190
3191 final T tmd = getDensity(HELIUM) .multiply(HE_MASS).
3192 add(getDensity(ATOMIC_OXYGEN) .multiply( O_MASS)).
3193 add(getDensity(MOLECULAR_NITROGEN).multiply(N2_MASS)).
3194 add(getDensity(MOLECULAR_OXYGEN) .multiply(O2_MASS)).
3195 add(getDensity(ARGON) .multiply(AR_MASS)).
3196 add(getDensity(HYDROGEN) .multiply( H_MASS)).
3197 add(getDensity(ATOMIC_NITROGEN) .multiply( N_MASS)).
3198 multiply(AMU);
3199 setDensity(TOTAL_MASS, tmd);
3200
3201 }
3202
3203
3204
3205
3206
3207
3208
3209
3210
3211
3212
3213
3214
3215
3216
3217
3218
3219
3220
3221
3222
3223
3224 void gtd7(final T alt) {
3225
3226
3227 final T altt = (alt.getReal() > ZN2[0]) ? alt : zero.newInstance(ZN2[0]);
3228 gts7(altt);
3229 if (alt.getReal() >= ZN2[0]) {
3230 return;
3231 }
3232
3233
3234
3235
3236 final double r = PMA[2][0] * PAVGM[2];
3237 meso_tgn2[0] = meso_tgn1[1];
3238 meso_tn2[0] = meso_tn1[4];
3239 meso_tn2[1] = glob7s(PMA[0]).multiply(sw[20] ).negate().add(1).reciprocal().multiply(PMA[0][0] * PAVGM[0]);
3240 meso_tn2[2] = glob7s(PMA[1]).multiply(sw[20] ).negate().add(1).reciprocal().multiply(PMA[1][0] * PAVGM[1]);
3241 meso_tn2[3] = glob7s(PMA[2]).multiply(sw[20] * sw[22]).negate().add(1).reciprocal().multiply(PMA[2][0] * PAVGM[2]);
3242 meso_tgn2[1] = glob7s(PMA[9]).multiply(sw[20] * sw[22]).add(1).multiply(PMA[9][0] * PAVGM[8]).
3243 multiply(meso_tn2[3]).multiply(meso_tn2[3]).divide(r * r);
3244 meso_tn3[0] = meso_tn2[3];
3245
3246
3247
3248
3249 if (alt.getReal() <= ZN3[0]) {
3250 final double q = PMA[6][0] * PAVGM[6];
3251 meso_tgn3[0] = meso_tgn2[1];
3252 meso_tn3[1] = glob7s(PMA[3]).multiply(sw[22]).negate().add(1).reciprocal().multiply(PMA[3][0] * PAVGM[3]);
3253 meso_tn3[2] = glob7s(PMA[4]).multiply(sw[22]).negate().add(1).reciprocal().multiply(PMA[4][0] * PAVGM[4]);
3254 meso_tn3[3] = glob7s(PMA[5]).multiply(sw[22]).negate().add(1).reciprocal().multiply(PMA[5][0] * PAVGM[5]);
3255 meso_tn3[4] = glob7s(PMA[6]).multiply(sw[22]).negate().add(1).reciprocal().multiply(PMA[6][0] * PAVGM[6]);
3256 meso_tgn3[1] = glob7s(PMA[7]).multiply(sw[22]) .add(1).multiply(PMA[7][0] * PAVGM[7]).
3257 multiply(meso_tn3[4]).multiply(meso_tn3[4]).divide(q * q);
3258
3259 }
3260
3261
3262 final T dmc = (alt.getReal() > ZMIX) ?
3263 alt.subtract(ZN2[0]).divide(ZN2[0] - ZMIX).add(1) :
3264 zero;
3265 final T dz28 = getDensity(MOLECULAR_NITROGEN);
3266
3267
3268 final T dm28m = dm28.multiply(1.0e+06);
3269 T dmr = dz28.divide(dm28m).subtract(1);
3270 T dst = densm(alt, dm28m, PDM[2][4]).multiply(dmr.multiply(dmc).add(1));
3271 setDensity(MOLECULAR_NITROGEN, dst);
3272
3273
3274 dmr = getDensity(HELIUM).divide(dz28.multiply(PDM[0][1])).subtract(1);
3275 dst = getDensity(MOLECULAR_NITROGEN).multiply(PDM[0][1]).multiply(dmr.multiply(dmc).add(1));
3276 setDensity(HELIUM, dst);
3277
3278
3279 setDensity(ATOMIC_OXYGEN, zero);
3280 setDensity(ANOMALOUS_OXYGEN, zero);
3281
3282
3283 dmr = getDensity(MOLECULAR_OXYGEN).divide(dz28.multiply(PDM[3][1])).subtract(1);
3284 dst = getDensity(MOLECULAR_NITROGEN).multiply(PDM[3][1]).multiply(dmr.multiply(dmc).add(1));
3285 setDensity(MOLECULAR_OXYGEN, dst);
3286
3287
3288 dmr = getDensity(ARGON).divide(dz28.multiply(PDM[4][1])).subtract(1);
3289 dst = getDensity(MOLECULAR_NITROGEN).multiply(PDM[4][1]).multiply(dmr.multiply(dmc).add(1));
3290 setDensity(ARGON, dst);
3291
3292
3293 setDensity(HYDROGEN, zero);
3294
3295
3296 setDensity(ATOMIC_NITROGEN, zero);
3297
3298
3299 final T tmd = getDensity(HELIUM) .multiply(HE_MASS).
3300 add(getDensity(ATOMIC_OXYGEN) .multiply( O_MASS)).
3301 add(getDensity(MOLECULAR_NITROGEN).multiply(N2_MASS)).
3302 add(getDensity(MOLECULAR_OXYGEN) .multiply(O2_MASS)).
3303 add(getDensity(ARGON) .multiply(AR_MASS)).
3304 add(getDensity(HYDROGEN) .multiply( H_MASS)).
3305 add(getDensity(ATOMIC_NITROGEN) .multiply( N_MASS)).
3306 multiply(AMU);
3307 setDensity(TOTAL_MASS, tmd);
3308
3309
3310 setTemperature(ALTITUDE, densm(alt, field.getOne(), 0));
3311
3312 }
3313
3314
3315
3316
3317
3318
3319
3320
3321
3322
3323
3324
3325
3326
3327
3328
3329
3330
3331
3332
3333
3334
3335
3336 void gtd7d(final T alt) {
3337
3338
3339 gtd7(alt);
3340
3341
3342 final T dTot = getDensity(TOTAL_MASS).add(getDensity(ANOMALOUS_OXYGEN).multiply( AMU * O_MASS));
3343 setDensity(TOTAL_MASS, dTot);
3344
3345 }
3346
3347
3348
3349
3350
3351
3352
3353
3354
3355
3356
3357
3358
3359
3360
3361
3362 void setDensity(final int index, final T d) {
3363 densities[index] = d;
3364 }
3365
3366
3367
3368
3369
3370
3371
3372
3373
3374 void setTemperature(final int index, final T t) {
3375 temperatures[index] = t;
3376 }
3377
3378
3379
3380
3381
3382
3383
3384
3385
3386
3387
3388
3389
3390
3391
3392
3393 public T getDensity(final int index) {
3394 return densities[index];
3395 }
3396
3397
3398
3399
3400
3401 private T globe7(final double[] p) {
3402
3403 final T[] t = MathArrays.buildArray(field, 14);
3404 final double cd32 = FastMath.cos(DAY_TO_RAD * (doy - p[31]));
3405 final double cd18 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[17]));
3406 final double cd14 = FastMath.cos(DAY_TO_RAD * (doy - p[13]));
3407 final double cd39 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[38]));
3408
3409
3410 final double df = f107 - f107a;
3411 final double dfa = f107a - FLUX_REF;
3412 t[0] = zero.newInstance(p[19] * df * (1.0 + p[59] * dfa) +
3413 p[20] * df * df +
3414 p[21] * dfa +
3415 p[29] * dfa * dfa);
3416
3417 final double f1 = 1.0 + (p[47] * dfa + p[19] * df + p[20] * df * df) * swc[1];
3418 final double f2 = 1.0 + (p[49] * dfa + p[19] * df + p[20] * df * df) * swc[1];
3419
3420
3421 t[1] = plg[0][2].multiply(p[ 1]).
3422 add(plg[0][4].multiply(p[ 2])).
3423 add(plg[0][6].multiply(p[22])).
3424 add(plg[0][2].multiply(p[14] * dfa * swc[1])).
3425 add(plg[0][1].multiply(p[26]));
3426
3427
3428 t[2] = zero.newInstance(p[18] * cd32);
3429
3430
3431 t[3] = plg[0][2].multiply(p[16]).add(p[15]).multiply(cd18);
3432
3433
3434 t[4] = plg[0][1].multiply(p[9]).add(plg[0][3].multiply(p[10])).multiply(f1 * cd14);
3435
3436
3437 t[5] = plg[0][1].multiply(p[37] * cd39);
3438
3439
3440 if (sw[7] != 0) {
3441 final T t71 = plg[1][2].multiply(p[11] * cd14 * swc[5]);
3442 final T t72 = plg[1][2].multiply(p[12] * cd14 * swc[5]);
3443 t[6] = plg[1][1].multiply(p[3]).add(plg[1][3].multiply(p[4])).add(plg[1][5].multiply(p[27])).add(t71).multiply(ctloc).
3444 add(plg[1][1].multiply(p[6]).add(plg[1][3].multiply(p[7])).add(plg[1][5].multiply(p[28])).add(t72).multiply(stloc)).
3445 multiply(f2);
3446 }
3447
3448
3449 if (sw[8] != 0) {
3450 final T t81 = plg[2][3].multiply(p[23]).add(plg[2][5].multiply(p[35])).multiply(cd14 * swc[5]);
3451 final T t82 = plg[2][3].multiply(p[33]).add(plg[2][5].multiply(p[36])).multiply(cd14 * swc[5]);
3452 t[7] = plg[2][2].multiply(p[5]).add(plg[2][4].multiply(p[41])).add(t81).multiply(c2tloc).
3453 add(plg[2][2].multiply(p[8]).add(plg[2][4].multiply(p[42])).add(t82).multiply(s2tloc)).
3454 multiply(f2);
3455 }
3456
3457
3458 if (sw[14] != 0) {
3459 t[13] = plg[3][3].multiply(p[39]).add(plg[3][4].multiply(p[93]).add(plg[3][6].multiply(p[46])).multiply(cd14 * swc[5])).multiply(s3tloc).
3460 add(plg[3][3].multiply(p[40]).add(plg[3][4].multiply(p[94]).add(plg[3][6].multiply(p[48])).multiply(cd14 * swc[5])).multiply(c3tloc)).
3461 multiply(f2);
3462 }
3463
3464
3465 if (sw[9] == -1) {
3466 if (p[51] != 0) {
3467 final T exp1 = lat.abs().negate().add(LAT_REF).multiply(p[138]).add(1).
3468 reciprocal().multiply(-10800.0 * FastMath.abs(p[51])).
3469 exp();
3470 final double p24 = FastMath.max(p[24], 1.0e-4);
3471 apt = sg0(min(0.99999, exp1), p24, p[25]);
3472 t[8] = plg[0][2].multiply(p[96]).add(plg[0][4].multiply(p[54])).add(p[50]).
3473 add((plg[0][1].multiply(p[125]).add(plg[0][3].multiply(p[126])).add(plg[0][5].multiply(p[127]))).multiply(cd14 * swc[5])).
3474 add((plg[1][1].multiply(p[128]).add(plg[1][3].multiply(p[129])).add(plg[1][5].multiply(p[130]))).multiply(swc[7]).multiply(hl.subtract(p[131]).multiply(HOUR_TO_RAD).cos())).
3475 multiply(apt);
3476 }
3477 } else {
3478 final double apd = ap[0] - 4.0;
3479 final double p44 = (p[43] < 0.) ? 1.0E-5 : p[43];
3480 final double p45 = p[44];
3481 apdf = apd + (p45 - 1.0) * (apd + (FastMath.exp(-p44 * apd) - 1.0) / p44);
3482 if (sw[9] != 0) {
3483 t[8] = plg[0][2].multiply(p[45]).add(plg[0][4].multiply(p[34])).add(p[32]).
3484 add((plg[0][1].multiply(p[100]).add(plg[0][3].multiply(p[101])).add(plg[0][5].multiply(p[102]))).multiply(cd14 * swc[5])).
3485 add((plg[1][1].multiply(p[121]).add(plg[1][3].multiply(p[122])).add(plg[1][5].multiply(p[123]))).multiply(swc[7]).multiply(hl.subtract(p[124]).multiply(HOUR_TO_RAD).cos())).
3486 multiply(apdf);
3487 }
3488 }
3489
3490 if (sw[10] != 0) {
3491 final T lonr = lon.multiply(DEG_TO_RAD);
3492 final FieldSinCos<T> scLonr = FastMath.sinCos(lonr);
3493
3494 if (sw[11] != 0) {
3495 t[10] = plg[1][2].multiply(p[ 64]) .add(plg[1][4].multiply(p[ 65])).add(plg[1][6].multiply(p[ 66])).
3496 add(plg[1][1].multiply(p[103])).add(plg[1][3].multiply(p[104])).add(plg[1][5].multiply(p[105])).
3497 add((plg[1][1].multiply(p[109])).add(plg[1][3].multiply(p[110])).add(plg[1][5].multiply(p[111])).multiply(swc[5] * cd14)).
3498 multiply(scLonr.cos()).
3499 add( plg[1][2].multiply(p[ 90]) .add(plg[1][4].multiply(p[ 91])).add(plg[1][6].multiply(p[ 92])).
3500 add(plg[1][1].multiply(p[106])).add(plg[1][3].multiply(p[107])).add(plg[1][5].multiply(p[108])).
3501 add((plg[1][1].multiply(p[112])).add(plg[1][3].multiply(p[113])).add(plg[1][5].multiply(p[114])).multiply(swc[5] * cd14)).
3502 multiply(scLonr.sin())).
3503 multiply(1.0 + p[80] * dfa * swc[1]);
3504 }
3505
3506
3507 if (sw[12] != 0) {
3508 t[11] = plg[0][1].multiply(p[95]).add(1).multiply(1.0 + p[81] * dfa * swc[1]).
3509 multiply(plg[0][1].multiply(p[119] * swc[5] * cd14).add(1)).
3510 multiply(plg[0][1].multiply(p[68]).add(plg[0][3].multiply(p[69])).add(plg[0][5].multiply(p[70]))).
3511 multiply(sec.subtract(p[71]).multiply(SEC_TO_RAD).cos());
3512 t[11] = t[11].
3513 add(plg[2][3].multiply(p[76]).add(plg[2][5].multiply(p[77])).add(plg[2][7].multiply(p[78])).
3514 multiply(swc[11] * (1.0 + p[137] * dfa * swc[1])).
3515 multiply(sec.subtract(p[79]).multiply(SEC_TO_RAD).add(lonr.multiply(2)).cos()));
3516 }
3517
3518
3519 if (sw[13] != 0) {
3520 if (sw[9] == -1) {
3521 if (p[51] != 0.) {
3522 t[12] = apt.multiply(swc[11]).multiply(plg[0][1].multiply(p[132]).add(1)).
3523 multiply(plg[1][2].multiply(p[52]).add(plg[1][4].multiply(p[98])).add(plg[1][6].multiply(p[67]))).
3524 multiply(lon.subtract(p[97]).multiply(DEG_TO_RAD).cos()).
3525 add(apt.multiply(swc[11] * swc[5] * cd14).
3526 multiply(plg[1][1].multiply(p[133]).add(plg[1][3].multiply(p[134])).add(plg[1][5].multiply(p[135]))).
3527 multiply(lon.subtract(p[136]).multiply(DEG_TO_RAD).cos())).
3528 add(apt.multiply(swc[12]).
3529 multiply(plg[0][1].multiply(p[55]).add(plg[0][3].multiply(p[56])).add(plg[0][5].multiply(p[57]))).
3530 multiply(sec.subtract(p[58]).multiply(SEC_TO_RAD).cos()));
3531 }
3532 } else {
3533 t[12] = plg[0][1].multiply(p[120]).add(1).multiply(apdf * swc[11]).
3534 multiply(plg[1][2].multiply(p[60]).add(plg[1][4].multiply(p[61])).add(plg[1][6].multiply(p[62]))).
3535 multiply(lon.subtract(p[63]).multiply(DEG_TO_RAD).cos()).
3536 add(plg[1][1].multiply(p[115]).add(plg[1][3].multiply(p[116])).add(plg[1][5].multiply(p[117])).
3537 multiply(apdf * swc[11] * swc[5] * cd14).
3538 multiply(lon.subtract(p[118]).multiply(DEG_TO_RAD).cos())).
3539 add(plg[0][1].multiply(p[83]).add(plg[0][3].multiply(p[84])).add(plg[0][5].multiply(p[85])).
3540 multiply(apdf * swc[12]).
3541 multiply(sec.subtract(p[75]).multiply(SEC_TO_RAD).cos()));
3542 }
3543 }
3544 }
3545
3546
3547 T tinf = zero.newInstance(p[30]);
3548 for (int i = 0; i < 14; i++) {
3549 tinf = tinf.add(t[i].multiply(FastMath.abs(sw[i + 1])));
3550 }
3551
3552
3553 return tinf;
3554
3555 }
3556
3557
3558
3559
3560
3561 private T glob7s(final double[] p) {
3562
3563 final T[] t = MathArrays.buildArray(field, 14);
3564 final double cd32 = FastMath.cos(DAY_TO_RAD * (doy - p[31]));
3565 final double cd18 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[17]));
3566 final double cd14 = FastMath.cos(DAY_TO_RAD * (doy - p[13]));
3567 final double cd39 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[38]));
3568
3569
3570 t[0] = zero.newInstance(p[21] * (f107a - FLUX_REF));
3571
3572
3573 t[1] = plg[0][2].multiply(p[1]).
3574 add(plg[0][4].multiply(p[2])).
3575 add(plg[0][6].multiply(p[22])).
3576 add(plg[0][1].multiply(p[26])).
3577 add(plg[0][3].multiply(p[14])).
3578 add(plg[0][5].multiply(p[59]));
3579
3580
3581 t[2] = plg[0][2].multiply(p[47]).add(plg[0][4].multiply(p[29])).add(p[18]).multiply(cd32);
3582
3583
3584 t[3] = plg[0][2].multiply(p[16]).add(plg[0][4].multiply(p[30])).add(p[15]).multiply(cd18);
3585
3586
3587 t[4] = plg[0][1].multiply(p[9]).add(plg[0][3].multiply(p[10])).add(plg[0][5].multiply(p[20])).multiply(cd14);
3588
3589
3590 t[5] = plg[0][1].multiply(p[37]).multiply(cd39);
3591
3592
3593 if (sw[7] != 0) {
3594 final T t71 = plg[1][2].multiply(p[11]).multiply(cd14 * swc[5]);
3595 final T t72 = plg[1][2].multiply(p[12]).multiply(cd14 * swc[5]);
3596 t[6] = plg[1][1].multiply(p[3]).add(plg[1][3].multiply(p[4])).add(t71).multiply(ctloc).
3597 add(plg[1][1].multiply(p[6]).add(plg[1][3].multiply(p[7])).add(t72).multiply(stloc));
3598 }
3599
3600
3601 if (sw[8] != 0) {
3602 final T t81 = plg[2][3].multiply(p[23]).add(plg[2][5].multiply(p[35])).multiply(cd14 * swc[5]);
3603 final T t82 = plg[2][3].multiply(p[33]).add(plg[2][5].multiply(p[36])).multiply(cd14 * swc[5]);
3604 t[7] = plg[2][2].multiply(p[5]).add(plg[2][4].multiply(p[41])).add(t81).multiply(c2tloc).
3605 add(plg[2][2].multiply(p[8]).add(plg[2][4].multiply(p[42])).add(t82).multiply(s2tloc));
3606 }
3607
3608
3609 if (sw[14] != 0) {
3610 t[13] = plg[3][3].multiply(p[39]).multiply(s3tloc).add(plg[3][3].multiply(p[40]).multiply(c3tloc));
3611 }
3612
3613
3614 if (sw[9] == 1) {
3615 t[8] = plg[0][2].multiply(p[45] * swc[2]).add(p[32]).multiply(apdf);
3616 } else if (sw[9] == -1) {
3617 t[8] = plg[0][2].multiply(p[96] * swc[2]).add(p[50]).multiply(apt);
3618 }
3619
3620
3621 if (!(sw[10] == 0 || sw[11] == 0)) {
3622 final T lonr = lon.multiply(DEG_TO_RAD);
3623 final FieldSinCos<T> scLonr = FastMath.sinCos(lonr);
3624 t[10] = plg[0][1].multiply(p[80] * swc[5] * FastMath.cos(DAY_TO_RAD * (doy - p[81])) +
3625 p[85] * swc[6] * FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[86]))).
3626 add(1.0 +
3627 p[83] * swc[3] * FastMath.cos(DAY_TO_RAD * (doy - p[84])) +
3628 p[87] * swc[4] * FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[88]))).
3629 multiply( plg[1][2].multiply(p[64]).
3630 add(plg[1][4].multiply(p[65])).
3631 add(plg[1][6].multiply(p[66])).
3632 add(plg[1][1].multiply(p[74])).
3633 add(plg[1][3].multiply(p[75])).
3634 add(plg[1][5].multiply(p[76])).multiply(scLonr.cos()).
3635 add( plg[1][2].multiply(p[90]).
3636 add(plg[1][4].multiply(p[91])).
3637 add(plg[1][6].multiply(p[92])).
3638 add(plg[1][1].multiply(p[77])).
3639 add(plg[1][3].multiply(p[78])).
3640 add(plg[1][5].multiply(p[79])).multiply(scLonr.sin())));
3641 }
3642
3643
3644 T gl = zero;
3645 for (int i = 0; i < 14; i++) {
3646 gl = gl.add(t[i].multiply(FastMath.abs(sw[i + 1])));
3647 }
3648
3649
3650 return gl;
3651 }
3652
3653
3654
3655
3656
3657
3658
3659 private T sg0(final T ex, final double p24, final double p25) {
3660 final double g01 = g0(ap[1], p24, p25);
3661 final double g02 = g0(ap[2], p24, p25);
3662 final double g03 = g0(ap[3], p24, p25);
3663 final double g04 = g0(ap[4], p24, p25);
3664 final double g05 = g0(ap[5], p24, p25);
3665 final double g06 = g0(ap[6], p24, p25);
3666 final T ex2 = ex.square();
3667 final T ex3 = ex.multiply(ex2);
3668 final T ex4 = ex2.square();
3669 final T ex8 = ex4.square();
3670 final T ex12 = ex4.multiply(ex8);
3671 final T g234 = ex.multiply(g02).add(ex2.multiply(g03)).add(ex3.multiply(g04));
3672 final T g56 = ex4.multiply(g05).add(ex12.multiply(g06));
3673 final T ex19 = ex3.multiply(ex4).multiply(ex12);
3674 final T omex = ex.negate().add(1);
3675 final T sumex = ex19.negate().add(1).divide(omex).multiply(ex.sqrt()).add(1);
3676 return ex8.negate().add(1).multiply(g56).divide(omex).add(g234).add(g01).divide(sumex);
3677 }
3678
3679
3680
3681
3682
3683
3684
3685 private double g0(final double apI, final double p24, final double p25) {
3686 final double am4 = apI - 4.0;
3687 return am4 + (p25 - 1.0) * (am4 + (FastMath.exp(-p24 * am4) - 1.0) / p24);
3688 }
3689
3690
3691
3692
3693
3694
3695
3696
3697 private T ccor(final T alt, final T r, final double h1, final double zh) {
3698 final T e = alt.subtract(zh).divide(h1);
3699 if (e.getReal() > 70.) {
3700 return field.getOne();
3701 } else if (e.getReal() < -70.) {
3702 return r.exp();
3703 } else {
3704 return r.divide(e.exp().add(1)).exp();
3705 }
3706 }
3707
3708
3709
3710
3711
3712
3713
3714
3715
3716
3717 private T ccor2(final T alt, final double r, final double h1, final double zh, final double h2) {
3718 final T e1 = alt.subtract(zh).divide(h1);
3719 final T e2 = alt.subtract(zh).divide(h2);
3720 if (e1.getReal() > 70. || e2.getReal() > 70.) {
3721 return field.getOne();
3722 } else if (e1.getReal() < -70. && e2.getReal() < -70.) {
3723 return zero.newInstance(FastMath.exp(r));
3724 } else {
3725 final T ex1 = e1.exp();
3726 final T ex2 = e2.exp();
3727 return ex1.add(ex2).multiply(0.5).add(1).reciprocal().multiply(r).exp();
3728 }
3729 }
3730
3731
3732
3733
3734
3735
3736
3737 private T scalh(final double alt, final double xm, final double temp) {
3738
3739 final T denom = rlat.reciprocal().multiply(alt).add(1);
3740 final T galt = glat.divide(denom.square());
3741 return galt.reciprocal().multiply(R_GAS * temp / xm);
3742 }
3743
3744
3745
3746
3747
3748
3749
3750
3751
3752 private T dnet(final T dd, final T dm, final double zhm, final double xmm, final double xm) {
3753 if (!(dm.getReal() > 0 && dd.getReal() > 0)) {
3754 T ddd = dd;
3755 if (dd.getReal() == 0 && dm.getReal() == 0) {
3756 ddd = field.getOne();
3757 }
3758 if (dm.getReal() == 0) {
3759 return ddd;
3760 }
3761 if (dd.getReal() == 0) {
3762 return dm;
3763 }
3764 }
3765
3766 final double a = zhm / (xmm - xm);
3767 final T ylog = dm.divide(dd).log().multiply(a);
3768 if (ylog.getReal() < -10.) {
3769 return dd;
3770 } else if (ylog.getReal() > 10.) {
3771 return dm;
3772 } else {
3773 return ylog.exp().add(1).pow(1.0 / a).multiply(dd);
3774 }
3775 }
3776
3777
3778
3779
3780
3781
3782
3783
3784
3785 private T splini(final T[] xa, final T[] ya, final T[] y2a, final T x) {
3786 final int n = xa.length;
3787 T yi = zero;
3788 int klo = 0;
3789 int khi = 1;
3790 while (x.getReal() > xa[klo].getReal() && khi < n) {
3791 T xx = x;
3792 if (khi < n - 1) {
3793 xx = (x.getReal() < xa[khi].getReal()) ? x : xa[khi];
3794 }
3795 final T h = xa[khi].subtract(xa[klo]);
3796 final T a = xa[khi].subtract(xx).divide(h);
3797 final T b = xx.subtract(xa[klo]).divide(h);
3798 final T a2 = a.square();
3799 final T b2 = b.square();
3800
3801 final T z =
3802 a2.divide(2).subtract(a2.square().add(1).divide(4)).multiply(y2a[klo]).
3803 add(b2.multiply(b2).divide(4).subtract(b2.divide(2)).multiply(y2a[khi]));
3804 yi = yi.add( a2.negate().add(1).multiply(ya[klo]).divide(2).
3805 add(b2.multiply(ya[khi]).divide(2)).
3806 add(z.multiply(h).multiply(h).divide(6)).
3807 multiply(h));
3808 klo++;
3809 khi++;
3810 }
3811 return yi;
3812 }
3813
3814
3815
3816
3817
3818
3819
3820
3821
3822 private T splint(final T[] xa, final T[] ya, final T[] y2a, final T x) {
3823 final int n = xa.length;
3824 int klo = 0;
3825 int khi = n - 1;
3826 while (khi - klo > 1) {
3827 final int k = (khi + klo) >>> 1;
3828 if (xa[k].getReal() > x.getReal()) {
3829 khi = k;
3830 } else {
3831 klo = k;
3832 }
3833 }
3834 final T h = xa[khi].subtract(xa[klo]);
3835 final T a = xa[khi].subtract(x).divide(h);
3836 final T b = x.subtract(xa[klo]).divide(h);
3837 return a.multiply(ya[klo]).add(b.multiply(ya[khi])).
3838 add(( a.square().multiply(a).subtract(a).multiply(y2a[klo]).
3839 add(b.multiply(b).multiply(b).subtract(b).multiply(y2a[khi]))
3840 ).multiply(h).multiply(h).divide(6));
3841 }
3842
3843
3844
3845
3846
3847
3848
3849
3850
3851 private T[] spline(final T[] x, final T[] y, final T yp1, final T ypn) {
3852 final int n = x.length;
3853 final T[] y2 = MathArrays.buildArray(field, n);
3854 final T[] u = MathArrays.buildArray(field, n);
3855
3856 if (yp1.getReal() < 1e+30) {
3857 y2[0] = zero.newInstance(-0.5);
3858 final T dx = x[1].subtract(x[0]);
3859 final T dy = y[1].subtract(y[0]);
3860 u[0] = dx.reciprocal().multiply(3.0).multiply(dy.divide(dx).subtract(yp1));
3861 }
3862 for (int i = 1; i < n - 1; i++) {
3863 final T dx0m = x[i].subtract(x[i - 1]);
3864 final T dy0m = y[i].subtract(y[i - 1]);
3865 final T dxpm = x[i + 1].subtract(x[i - 1]);
3866 final T dxp0 = x[i + 1].subtract(x[i]);
3867 final T dyp0 = y[i + 1].subtract(y[i]);
3868 final T sig = dx0m.divide(dxpm);
3869 final T p = sig.multiply(y2[i - 1]).add(2.0);
3870 y2[i] = sig.subtract(1.0).divide(p);
3871 u[i] = dyp0.divide(dxp0).subtract(dy0m.divide(dx0m)).multiply(6).divide(dxpm).subtract(sig.multiply(u[i - 1])).divide(p);
3872 }
3873
3874 double qn = 0;
3875 T un = zero;
3876 if (ypn.getReal() < 1e+30) {
3877 final T dx12 = x[n - 1].subtract(x[n - 2]);
3878 final T dy12 = y[n - 1].subtract(y[n - 2]);
3879 qn = 0.5;
3880 un = dx12.reciprocal().multiply(3.0).multiply(ypn.subtract(dy12.divide(dx12)));
3881 }
3882
3883 y2[n - 1] = un.subtract(u[n - 2].multiply(qn)).divide(y2[n - 2].multiply(qn).add(1.0));
3884 for (int k = n - 2; k >= 0; k--) {
3885 y2[k] = y2[k].multiply(y2[k + 1]).add(u[k]);
3886 }
3887
3888 return y2;
3889
3890 }
3891
3892
3893
3894
3895
3896
3897
3898 private T densm(final T alt, final T d0, final double xm) {
3899
3900 T densm = d0;
3901
3902
3903 int mn = ZN2.length;
3904 T z = (alt.getReal() > ZN2[mn - 1]) ? alt : zero.newInstance(ZN2[mn - 1]);
3905
3906 double z1 = ZN2[0];
3907 double z2 = ZN2[mn - 1];
3908 T t1 = meso_tn2[0];
3909 T t2 = meso_tn2[mn - 1];
3910 T zg = zeta(z, z1);
3911 T zgdif = zeta(zero.newInstance(z2), z1);
3912
3913
3914 T[] xs = MathArrays.buildArray(field, mn);
3915 T[] ys = MathArrays.buildArray(field, mn);
3916 for (int k = 0; k < mn; k++) {
3917 xs[k] = zeta(zero.newInstance(ZN2[k]), z1).divide(zgdif);
3918 ys[k] = meso_tn2[k].reciprocal();
3919 }
3920 final T qSM = rlat.add(z2).divide(rlat.add(z1));
3921 T yd1 = meso_tgn2[0].negate().divide(t1.square()).multiply(zgdif);
3922 T yd2 = meso_tgn2[1].negate().divide(t2.square()).multiply(zgdif).multiply(qSM.square());
3923
3924
3925 T[] y2out = spline(xs, ys, yd1, yd2);
3926 T x = zg.divide(zgdif);
3927 T y = splint(xs, ys, y2out, x);
3928
3929
3930 T tz = y.reciprocal();
3931
3932 if (xm != 0.0) {
3933
3934 final T glb = galt(zero.newInstance(z1));
3935 final T gamm = glb.multiply(zgdif).multiply(xm / R_GAS);
3936
3937
3938 final T yi = splini(xs, ys, y2out, x);
3939 final T expl = min(MIN_TEMP, gamm.multiply(yi));
3940
3941
3942 densm = densm.multiply(t1.divide(tz).multiply(expl.negate().exp()));
3943 }
3944
3945 if (alt.getReal() > ZN3[0]) {
3946 return (xm == 0.0) ? tz : densm;
3947 }
3948
3949
3950 z = alt;
3951 mn = ZN3.length;
3952 z1 = ZN3[0];
3953 z2 = ZN3[mn - 1];
3954 t1 = meso_tn3[0];
3955 t2 = meso_tn3[mn - 1];
3956 zg = zeta(z, z1);
3957 zgdif = zeta(zero.newInstance(z2), z1);
3958
3959
3960 xs = MathArrays.buildArray(field, mn);
3961 ys = MathArrays.buildArray(field, mn);
3962 for (int k = 0; k < mn; k++) {
3963 xs[k] = zeta(zero.newInstance(ZN3[k]), z1).divide(zgdif);
3964 ys[k] = meso_tn3[k].reciprocal();
3965 }
3966 final T qTS = rlat.add(z2) .divide(rlat.add(z1));
3967 yd1 = meso_tgn3[0].negate().divide(t1.multiply(t1)).multiply(zgdif);
3968 yd2 = meso_tgn3[1].negate().divide(t2.multiply(t2)).multiply(zgdif).multiply(qTS).multiply(qTS);
3969
3970
3971 y2out = spline(xs, ys, yd1, yd2);
3972 x = zg.divide(zgdif);
3973 y = splint(xs, ys, y2out, x);
3974
3975
3976 tz = y.reciprocal();
3977
3978 if (xm != 0.0) {
3979
3980 final T glb = galt(zero.newInstance(z1));
3981 final T gamm = glb.multiply(zgdif).multiply(xm / R_GAS);
3982
3983
3984 final T yi = splini(xs, ys, y2out, x);
3985 final T expl = min(MIN_TEMP, gamm.multiply(yi));
3986
3987
3988 densm = densm.multiply(t1.divide(tz).multiply(expl.negate().exp()));
3989 }
3990
3991 return (xm == 0.0) ? tz : densm;
3992 }
3993
3994
3995
3996
3997
3998
3999
4000
4001
4002
4003
4004
4005 private T densu(final T alt, final T dlb, final T tinf,
4006 final T tlb, final double xm, final double alpha,
4007 final double zlb, final T s2) {
4008
4009 T z = (alt.getReal() > ZN1[0]) ? alt : zero.newInstance(ZN1[0]);
4010
4011
4012 final T zg2 = zeta(z, zlb);
4013
4014
4015 final T tt = tinf.subtract(tinf.subtract(tlb).multiply(s2.negate().multiply(zg2).exp()));
4016 final T ta = tt;
4017 T tz = tt;
4018
4019 final int mn = ZN1.length;
4020 final T[] xs = MathArrays.buildArray(field, mn);
4021 final T[] ys = MathArrays.buildArray(field, mn);
4022 T x = zero;
4023 T[] y2out = MathArrays.buildArray(field, mn);
4024 T zgdif = zero;
4025 if (alt.getReal() < ZN1[0]) {
4026
4027
4028 final T p = rlat.add(zlb).divide(rlat.add(ZN1[0]));
4029 final T dta = tinf.subtract(ta).multiply(s2).multiply(p.square());
4030 meso_tgn1[0] = dta;
4031 meso_tn1[0] = ta;
4032 final T tzn1mn1 = zero.newInstance(ZN1[mn - 1]);
4033 z = (alt.getReal() > ZN1[mn - 1]) ? alt : tzn1mn1;
4034
4035 final T t1 = meso_tn1[0];
4036 final T t2 = meso_tn1[mn - 1];
4037
4038 final T zg = zeta(z, ZN1[0]);
4039 zgdif = zeta(tzn1mn1, ZN1[0]);
4040
4041 for (int k = 0; k < mn; k++) {
4042 xs[k] = zeta(zero.newInstance(ZN1[k]), ZN1[0]).divide(zgdif);
4043 ys[k] = meso_tn1[k].reciprocal();
4044 }
4045
4046 final T q = rlat.add(ZN1[mn - 1]).divide(rlat.add(ZN1[0]));
4047 final T yd1 = meso_tgn1[0].negate().divide(t1.square()).multiply(zgdif);
4048 final T yd2 = meso_tgn1[1].negate().divide(t2.square()).multiply(zgdif).multiply(q.square());
4049
4050 y2out = spline(xs, ys, yd1, yd2);
4051 x = zg.divide(zgdif);
4052 final T y = splint(xs, ys, y2out, x);
4053
4054 tz = y.reciprocal();
4055 }
4056
4057 if (xm == 0) {
4058 return tz;
4059 }
4060
4061
4062 T glb = galt(zero.newInstance(zlb));
4063 T gamma = glb.divide(s2.multiply(tinf)).multiply(xm / R_GAS);
4064 T expl = tt.getReal() <= 0 ?
4065 zero.newInstance(MIN_TEMP) :
4066 min(MIN_TEMP, s2.negate().multiply(gamma).multiply(zg2).exp());
4067 T densu = dlb.multiply(expl).multiply(tlb.divide(tt).pow(gamma.add(alpha + 1)));
4068
4069
4070 if (!Double.isFinite(densu.getReal())) {
4071 if (expl.getReal() < MIN_TEMP) {
4072 densu = dlb.multiply(FastMath.exp((FastMath.log(tlb.divide(tt)).multiply(gamma.add(alpha + 1))).
4073 subtract(s2.multiply(gamma).multiply(zg2))));
4074 } else {
4075 throw new OrekitException(OrekitMessages.INFINITE_NRLMSISE00_DENSITY);
4076 }
4077 }
4078
4079
4080 if (alt.getReal() < ZN1[0]) {
4081 glb = galt(zero.newInstance(ZN1[0]));
4082 gamma = glb.multiply(zgdif).multiply(xm / R_GAS);
4083
4084 expl = tz.getReal() <= 0 ?
4085 zero.newInstance(MIN_TEMP) :
4086 min(MIN_TEMP, gamma.multiply(splini(xs, ys, y2out, x)));
4087
4088 densu = densu.multiply(meso_tn1[0].divide(tz).pow(alpha + 1).multiply(expl.negate().exp()));
4089 }
4090
4091
4092 return densu;
4093 }
4094
4095
4096
4097
4098
4099
4100 private T min(final double d, final T f) {
4101 return (f.getReal() > d) ? zero.newInstance(d) : f;
4102 }
4103
4104
4105
4106
4107
4108 private T galt(final T alt) {
4109 final T r = alt.divide(rlat).add(1);
4110 return glat.divide(r.square());
4111 }
4112
4113
4114
4115
4116
4117
4118 private T zeta(final T zz, final double zl) {
4119 return zz.subtract(zl).multiply(rlat.add(zl)).divide(rlat.add(zz));
4120 }
4121
4122 }
4123
4124 }