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