1 /* Copyright 2002-2024 CS GROUP
2 * Licensed to CS GROUP (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.models.earth.ionosphere;
18
19 import org.hipparchus.util.FastMath;
20 import org.hipparchus.util.MathArrays;
21 import org.hipparchus.util.SinCos;
22 import org.orekit.time.DateComponents;
23 import org.orekit.time.DateTimeComponents;
24 import org.orekit.time.TimeComponents;
25
26 /**
27 * This class perfoms the computation of the parameters used by the NeQuick model.
28 *
29 * @author Bryan Cazabonne
30 *
31 * @see "European Union (2016). European GNSS (Galileo) Open Service-Ionospheric Correction
32 * Algorithm for Galileo Single Frequency Users. 1.2."
33 *
34 * @since 10.1
35 */
36 class NeQuickParameters {
37
38 /** Solar zenith angle at day night transition, degrees. */
39 private static final double X0 = 86.23292796211615;
40
41 /** F2 layer maximum density. */
42 private final double nmF2;
43
44 /** F2 layer maximum density height [km]. */
45 private final double hmF2;
46
47 /** F1 layer maximum density height [km]. */
48 private final double hmF1;
49
50 /** E layer maximum density height [km]. */
51 private final double hmE;
52
53 /** F2 layer bottom thickness parameter [km]. */
54 private final double b2Bot;
55
56 /** F1 layer top thickness parameter [km]. */
57 private final double b1Top;
58
59 /** F1 layer bottom thickness parameter [km]. */
60 private final double b1Bot;
61
62 /** E layer top thickness parameter [km]. */
63 private final double beTop;
64
65 /** E layer bottom thickness parameter [km]. */
66 private final double beBot;
67
68 /** topside thickness parameter [km]. */
69 private final double h0;
70
71 /** Layer amplitudes. */
72 private final double[] amplitudes;
73
74 /**
75 * Build a new instance.
76 * @param dateTime current date time components
77 * @param f2 F2 coefficients used by the F2 layer
78 * @param fm3 Fm3 coefficients used by the F2 layer
79 * @param latitude latitude of a point along the integration path, in radians
80 * @param longitude longitude of a point along the integration path, in radians
81 * @param alpha effective ionisation level coefficients
82 * @param modipGrip modip grid
83 */
84 NeQuickParameters(final DateTimeComponents dateTime, final double[][][] f2,
85 final double[][][] fm3, final double latitude, final double longitude,
86 final double[] alpha, final double[][] modipGrip) {
87
88 // MODIP in degrees
89 final double modip = computeMODIP(latitude, longitude, modipGrip);
90 // Effective ionisation level Az
91 final double az = computeAz(modip, alpha);
92 // Effective sunspot number (Eq. 19)
93 final double azr = FastMath.sqrt(167273.0 + (az - 63.7) * 1123.6) - 408.99;
94 // Date and Time components
95 final DateComponents date = dateTime.getDate();
96 final TimeComponents time = dateTime.getTime();
97 // Hours
98 final double hours = time.getSecondsInUTCDay() / 3600.0;
99 // Effective solar zenith angle in radians
100 final double xeff = computeEffectiveSolarAngle(date.getMonth(), hours, latitude, longitude);
101
102 // Coefficients for F2 layer parameters
103 // Compute the array of interpolated coefficients for foF2 (Eq. 44)
104 final double[][] af2 = new double[76][13];
105 for (int j = 0; j < 76; j++) {
106 for (int k = 0; k < 13; k++ ) {
107 af2[j][k] = f2[0][j][k] * (1.0 - (azr * 0.01)) + f2[1][j][k] * (azr * 0.01);
108 }
109 }
110
111 // Compute the array of interpolated coefficients for M(3000)F2 (Eq. 46)
112 final double[][] am3 = new double[49][9];
113 for (int j = 0; j < 49; j++) {
114 for (int k = 0; k < 9; k++ ) {
115 am3[j][k] = fm3[0][j][k] * (1.0 - (azr * 0.01)) + fm3[1][j][k] * (azr * 0.01);
116 }
117 }
118
119 // E layer maximum density height in km (Eq. 78)
120 this.hmE = 120.0;
121 // E layer critical frequency in MHz
122 final double foE = computefoE(date.getMonth(), az, xeff, latitude);
123 // E layer maximum density in 10^11 m-3 (Eq. 36)
124 final double nmE = 0.124 * foE * foE;
125
126 // Time argument (Eq. 49)
127 final double t = FastMath.toRadians(15 * hours) - FastMath.PI;
128 // Compute Fourier time series for foF2 and M(3000)F2
129 final double[] cf2 = computeCF2(af2, t);
130 final double[] cm3 = computeCm3(am3, t);
131 // F2 layer critical frequency in MHz
132 final double foF2 = computefoF2(modip, cf2, latitude, longitude);
133 // Maximum Usable Frequency factor
134 final double mF2 = computeMF2(modip, cm3, latitude, longitude);
135 // F2 layer maximum density in 10^11 m-3
136 this.nmF2 = 0.124 * foF2 * foF2;
137 // F2 layer maximum density height in km
138 this.hmF2 = computehmF2(foE, foF2, mF2);
139
140 // F1 layer critical frequency in MHz
141 final double foF1 = computefoF1(foE, foF2);
142 // F1 layer maximum density in 10^11 m-3
143 final double nmF1;
144 if (foF1 <= 0.0 && foE > 2.0) {
145 final double foEpopf = foE + 0.5;
146 nmF1 = 0.124 * foEpopf * foEpopf;
147 } else {
148 nmF1 = 0.124 * foF1 * foF1;
149 }
150 // F1 layer maximum density height in km
151 this.hmF1 = 0.5 * (hmF2 + hmE);
152
153 // Thickness parameters (Eq. 85 to 89)
154 final double a = 0.01 * clipExp(-3.467 + 0.857 * FastMath.log(foF2 * foF2) + 2.02 * FastMath.log(mF2));
155 this.b2Bot = 0.385 * nmF2 / a;
156 this.b1Top = 0.3 * (hmF2 - hmF1);
157 this.b1Bot = 0.5 * (hmF1 - hmE);
158 this.beTop = FastMath.max(b1Bot, 7.0);
159 this.beBot = 5.0;
160
161 // Layer amplitude coefficients
162 this.amplitudes = computeLayerAmplitudes(nmE, nmF1, foF1);
163
164 // Topside thickness parameter
165 this.h0 = computeH0(date.getMonth(), azr);
166 }
167
168 /**
169 * Get the F2 layer maximum density.
170 * @return nmF2
171 */
172 public double getNmF2() {
173 return nmF2;
174 }
175
176 /**
177 * Get the F2 layer maximum density height.
178 * @return hmF2 in km
179 */
180 public double getHmF2() {
181 return hmF2;
182 }
183
184 /**
185 * Get the F1 layer maximum density height.
186 * @return hmF1 in km
187 */
188 public double getHmF1() {
189 return hmF1;
190 }
191
192 /**
193 * Get the E layer maximum density height.
194 * @return hmE in km
195 */
196 public double getHmE() {
197 return hmE;
198 }
199
200 /**
201 * Get the F2 layer thickness parameter (bottom).
202 * @return B2Bot in km
203 */
204 public double getB2Bot() {
205 return b2Bot;
206 }
207
208 /**
209 * Get the F1 layer thickness parameter (top).
210 * @return B1Top in km
211 */
212 public double getB1Top() {
213 return b1Top;
214 }
215
216 /**
217 * Get the F1 layer thickness parameter (bottom).
218 * @return B1Bot in km
219 */
220 public double getB1Bot() {
221 return b1Bot;
222 }
223
224 /**
225 * Get the E layer thickness parameter (bottom).
226 * @return BeBot in km
227 */
228 public double getBEBot() {
229 return beBot;
230 }
231
232 /**
233 * Get the E layer thickness parameter (top).
234 * @return BeTop in km
235 */
236 public double getBETop() {
237 return beTop;
238 }
239
240 /**
241 * Get the F2, F1 and E layer amplitudes.
242 * <p>
243 * The resulting element is an array having the following form:
244 * <ul>
245 * <li>double[0] = A1 → F2 layer amplitude
246 * <li>double[1] = A2 → F1 layer amplitude
247 * <li>double[2] = A3 → E layer amplitude
248 * </ul>
249 * @return layer amplitudes
250 */
251 public double[] getLayerAmplitudes() {
252 return amplitudes.clone();
253 }
254
255 /**
256 * Get the topside thickness parameter H0.
257 * @return H0 in km
258 */
259 public double getH0() {
260 return h0;
261 }
262
263 /**
264 * Computes the value of the modified dip latitude (MODIP) for the
265 * given latitude and longitude.
266 *
267 * @param lat receiver latitude, radians
268 * @param lon receiver longitude, radians
269 * @param stModip modip grid
270 * @return the MODIP in degrees
271 */
272 private double computeMODIP(final double lat, final double lon, final double[][] stModip) {
273
274 // For the MODIP computation, latitude and longitude have to be converted in degrees
275 final double latitude = FastMath.toDegrees(lat);
276 final double longitude = FastMath.toDegrees(lon);
277
278 // Extreme cases
279 if (latitude == 90.0 || latitude == -90.0) {
280 return latitude;
281 }
282
283 // Auxiliary parameter l (Eq. 6 to 8)
284 final int lF = (int) ((longitude + 180) * 0.1);
285 int l = lF - 2;
286 if (l < -2) {
287 l += 36;
288 } else if (l > 33) {
289 l -= 36;
290 }
291
292 // Auxiliary parameter a (Eq. 9 to 11)
293 final double a = 0.2 * (latitude + 90) + 1.0;
294 final double aF = FastMath.floor(a);
295 // Eq. 10
296 final double x = a - aF;
297 // Eq. 11
298 final int i = (int) aF - 2;
299
300 // zi coefficients (Eq. 12 and 13)
301 final double z1 = interpolate(stModip[i + 1][l + 2], stModip[i + 2][l + 2], stModip[i + 3][l + 2], stModip[i + 4][l + 2], x);
302 final double z2 = interpolate(stModip[i + 1][l + 3], stModip[i + 2][l + 3], stModip[i + 3][l + 3], stModip[i + 4][l + 3], x);
303 final double z3 = interpolate(stModip[i + 1][l + 4], stModip[i + 2][l + 4], stModip[i + 3][l + 4], stModip[i + 4][l + 4], x);
304 final double z4 = interpolate(stModip[i + 1][l + 5], stModip[i + 2][l + 5], stModip[i + 3][l + 5], stModip[i + 4][l + 5], x);
305
306 // Auxiliary parameter b (Eq. 14 and 15)
307 final double b = (longitude + 180) * 0.1;
308 final double bF = FastMath.floor(b);
309 final double y = b - bF;
310
311 // MODIP (Ref Eq. 16)
312 final double modip = interpolate(z1, z2, z3, z4, y);
313
314 return modip;
315 }
316
317 /**
318 * This method computes the effective ionisation level Az.
319 * <p>
320 * This parameter is used for the computation of the Total Electron Content (TEC).
321 * </p>
322 * @param modip modified dip latitude (MODIP) in degrees
323 * @param alpha effective ionisation level coefficients
324 * @return the ionisation level Az
325 */
326 private double computeAz(final double modip, final double[] alpha) {
327 // Particular condition (Eq. 17)
328 if (alpha[0] == 0.0 && alpha[1] == 0.0 && alpha[2] == 0.0) {
329 return 63.7;
330 }
331 // Az = a0 + modip * a1 + modip^2 * a2 (Eq. 18)
332 double az = alpha[0] + modip * (alpha[1] + modip * alpha[2]);
333 // If Az < 0 -> Az = 0
334 az = FastMath.max(0.0, az);
335 // If Az > 400 -> Az = 400
336 az = FastMath.min(400.0, az);
337 return az;
338 }
339
340 /**
341 * This method computes the effective solar zenith angle.
342 * <p>
343 * The effective solar zenith angle is compute as a function of the
344 * solar zenith angle and the solar zenith angle at day night transition.
345 * </p>
346 * @param month current month of the year
347 * @param hours universal time (hours)
348 * @param latitude in radians
349 * @param longitude in radians
350 * @return the effective solar zenith angle, radians
351 */
352 private double computeEffectiveSolarAngle(final int month,
353 final double hours,
354 final double latitude,
355 final double longitude) {
356 // Local time (Eq.4)
357 final double lt = hours + longitude / FastMath.toRadians(15.0);
358 // Day of year at the middle of the month (Eq. 20)
359 final double dy = 30.5 * month - 15.0;
360 // Time (Eq. 21)
361 final double t = dy + (18 - hours) / 24;
362 // Arguments am and al (Eq. 22 and 23)
363 final double am = FastMath.toRadians(0.9856 * t - 3.289);
364 final double al = am + FastMath.toRadians(1.916 * FastMath.sin(am) + 0.020 * FastMath.sin(2.0 * am) + 282.634);
365 // Sine and cosine of solar declination (Eq. 24 and 25)
366 final double sDec = 0.39782 * FastMath.sin(al);
367 final double cDec = FastMath.sqrt(1. - sDec * sDec);
368 // Solar zenith angle, deg (Eq. 26 and 27)
369 final SinCos scLat = FastMath.sinCos(latitude);
370 final double coef = (FastMath.PI / 12) * (12 - lt);
371 final double cZenith = scLat.sin() * sDec + scLat.cos() * cDec * FastMath.cos(coef);
372 final double angle = FastMath.atan2(FastMath.sqrt(1.0 - cZenith * cZenith), cZenith);
373 final double x = FastMath.toDegrees(angle);
374 // Effective solar zenith angle (Eq. 28)
375 final double xeff = join(90.0 - 0.24 * clipExp(20.0 - 0.2 * x), x, 12.0, x - X0);
376 return FastMath.toRadians(xeff);
377 }
378
379 /**
380 * This method computes the E layer critical frequency at a given location.
381 * @param month current month
382 * @param az ffective ionisation level
383 * @param xeff effective solar zenith angle in radians
384 * @param latitude latitude in radians
385 * @return the E layer critical frequency at a given location in MHz
386 */
387 private double computefoE(final int month, final double az,
388 final double xeff, final double latitude) {
389 // The latitude has to be converted in degrees
390 final double lat = FastMath.toDegrees(latitude);
391 // Square root of the effective ionisation level
392 final double sqAz = FastMath.sqrt(az);
393 // seas parameter (Eq. 30 to 32)
394 final int seas;
395 if (month == 1 || month == 2 || month == 11 || month == 12) {
396 seas = -1;
397 } else if (month == 3 || month == 4 || month == 9 || month == 10) {
398 seas = 0;
399 } else {
400 seas = 1;
401 }
402 // Latitudinal dependence (Eq. 33 and 34)
403 final double ee = clipExp(0.3 * lat);
404 final double seasp = seas * ((ee - 1.0) / (ee + 1.0));
405 // Critical frequency (Eq. 35)
406 final double coef = 1.112 - 0.019 * seasp;
407 final double foE = FastMath.sqrt(coef * coef * sqAz * FastMath.pow(FastMath.cos(xeff), 0.6) + 0.49);
408 return foE;
409 }
410
411 /**
412 * Computes the F2 layer height of maximum electron density.
413 * @param foE E layer layer critical frequency in MHz
414 * @param foF2 F2 layer layer critical frequency in MHz
415 * @param mF2 maximum usable frequency factor
416 * @return hmF2 in km
417 */
418 private double computehmF2(final double foE, final double foF2, final double mF2) {
419 // Ratio
420 final double fo = foF2 / foE;
421 final double ratio = join(fo, 1.75, 20.0, fo - 1.75);
422
423 // deltaM parameter
424 double deltaM = -0.012;
425 if (foE >= 1e-30) {
426 deltaM += 0.253 / (ratio - 1.215);
427 }
428
429 // hmF2 Eq. 80
430 final double mF2Sq = mF2 * mF2;
431 final double temp = FastMath.sqrt((0.0196 * mF2Sq + 1) / (1.2967 * mF2Sq - 1.0));
432 final double height = ((1490.0 * mF2 * temp) / (mF2 + deltaM)) - 176.0;
433 return height;
434 }
435
436 /**
437 * Computes cf2 coefficients.
438 * @param af2 interpolated coefficients for foF2
439 * @param t time argument
440 * @return the cf2 coefficients array
441 */
442 private double[] computeCF2(final double[][] af2, final double t) {
443 // Eq. 50
444 final double[] cf2 = new double[76];
445 for (int i = 0; i < cf2.length; i++) {
446 double sum = 0.0;
447 for (int k = 0; k < 6; k++) {
448 final SinCos sc = FastMath.sinCos((k + 1) * t);
449 sum += af2[i][2 * k + 1] * sc.sin() + af2[i][2 * (k + 1)] * sc.cos();
450 }
451 cf2[i] = af2[i][0] + sum;
452 }
453 return cf2;
454 }
455
456 /**
457 * Computes Cm3 coefficients.
458 * @param am3 interpolated coefficients for foF2
459 * @param t time argument
460 * @return the Cm3 coefficients array
461 */
462 private double[] computeCm3(final double[][] am3, final double t) {
463 // Eq. 51
464 final double[] cm3 = new double[49];
465 for (int i = 0; i < cm3.length; i++) {
466 double sum = 0.0;
467 for (int k = 0; k < 4; k++) {
468 final SinCos sc = FastMath.sinCos((k + 1) * t);
469 sum += am3[i][2 * k + 1] * sc.sin() + am3[i][2 * (k + 1)] * sc.cos();
470 }
471 cm3[i] = am3[i][0] + sum;
472 }
473 return cm3;
474 }
475
476 /**
477 * This method computes the F2 layer critical frequency.
478 * @param modip modified DIP latitude, in degrees
479 * @param cf2 Fourier time series for foF2
480 * @param latitude latitude in radians
481 * @param longitude longitude in radians
482 * @return the F2 layer critical frequency, MHz
483 */
484 private double computefoF2(final double modip, final double[] cf2,
485 final double latitude, final double longitude) {
486
487 // Legendre grades (Eq. 63)
488 final int[] q = new int[] {
489 12, 12, 9, 5, 2, 1, 1, 1, 1
490 };
491
492 // Array for geographic terms
493 final double[] g = new double[cf2.length];
494 g[0] = 1.0;
495
496 // MODIP coefficients Eq. 57
497 final double sinMODIP = FastMath.sin(FastMath.toRadians(modip));
498 final double[] m = new double[12];
499 m[0] = 1.0;
500 for (int i = 1; i < q[0]; i++) {
501 m[i] = sinMODIP * m[i - 1];
502 g[i] = m[i];
503 }
504
505 // Latitude coefficients (Eq. 58)
506 final double cosLat = FastMath.cos(latitude);
507 final double[] p = new double[8];
508 p[0] = cosLat;
509 for (int n = 2; n < 9; n++) {
510 p[n - 1] = cosLat * p[n - 2];
511 }
512
513 // latitude and longitude terms
514 int index = 12;
515 for (int i = 1; i < q.length; i++) {
516 final SinCos sc = FastMath.sinCos(i * longitude);
517 for (int j = 0; j < q[i]; j++) {
518 g[index++] = m[j] * p[i - 1] * sc.cos();
519 g[index++] = m[j] * p[i - 1] * sc.sin();
520 }
521 }
522
523 // Compute foF2 by linear combination
524 final double frequency = MathArrays.linearCombination(cf2, g);
525 return frequency;
526 }
527
528 /**
529 * This method computes the Maximum Usable Frequency factor.
530 * @param modip modified DIP latitude, in degrees
531 * @param cm3 Fourier time series for M(3000)F2
532 * @param latitude latitude in radians
533 * @param longitude longitude in radians
534 * @return the Maximum Usable Frequency factor
535 */
536 private double computeMF2(final double modip, final double[] cm3,
537 final double latitude, final double longitude) {
538
539 // Legendre grades (Eq. 71)
540 final int[] r = new int[] {
541 7, 8, 6, 3, 2, 1, 1
542 };
543
544 // Array for geographic terms
545 final double[] g = new double[cm3.length];
546 g[0] = 1.0;
547
548 // MODIP coefficients Eq. 57
549 final double sinMODIP = FastMath.sin(FastMath.toRadians(modip));
550 final double[] m = new double[12];
551 m[0] = 1.0;
552 for (int i = 1; i < 12; i++) {
553 m[i] = sinMODIP * m[i - 1];
554 if (i < 7) {
555 g[i] = m[i];
556 }
557 }
558
559 // Latitude coefficients (Eq. 58)
560 final double cosLat = FastMath.cos(latitude);
561 final double[] p = new double[8];
562 p[0] = cosLat;
563 for (int n = 2; n < 9; n++) {
564 p[n - 1] = cosLat * p[n - 2];
565 }
566
567 // latitude and longitude terms
568 int index = 7;
569 for (int i = 1; i < r.length; i++) {
570 final SinCos sc = FastMath.sinCos(i * longitude);
571 for (int j = 0; j < r[i]; j++) {
572 g[index++] = m[j] * p[i - 1] * sc.cos();
573 g[index++] = m[j] * p[i - 1] * sc.sin();
574 }
575 }
576
577 // Compute m3000 by linear combination
578 final double m3000 = MathArrays.linearCombination(g, cm3);
579 return m3000;
580 }
581
582 /**
583 * This method computes the F1 layer critical frequency.
584 * <p>
585 * This computation performs the algorithm exposed in Annex F
586 * of the reference document.
587 * </p>
588 * @param foE the E layer critical frequency, MHz
589 * @return the F1 layer critical frequency, MHz
590 * @param foF2 the F2 layer critical frequency, MHz
591 */
592 private double computefoF1(final double foE, final double foF2) {
593 final double temp = join(1.4 * foE, 0.0, 1000.0, foE - 2.0);
594 final double temp2 = join(0.0, temp, 1000.0, foE - temp);
595 final double value = join(temp2, 0.85 * temp2, 60.0, 0.85 * foF2 - temp2);
596 if (value < 1.0E-6) {
597 return 0.0;
598 } else {
599 return value;
600 }
601 }
602
603 /**
604 * This method allows the computation of the F2, F1 and E layer amplitudes.
605 * <p>
606 * The resulting element is an array having the following form:
607 * <ul>
608 * <li>double[0] = A1 → F2 layer amplitude
609 * <li>double[1] = A2 → F1 layer amplitude
610 * <li>double[2] = A3 → E layer amplitude
611 * </ul>
612 * </p>
613 * @param nmE E layer maximum density in 10^11 m-3
614 * @param nmF1 F1 layer maximum density in 10^11 m-3
615 * @param foF1 F1 layer critical frequency in MHz
616 * @return a three components array containing the layer amplitudes
617 */
618 private double[] computeLayerAmplitudes(final double nmE, final double nmF1, final double foF1) {
619 // Initialize array
620 final double[] amplitude = new double[3];
621
622 // F2 layer amplitude (Eq. 90)
623 final double a1 = 4.0 * nmF2;
624 amplitude[0] = a1;
625
626 // F1 and E layer amplitudes (Eq. 91 to 98)
627 if (foF1 < 0.5) {
628 amplitude[1] = 0.0;
629 amplitude[2] = 4.0 * (nmE - epst(a1, hmF2, b2Bot, hmE));
630 } else {
631 double a2a = 0.0;
632 double a3a = 4.0 * nmE;
633 for (int i = 0; i < 5; i++) {
634 a2a = 4.0 * (nmF1 - epst(a1, hmF2, b2Bot, hmF1) - epst(a3a, hmE, beTop, hmF1));
635 a2a = join(a2a, 0.8 * nmF1, 1.0, a2a - 0.8 * nmF1);
636 a3a = 4.0 * (nmE - epst(a2a, hmF1, b1Bot, hmE) - epst(a1, hmF2, b2Bot, hmE));
637 }
638 amplitude[1] = a2a;
639 amplitude[2] = join(a3a, 0.05, 60.0, a3a - 0.005);
640 }
641
642 return amplitude;
643 }
644
645 /**
646 * This method computes the topside thickness parameter H0.
647 *
648 * @param month current month
649 * @param azr effective sunspot number
650 * @return H0 in km
651 */
652 private double computeH0(final int month, final double azr) {
653
654 // Auxiliary parameter ka (Eq. 99 and 100)
655 final double ka;
656 if (month > 3 && month < 10) {
657 // month = 4,5,6,7,8,9
658 ka = 6.705 - 0.014 * azr - 0.008 * hmF2;
659 } else {
660 // month = 1,2,3,10,11,12
661 final double ratio = hmF2 / b2Bot;
662 ka = -7.77 + 0.097 * ratio * ratio + 0.153 * nmF2;
663 }
664
665 // Auxiliary parameter kb (Eq. 101 and 102)
666 double kb = join(ka, 2.0, 1.0, ka - 2.0);
667 kb = join(8.0, kb, 1.0, kb - 8.0);
668
669 // Auxiliary parameter Ha (Eq. 103)
670 final double hA = kb * b2Bot;
671
672 // Auxiliary parameters x and v (Eq. 104 and 105)
673 final double x = 0.01 * (hA - 150.0);
674 final double v = (0.041163 * x - 0.183981) * x + 1.424472;
675
676 // Topside thickness parameter (Eq. 106)
677 final double h = hA / v;
678 return h;
679 }
680
681 /**
682 * A clipped exponential function.
683 * <p>
684 * This function, describe in section F.2.12.2 of the reference document, is
685 * recommanded for the computation of exponential values.
686 * </p>
687 * @param power power for exponential function
688 * @return clipped exponential value
689 */
690 private double clipExp(final double power) {
691 if (power > 80.0) {
692 return 5.5406E34;
693 } else if (power < -80) {
694 return 1.8049E-35;
695 } else {
696 return FastMath.exp(power);
697 }
698 }
699
700 /**
701 * This method provides a third order interpolation function
702 * as recommended in the reference document (Ref Eq. 128 to Eq. 138)
703 *
704 * @param z1 z1 coefficient
705 * @param z2 z2 coefficient
706 * @param z3 z3 coefficient
707 * @param z4 z4 coefficient
708 * @param x position
709 * @return a third order interpolation
710 */
711 private double interpolate(final double z1, final double z2,
712 final double z3, final double z4,
713 final double x) {
714
715 if (FastMath.abs(2.0 * x) < 1e-10) {
716 return z2;
717 }
718
719 final double delta = 2.0 * x - 1.0;
720 final double g1 = z3 + z2;
721 final double g2 = z3 - z2;
722 final double g3 = z4 + z1;
723 final double g4 = (z4 - z1) / 3.0;
724 final double a0 = 9.0 * g1 - g3;
725 final double a1 = 9.0 * g2 - g4;
726 final double a2 = g3 - g1;
727 final double a3 = g4 - g2;
728 final double zx = 0.0625 * (a0 + delta * (a1 + delta * (a2 + delta * a3)));
729
730 return zx;
731 }
732
733 /**
734 * Allows smooth joining of functions f1 and f2
735 * (i.e. continuous first derivatives) at origin.
736 * <p>
737 * This function, describe in section F.2.12.1 of the reference document, is
738 * recommanded for computational efficiency.
739 * </p>
740 * @param dF1 first function
741 * @param dF2 second function
742 * @param dA width of transition region
743 * @param dX x value
744 * @return the computed value
745 */
746 private double join(final double dF1, final double dF2,
747 final double dA, final double dX) {
748 final double ee = clipExp(dA * dX);
749 return (dF1 * ee + dF2) / (ee + 1.0);
750 }
751
752 /**
753 * The Epstein function.
754 * <p>
755 * This function, describe in section 2.5.1 of the reference document, is used
756 * as a basis analytical function in NeQuick for the construction of the ionospheric layers.
757 * </p>
758 * @param x x parameter
759 * @param y y parameter
760 * @param z z parameter
761 * @param w w parameter
762 * @return value of the epstein function
763 */
764 private double epst(final double x, final double y,
765 final double z, final double w) {
766 final double ex = clipExp((w - y) / z);
767 final double opex = 1.0 + ex;
768 final double epst = x * ex / (opex * opex);
769 return epst;
770 }
771
772 }