1 /* Copyright 2002-2025 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.nequick;
18
19 import org.hipparchus.util.FastMath;
20 import org.hipparchus.util.SinCos;
21 import org.orekit.time.DateComponents;
22 import org.orekit.time.DateTimeComponents;
23 import org.orekit.time.TimeComponents;
24
25 /**
26 * This class performs the computation of the parameters used by the NeQuick model.
27 *
28 * @author Bryan Cazabonne
29 *
30 * @see "European Union (2016). European GNSS (Galileo) Open Service-Ionospheric Correction
31 * Algorithm for Galileo Single Frequency Users. 1.2."
32 * @see <a href="https://www.itu.int/rec/R-REC-P.531/en">ITU-R P.531</a>
33 *
34 * @since 10.1
35 */
36 public class NeQuickParameters {
37
38 /** Solar zenith angle at day night transition, degrees. */
39 private static final double X0 = 86.23292796211615;
40
41 /** Current date time components.
42 * @since 13.0
43 */
44 private final DateTimeComponents dateTime;
45
46 /** Effective sunspot number.
47 * @since 13.0
48 */
49 private final double azr;
50
51 /** F2 layer critical frequency.
52 * @since 13.0
53 */
54 private final double foF2;
55
56 /** F2 layer maximum density. */
57 private final double nmF2;
58
59 /** F2 layer maximum density height [km]. */
60 private final double hmF2;
61
62 /** F1 layer maximum density height [km]. */
63 private final double hmF1;
64
65 /** E layer maximum density height [km]. */
66 private final double hmE;
67
68 /** F2 layer bottom thickness parameter [km]. */
69 private final double b2Bot;
70
71 /** F1 layer top thickness parameter [km]. */
72 private final double b1Top;
73
74 /** F1 layer bottom thickness parameter [km]. */
75 private final double b1Bot;
76
77 /** E layer top thickness parameter [km]. */
78 private final double beTop;
79
80 /** E layer bottom thickness parameter [km]. */
81 private final double beBot;
82
83 /** Layer amplitudes. */
84 private final double[] amplitudes;
85
86 /**
87 * Build a new instance.
88 * @param dateTime current date time components
89 * @param flattenF2 F2 coefficients used by the F2 layer (flatten array)
90 * @param flattenFm3 Fm3 coefficients used by the M(3000)F2 layer (flatten array)
91 * @param latitude latitude of a point along the integration path, in radians
92 * @param longitude longitude of a point along the integration path, in radians
93 * @param az effective ionisation level
94 * @param modip modip
95 */
96 public NeQuickParameters(final DateTimeComponents dateTime, final double[] flattenF2, final double[] flattenFm3,
97 final double latitude, final double longitude, final double az, final double modip) {
98
99 this.dateTime = dateTime;
100
101 // Effective sunspot number (Eq. 19)
102 azr = FastMath.sqrt(167273.0 + (az - 63.7) * 1123.6) - 408.99;
103
104 // Date and Time components
105 final DateComponents date = dateTime.getDate();
106 final TimeComponents time = dateTime.getTime();
107 // Hours
108 final double hours = time.getSecondsInUTCDay() / 3600.0;
109 // Effective solar zenith angle in radians
110 final double xeff = computeEffectiveSolarAngle(date.getMonth(), hours, latitude, longitude);
111
112 // E layer maximum density height in km (Eq. 78)
113 this.hmE = 120.0;
114 // E layer critical frequency in MHz
115 final double foE = computefoE(date.getMonth(), az, xeff, latitude);
116 // E layer maximum density in 10^11 m⁻³ (Eq. 36)
117 final double nmE = 0.124 * foE * foE;
118
119 // Time argument (Eq. 49)
120 final double t = FastMath.toRadians(15 * hours) - FastMath.PI;
121 // Compute Fourier time series for foF2 and M(3000)F2
122 final double[] scT = sinCos(t, 6);
123 final double[] cf2 = computeCF2(flattenF2, scT);
124 final double[] cm3 = computeCm3(flattenFm3, scT);
125 // F2 layer critical frequency in MHz
126 final double[] scL = sinCos(longitude, 8);
127 this.foF2 = computefoF2(modip, cf2, latitude, scL);
128 // Maximum Usable Frequency factor
129 final double mF2 = computeMF2(modip, cm3, latitude, scL);
130 // F2 layer maximum density in 10^11 m⁻³
131 this.nmF2 = 0.124 * foF2 * foF2;
132 // F2 layer maximum density height in km
133 this.hmF2 = computehmF2(foE, mF2);
134
135 // F1 layer critical frequency in MHz
136 final double foF1 = computefoF1(foE);
137 // F1 layer maximum density in 10^11 m⁻³
138 final double nmF1;
139 if (foF1 <= 0.0 && foE > 2.0) {
140 final double foEpopf = foE + 0.5;
141 nmF1 = 0.124 * foEpopf * foEpopf;
142 } else {
143 nmF1 = 0.124 * foF1 * foF1;
144 }
145 // F1 layer maximum density height in km
146 this.hmF1 = 0.5 * (hmF2 + hmE);
147
148 // Thickness parameters (Eq. 85 to 89)
149 final double a = 0.01 * clipExp(-3.467 + 0.857 * FastMath.log(foF2 * foF2) + 2.02 * FastMath.log(mF2));
150 this.b2Bot = 0.385 * nmF2 / a;
151 this.b1Top = 0.3 * (hmF2 - hmF1);
152 this.b1Bot = 0.5 * (hmF1 - hmE);
153 this.beTop = FastMath.max(b1Bot, 7.0);
154 this.beBot = 5.0;
155
156 // Layer amplitude coefficients
157 this.amplitudes = computeLayerAmplitudes(nmE, nmF1, foF1);
158
159 }
160
161 /**
162 * Get current date time components.
163 * @return current date time components
164 * @since 13.0
165 */
166 public DateTimeComponents getDateTime() {
167 return dateTime;
168 }
169
170 /**
171 * Get effective sunspot number.
172 * @return effective sunspot number
173 * @since 13.0
174 */
175 public double getAzr() {
176 return azr;
177 }
178
179 /**
180 * Get F2 layer critical frequency.
181 * @return F2 layer critical frequency
182 * @since 13.0
183 */
184 public double getFoF2() {
185 return foF2;
186 }
187
188 /**
189 * Get the F2 layer maximum density.
190 * @return nmF2
191 */
192 public double getNmF2() {
193 return nmF2;
194 }
195
196 /**
197 * Get the F2 layer maximum density height.
198 * @return hmF2 in km
199 */
200 public double getHmF2() {
201 return hmF2;
202 }
203
204 /**
205 * Get the F1 layer maximum density height.
206 * @return hmF1 in km
207 */
208 public double getHmF1() {
209 return hmF1;
210 }
211
212 /**
213 * Get the E layer maximum density height.
214 * @return hmE in km
215 */
216 public double getHmE() {
217 return hmE;
218 }
219
220 /**
221 * Get the F2 layer thickness parameter (bottom).
222 * @return B2Bot in km
223 */
224 public double getB2Bot() {
225 return b2Bot;
226 }
227
228 /**
229 * Get the F1 layer thickness parameter.
230 * @param h current height (km)
231 * @return B1 in km
232 * @since 13.0
233 */
234 public double getBF1(final double h) {
235 // Eq. 110
236 return (h > hmF1) ? b1Top : b1Bot;
237 }
238
239 /**
240 * Get the E layer thickness parameter.
241 * @param h current height (km)
242 * @return Be in km
243 * @since 13.0
244 */
245 public double getBE(final double h) {
246 // Eq. 109
247 return (h > hmE) ? beTop : beBot;
248 }
249
250 /**
251 * Get the F2, F1 and E layer amplitudes.
252 * <p>
253 * The resulting element is an array having the following form:
254 * <ul>
255 * <li>double[0] = A1 → F2 layer amplitude
256 * <li>double[1] = A2 → F1 layer amplitude
257 * <li>double[2] = A3 → E layer amplitude
258 * </ul>
259 * @return layer amplitudes
260 */
261 public double[] getLayerAmplitudes() {
262 return amplitudes.clone();
263 }
264
265 /**
266 * This method computes the effective solar zenith angle.
267 * <p>
268 * The effective solar zenith angle is compute as a function of the
269 * solar zenith angle and the solar zenith angle at day night transition.
270 * </p>
271 * @param month current month of the year
272 * @param hours universal time (hours)
273 * @param latitude in radians
274 * @param longitude in radians
275 * @return the effective solar zenith angle, radians
276 */
277 private double computeEffectiveSolarAngle(final int month,
278 final double hours,
279 final double latitude,
280 final double longitude) {
281 // Local time (Eq.4)
282 final double lt = hours + longitude / FastMath.toRadians(15.0);
283 // Day of year at the middle of the month (Eq. 20)
284 final double dy = 30.5 * month - 15.0;
285 // Time (Eq. 21)
286 final double t = dy + (18 - hours) / 24;
287 // Arguments am and al (Eq. 22 and 23)
288 final double am = FastMath.toRadians(0.9856 * t - 3.289);
289 final double al = am + FastMath.toRadians(1.916 * FastMath.sin(am) + 0.020 * FastMath.sin(2.0 * am) + 282.634);
290 // Sine and cosine of solar declination (Eq. 24 and 25)
291 final double sDec = 0.39782 * FastMath.sin(al);
292 final double cDec = FastMath.sqrt(1. - sDec * sDec);
293 // Solar zenith angle, deg (Eq. 26 and 27)
294 final SinCos scLat = FastMath.sinCos(latitude);
295 final double coef = (FastMath.PI / 12) * (12 - lt);
296 final double cZenith = scLat.sin() * sDec + scLat.cos() * cDec * FastMath.cos(coef);
297 final double angle = FastMath.atan2(FastMath.sqrt(1.0 - cZenith * cZenith), cZenith);
298 final double x = FastMath.toDegrees(angle);
299 // Effective solar zenith angle (Eq. 28)
300 final double xeff = join(90.0 - 0.24 * clipExp(20.0 - 0.2 * x), x, 12.0, x - X0);
301 return FastMath.toRadians(xeff);
302 }
303
304 /**
305 * This method computes the E layer critical frequency at a given location.
306 * @param month current month
307 * @param az ffective ionisation level
308 * @param xeff effective solar zenith angle in radians
309 * @param latitude latitude in radians
310 * @return the E layer critical frequency at a given location in MHz
311 */
312 private double computefoE(final int month, final double az,
313 final double xeff, final double latitude) {
314 // The latitude has to be converted in degrees
315 final double lat = FastMath.toDegrees(latitude);
316 // Square root of the effective ionisation level
317 final double sqAz = FastMath.sqrt(az);
318 // seas parameter (Eq. 30 to 32)
319 final int seas;
320 if (month == 1 || month == 2 || month == 11 || month == 12) {
321 seas = -1;
322 } else if (month == 3 || month == 4 || month == 9 || month == 10) {
323 seas = 0;
324 } else {
325 seas = 1;
326 }
327 // Latitudinal dependence (Eq. 33 and 34)
328 final double ee = clipExp(0.3 * lat);
329 final double seasp = seas * ((ee - 1.0) / (ee + 1.0));
330 // Critical frequency (Eq. 35)
331 final double coef = 1.112 - 0.019 * seasp;
332 return FastMath.sqrt(coef * coef * sqAz * FastMath.pow(FastMath.cos(xeff), 0.6) + 0.49);
333
334 }
335
336 /**
337 * Computes the F2 layer height of maximum electron density.
338 * @param foE E layer layer critical frequency in MHz
339 * @param mF2 maximum usable frequency factor
340 * @return hmF2 in km
341 */
342 private double computehmF2(final double foE, final double mF2) {
343 // Ratio
344 final double fo = foF2 / foE;
345 final double ratio = join(fo, 1.75, 20.0, fo - 1.75);
346
347 // deltaM parameter
348 double deltaM = -0.012;
349 if (foE >= 1e-30) {
350 deltaM += 0.253 / (ratio - 1.215);
351 }
352
353 // hmF2 Eq. 80
354 final double mF2Sq = mF2 * mF2;
355 final double temp = FastMath.sqrt((0.0196 * mF2Sq + 1) / (1.2967 * mF2Sq - 1.0));
356 return ((1490.0 * mF2 * temp) / (mF2 + deltaM)) - 176.0;
357
358 }
359
360 /** Compute sines and cosines.
361 * @param a argument
362 * @param n number of terms
363 * @return sin(a), cos(a), sin(2a), cos(2a) … sin(n a), cos(n a) array
364 * @since 12.1.3
365 */
366 private double[] sinCos(final double a, final int n) {
367
368 final SinCos sc0 = FastMath.sinCos(a);
369 SinCos sci = sc0;
370 final double[] sc = new double[2 * n];
371 int isc = 0;
372 sc[isc++] = sci.sin();
373 sc[isc++] = sci.cos();
374 for (int i = 1; i < n; i++) {
375 sci = SinCos.sum(sc0, sci);
376 sc[isc++] = sci.sin();
377 sc[isc++] = sci.cos();
378 }
379
380 return sc;
381
382 }
383
384 /**
385 * Computes cf2 coefficients.
386 * @param flattenF2 F2 coefficients used by the F2 layer (flatten array)
387 * @param scT sines/cosines array of time argument
388 * @return the cf2 coefficients array
389 */
390 private double[] computeCF2(final double[] flattenF2, final double[] scT) {
391
392 // interpolation coefficients for effective spot number
393 final double azr01 = azr * 0.01;
394 final double omazr01 = 1 - azr01;
395
396 // Eq. 44 and Eq. 50 merged into one loop
397 final double[] cf2 = new double[76];
398 int index = 0;
399 for (int i = 0; i < cf2.length; i++) {
400 cf2[i] = omazr01 * flattenF2[index ] + azr01 * flattenF2[index + 1] +
401 (omazr01 * flattenF2[index + 2] + azr01 * flattenF2[index + 3]) * scT[ 0] +
402 (omazr01 * flattenF2[index + 4] + azr01 * flattenF2[index + 5]) * scT[ 1] +
403 (omazr01 * flattenF2[index + 6] + azr01 * flattenF2[index + 7]) * scT[ 2] +
404 (omazr01 * flattenF2[index + 8] + azr01 * flattenF2[index + 9]) * scT[ 3] +
405 (omazr01 * flattenF2[index + 10] + azr01 * flattenF2[index + 11]) * scT[ 4] +
406 (omazr01 * flattenF2[index + 12] + azr01 * flattenF2[index + 13]) * scT[ 5] +
407 (omazr01 * flattenF2[index + 14] + azr01 * flattenF2[index + 15]) * scT[ 6] +
408 (omazr01 * flattenF2[index + 16] + azr01 * flattenF2[index + 17]) * scT[ 7] +
409 (omazr01 * flattenF2[index + 18] + azr01 * flattenF2[index + 19]) * scT[ 8] +
410 (omazr01 * flattenF2[index + 20] + azr01 * flattenF2[index + 21]) * scT[ 9] +
411 (omazr01 * flattenF2[index + 22] + azr01 * flattenF2[index + 23]) * scT[10] +
412 (omazr01 * flattenF2[index + 24] + azr01 * flattenF2[index + 25]) * scT[11];
413 index += 26;
414 }
415 return cf2;
416 }
417
418 /**
419 * Computes Cm3 coefficients.
420 * @param flattenFm3 Fm3 coefficients used by the M(3000)F2 layer (flatten array)
421 * @param scT sines/cosines array of time argument
422 * @return the Cm3 coefficients array
423 */
424 private double[] computeCm3(final double[] flattenFm3, final double[] scT) {
425
426 // interpolation coefficients for effective spot number
427 final double azr01 = azr * 0.01;
428 final double omazr01 = 1 - azr01;
429
430 // Eq. 44 and Eq. 51 merged into one loop
431 final double[] cm3 = new double[49];
432 int index = 0;
433 for (int i = 0; i < cm3.length; i++) {
434 cm3[i] = omazr01 * flattenFm3[index ] + azr01 * flattenFm3[index + 1] +
435 (omazr01 * flattenFm3[index + 2] + azr01 * flattenFm3[index + 3]) * scT[ 0] +
436 (omazr01 * flattenFm3[index + 4] + azr01 * flattenFm3[index + 5]) * scT[ 1] +
437 (omazr01 * flattenFm3[index + 6] + azr01 * flattenFm3[index + 7]) * scT[ 2] +
438 (omazr01 * flattenFm3[index + 8] + azr01 * flattenFm3[index + 9]) * scT[ 3] +
439 (omazr01 * flattenFm3[index + 10] + azr01 * flattenFm3[index + 11]) * scT[ 4] +
440 (omazr01 * flattenFm3[index + 12] + azr01 * flattenFm3[index + 13]) * scT[ 5] +
441 (omazr01 * flattenFm3[index + 14] + azr01 * flattenFm3[index + 15]) * scT[ 6] +
442 (omazr01 * flattenFm3[index + 16] + azr01 * flattenFm3[index + 17]) * scT[ 7];
443 index += 18;
444 }
445 return cm3;
446 }
447
448 /**
449 * This method computes the F2 layer critical frequency.
450 * @param modip modified DIP latitude, in degrees
451 * @param cf2 Fourier time series for foF2
452 * @param latitude latitude in radians
453 * @param scL sines/cosines array of longitude argument
454 * @return the F2 layer critical frequency, MHz
455 */
456 private double computefoF2(final double modip, final double[] cf2,
457 final double latitude, final double[] scL) {
458
459 // Legendre grades (Eq. 63)
460 final int[] q = new int[] {
461 12, 12, 9, 5, 2, 1, 1, 1, 1
462 };
463
464 double frequency = cf2[0];
465
466 // ModipGrid coefficients Eq. 57
467 final double sinMODIP = FastMath.sin(FastMath.toRadians(modip));
468 final double[] m = new double[12];
469 m[0] = 1.0;
470 for (int i = 1; i < q[0]; i++) {
471 m[i] = sinMODIP * m[i - 1];
472 frequency += m[i] * cf2[i];
473 }
474
475 // latitude and longitude terms
476 int index = 12;
477 final double cosLat1 = FastMath.cos(latitude);
478 double cosLatI = cosLat1;
479 for (int i = 1; i < q.length; i++) {
480 final double c = cosLatI * scL[2 * i - 1];
481 final double s = cosLatI * scL[2 * i - 2];
482 for (int j = 0; j < q[i]; j++) {
483 frequency += m[j] * c * cf2[index++];
484 frequency += m[j] * s * cf2[index++];
485 }
486 cosLatI *= cosLat1; // Eq. 58
487 }
488
489 return frequency;
490
491 }
492
493 /**
494 * This method computes the Maximum Usable Frequency factor.
495 * @param modip modified DIP latitude, in degrees
496 * @param cm3 Fourier time series for M(3000)F2
497 * @param latitude latitude in radians
498 * @param scL sines/cosines array of longitude argument
499 * @return the Maximum Usable Frequency factor
500 */
501 private double computeMF2(final double modip, final double[] cm3,
502 final double latitude, final double[] scL) {
503
504 // Legendre grades (Eq. 71)
505 final int[] r = new int[] {
506 7, 8, 6, 3, 2, 1, 1
507 };
508
509 double m3000 = cm3[0];
510
511 // ModipGrid coefficients Eq. 57
512 final double sinMODIP = FastMath.sin(FastMath.toRadians(modip));
513 final double[] m = new double[12];
514 m[0] = 1.0;
515 for (int i = 1; i < 12; i++) {
516 m[i] = sinMODIP * m[i - 1];
517 if (i < 7) {
518 m3000 += m[i] * cm3[i];
519 }
520 }
521
522 // latitude and longitude terms
523 int index = 7;
524 final double cosLat1 = FastMath.cos(latitude);
525 double cosLatI = cosLat1;
526 for (int i = 1; i < r.length; i++) {
527 final double c = cosLatI * scL[2 * i - 1];
528 final double s = cosLatI * scL[2 * i - 2];
529 for (int j = 0; j < r[i]; j++) {
530 m3000 += m[j] * c * cm3[index++];
531 m3000 += m[j] * s * cm3[index++];
532 }
533 cosLatI *= cosLat1; // Eq. 58
534 }
535
536 return m3000;
537
538 }
539
540 /**
541 * This method computes the F1 layer critical frequency.
542 * <p>
543 * This computation performs the algorithm exposed in Annex F
544 * of the reference document.
545 * </p>
546 * @param foE the E layer critical frequency, MHz
547 * @return the F1 layer critical frequency, MHz
548 */
549 private double computefoF1(final double foE) {
550 final double temp = join(1.4 * foE, 0.0, 1000.0, foE - 2.0);
551 final double temp2 = join(0.0, temp, 1000.0, foE - temp);
552 final double value = join(temp2, 0.85 * temp2, 60.0, 0.85 * foF2 - temp2);
553 if (value < 1.0E-6) {
554 return 0.0;
555 } else {
556 return value;
557 }
558 }
559
560 /**
561 * This method allows the computation of the F2, F1 and E layer amplitudes.
562 * <p>
563 * The resulting element is an array having the following form:
564 * <ul>
565 * <li>double[0] = A1 → F2 layer amplitude
566 * <li>double[1] = A2 → F1 layer amplitude
567 * <li>double[2] = A3 → E layer amplitude
568 * </ul>
569 * </p>
570 * @param nmE E layer maximum density in 10^11 m⁻³
571 * @param nmF1 F1 layer maximum density in 10^11 m⁻³
572 * @param foF1 F1 layer critical frequency in MHz
573 * @return a three components array containing the layer amplitudes
574 */
575 private double[] computeLayerAmplitudes(final double nmE, final double nmF1, final double foF1) {
576 // Initialize array
577 final double[] amplitude = new double[3];
578
579 // F2 layer amplitude (Eq. 90)
580 final double a1 = 4.0 * nmF2;
581 amplitude[0] = a1;
582
583 // F1 and E layer amplitudes (Eq. 91 to 98)
584 if (foF1 < 0.5) {
585 amplitude[1] = 0.0;
586 amplitude[2] = 4.0 * (nmE - epst(a1, hmF2, b2Bot, hmE));
587 } else {
588 double a2a = 0.0;
589 double a3a = 4.0 * nmE;
590 for (int i = 0; i < 5; i++) {
591 a2a = 4.0 * (nmF1 - epst(a1, hmF2, b2Bot, hmF1) - epst(a3a, hmE, beTop, hmF1));
592 a2a = join(a2a, 0.8 * nmF1, 1.0, a2a - 0.8 * nmF1);
593 a3a = 4.0 * (nmE - epst(a2a, hmF1, b1Bot, hmE) - epst(a1, hmF2, b2Bot, hmE));
594 }
595 amplitude[1] = a2a;
596 amplitude[2] = join(a3a, 0.05, 60.0, a3a - 0.005);
597 }
598
599 return amplitude;
600 }
601
602 /**
603 * A clipped exponential function.
604 * <p>
605 * This function, describe in section F.2.12.2 of the reference document, is
606 * recommanded for the computation of exponential values.
607 * </p>
608 * @param power power for exponential function
609 * @return clipped exponential value
610 */
611 private double clipExp(final double power) {
612 if (power > 80.0) {
613 return 5.5406E34;
614 } else if (power < -80) {
615 return 1.8049E-35;
616 } else {
617 return FastMath.exp(power);
618 }
619 }
620
621 /**
622 * Allows smooth joining of functions f1 and f2
623 * (i.e. continuous first derivatives) at origin.
624 * <p>
625 * This function, describe in section F.2.12.1 of the reference document, is
626 * recommended for computational efficiency.
627 * </p>
628 * @param dF1 first function
629 * @param dF2 second function
630 * @param dA width of transition region
631 * @param dX x value
632 * @return the computed value
633 */
634 double join(final double dF1, final double dF2, final double dA, final double dX) {
635 final double ee = clipExp(dA * dX);
636 return (dF1 * ee + dF2) / (ee + 1.0);
637 }
638
639 /**
640 * The Epstein function.
641 * <p>
642 * This function, describe in section 2.5.1 of the reference document, is used
643 * as a basis analytical function in NeQuick for the construction of the ionospheric layers.
644 * </p>
645 * @param x x parameter
646 * @param y y parameter
647 * @param z z parameter
648 * @param w w parameter
649 * @return value of the epstein function
650 */
651 private double epst(final double x, final double y,
652 final double z, final double w) {
653 final double ex = clipExp((w - y) / z);
654 final double opex = 1.0 + ex;
655 return x * ex / (opex * opex);
656
657 }
658
659 }