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