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