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