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