NeQuickParameters.java

/* Copyright 2002-2023 CS GROUP
 * Licensed to CS GROUP (CS) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * CS licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.orekit.models.earth.ionosphere;

import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathArrays;
import org.hipparchus.util.SinCos;
import org.orekit.time.DateComponents;
import org.orekit.time.DateTimeComponents;
import org.orekit.time.TimeComponents;

/**
 * This class perfoms the computation of the parameters used by the NeQuick model.
 *
 * @author Bryan Cazabonne
 *
 * @see "European Union (2016). European GNSS (Galileo) Open Service-Ionospheric Correction
 *       Algorithm for Galileo Single Frequency Users. 1.2."
 *
 * @since 10.1
 */
class NeQuickParameters {

    /** Solar zenith angle at day night transition, degrees. */
    private static final double X0 = 86.23292796211615;

    /** F2 layer maximum density. */
    private final double nmF2;

    /** F2 layer maximum density height [km]. */
    private final double hmF2;

    /** F1 layer maximum density height [km]. */
    private final double hmF1;

    /** E layer maximum density height [km]. */
    private final double hmE;

    /** F2 layer bottom thickness parameter [km]. */
    private final double b2Bot;

    /** F1 layer top thickness parameter [km]. */
    private final double b1Top;

    /** F1 layer bottom thickness parameter [km]. */
    private final double b1Bot;

    /** E layer top thickness parameter [km]. */
    private final double beTop;

    /** E layer bottom thickness parameter [km]. */
    private final double beBot;

    /** topside thickness parameter [km]. */
    private final double h0;

    /** Layer amplitudes. */
    private final double[] amplitudes;

    /**
     * Build a new instance.
     * @param dateTime current date time components
     * @param f2 F2 coefficients used by the F2 layer
     * @param fm3 Fm3 coefficients used by the F2 layer
     * @param latitude latitude of a point along the integration path, in radians
     * @param longitude longitude of a point along the integration path, in radians
     * @param alpha effective ionisation level coefficients
     * @param modipGrip modip grid
     */
    NeQuickParameters(final DateTimeComponents dateTime, final double[][][] f2,
                      final double[][][] fm3, final double latitude, final double longitude,
                      final double[] alpha, final double[][] modipGrip) {

        // MODIP in degrees
        final double modip = computeMODIP(latitude, longitude, modipGrip);
        // Effective ionisation level Az
        final double az = computeAz(modip, alpha);
        // Effective sunspot number (Eq. 19)
        final double azr = FastMath.sqrt(167273.0 + (az - 63.7) * 1123.6) - 408.99;
        // Date and Time components
        final DateComponents date = dateTime.getDate();
        final TimeComponents time = dateTime.getTime();
        // Hours
        final double hours  = time.getSecondsInUTCDay() / 3600.0;
        // Effective solar zenith angle in radians
        final double xeff = computeEffectiveSolarAngle(date.getMonth(), hours, latitude, longitude);

        // Coefficients for F2 layer parameters
        // Compute the array of interpolated coefficients for foF2 (Eq. 44)
        final double[][] af2 = new double[76][13];
        for (int j = 0; j < 76; j++) {
            for (int k = 0; k < 13; k++ ) {
                af2[j][k] = f2[0][j][k] * (1.0 - (azr * 0.01)) + f2[1][j][k] * (azr * 0.01);
            }
        }

        // Compute the array of interpolated coefficients for M(3000)F2 (Eq. 46)
        final  double[][] am3 = new double[49][9];
        for (int j = 0; j < 49; j++) {
            for (int k = 0; k < 9; k++ ) {
                am3[j][k] = fm3[0][j][k] * (1.0 - (azr * 0.01)) + fm3[1][j][k] * (azr * 0.01);
            }
        }

        // E layer maximum density height in km (Eq. 78)
        this.hmE = 120.0;
        // E layer critical frequency in MHz
        final double foE = computefoE(date.getMonth(), az, xeff, latitude);
        // E layer maximum density in 10^11 m-3 (Eq. 36)
        final double nmE = 0.124 * foE * foE;

        // Time argument (Eq. 49)
        final double t = FastMath.toRadians(15 * hours) - FastMath.PI;
        // Compute Fourier time series for foF2 and M(3000)F2
        final double[] cf2 = computeCF2(af2, t);
        final double[] cm3 = computeCm3(am3, t);
        // F2 layer critical frequency in MHz
        final double foF2 = computefoF2(modip, cf2, latitude, longitude);
        // Maximum Usable Frequency factor
        final double mF2  = computeMF2(modip, cm3, latitude, longitude);
        // F2 layer maximum density in 10^11 m-3
        this.nmF2 = 0.124 * foF2 * foF2;
        // F2 layer maximum density height in km
        this.hmF2 = computehmF2(foE, foF2, mF2);

        // F1 layer critical frequency in MHz
        final double foF1 = computefoF1(foE, foF2);
        // F1 layer maximum density in 10^11 m-3
        final double nmF1;
        if (foF1 <= 0.0 && foE > 2.0) {
            final double foEpopf = foE + 0.5;
            nmF1 = 0.124 * foEpopf * foEpopf;
        } else {
            nmF1 = 0.124 * foF1 * foF1;
        }
        // F1 layer maximum density height in km
        this.hmF1 = 0.5 * (hmF2 + hmE);

        // Thickness parameters (Eq. 85 to 89)
        final double a = 0.01 * clipExp(-3.467 + 0.857 * FastMath.log(foF2 * foF2) + 2.02 * FastMath.log(mF2));
        this.b2Bot = 0.385 * nmF2 / a;
        this.b1Top = 0.3 * (hmF2 - hmF1);
        this.b1Bot = 0.5 * (hmF1 - hmE);
        this.beTop = FastMath.max(b1Bot, 7.0);
        this.beBot = 5.0;

        // Layer amplitude coefficients
        this.amplitudes = computeLayerAmplitudes(nmE, nmF1, foF1);

        // Topside thickness parameter
        this.h0 = computeH0(date.getMonth(), azr);
    }

    /**
     * Get the F2 layer maximum density.
     * @return nmF2
     */
    public double getNmF2() {
        return nmF2;
    }

    /**
     * Get the F2 layer maximum density height.
     * @return hmF2 in km
     */
    public double getHmF2() {
        return hmF2;
    }

    /**
     * Get the F1 layer maximum density height.
     * @return hmF1 in km
     */
    public double getHmF1() {
        return hmF1;
    }

    /**
     * Get the E layer maximum density height.
     * @return hmE in km
     */
    public double getHmE() {
        return hmE;
    }

    /**
     * Get the F2 layer thickness parameter (bottom).
     * @return B2Bot in km
     */
    public double getB2Bot() {
        return b2Bot;
    }

    /**
     * Get the F1 layer thickness parameter (top).
     * @return B1Top in km
     */
    public double getB1Top() {
        return b1Top;
    }

    /**
     * Get the F1 layer thickness parameter (bottom).
     * @return B1Bot in km
     */
    public double getB1Bot() {
        return b1Bot;
    }

    /**
     * Get the E layer thickness parameter (bottom).
     * @return BeBot in km
     */
    public double getBEBot() {
        return beBot;
    }

    /**
     * Get the E layer thickness parameter (top).
     * @return BeTop in km
     */
    public double getBETop() {
        return beTop;
    }

    /**
     * Get the F2, F1 and E layer amplitudes.
     * <p>
     * The resulting element is an array having the following form:
     * <ul>
     * <li>double[0] = A1 → F2 layer amplitude
     * <li>double[1] = A2 → F1 layer amplitude
     * <li>double[2] = A3 → E  layer amplitude
     * </ul>
     * @return layer amplitudes
     */
    public double[] getLayerAmplitudes() {
        return amplitudes.clone();
    }

    /**
     * Get the topside thickness parameter H0.
     * @return H0 in km
     */
    public double getH0() {
        return h0;
    }

    /**
     * Computes the value of the modified dip latitude (MODIP) for the
     * given latitude and longitude.
     *
     * @param lat receiver latitude, radians
     * @param lon receiver longitude, radians
     * @param stModip modip grid
     * @return the MODIP in degrees
     */
    private double computeMODIP(final double lat, final double lon, final double[][] stModip) {

        // For the MODIP computation, latitude and longitude have to be converted in degrees
        final double latitude = FastMath.toDegrees(lat);
        final double longitude = FastMath.toDegrees(lon);

        // Extreme cases
        if (latitude == 90.0 || latitude == -90.0) {
            return latitude;
        }

        // Auxiliary parameter l (Eq. 6 to 8)
        final int lF = (int) ((longitude + 180) * 0.1);
        int l = lF - 2;
        if (l < -2) {
            l += 36;
        } else if (l > 33) {
            l -= 36;
        }

        // Auxiliary parameter a (Eq. 9 to 11)
        final double a  = 0.2 * (latitude + 90) + 1.0;
        final double aF = FastMath.floor(a);
        // Eq. 10
        final double x = a - aF;
        // Eq. 11
        final int i = (int) aF - 2;

        // zi coefficients (Eq. 12 and 13)
        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);
        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);
        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);
        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);

        // Auxiliary parameter b (Eq. 14 and 15)
        final double b  = (longitude + 180) * 0.1;
        final double bF = FastMath.floor(b);
        final double y  = b - bF;

        // MODIP (Ref Eq. 16)
        final double modip = interpolate(z1, z2, z3, z4, y);

        return modip;
    }

    /**
     * This method computes the effective ionisation level Az.
     * <p>
     * This parameter is used for the computation of the Total Electron Content (TEC).
     * </p>
     * @param modip modified dip latitude (MODIP) in degrees
     * @param alpha effective ionisation level coefficients
     * @return the ionisation level Az
     */
    private double computeAz(final double modip, final double[] alpha) {
        // Particular condition (Eq. 17)
        if (alpha[0] == 0.0 && alpha[1] == 0.0 && alpha[2] == 0.0) {
            return 63.7;
        }
        // Az = a0 + modip * a1 + modip^2 * a2 (Eq. 18)
        double az = alpha[0] + modip * (alpha[1] + modip * alpha[2]);
        // If Az < 0 -> Az = 0
        az = FastMath.max(0.0, az);
        // If Az > 400 -> Az = 400
        az = FastMath.min(400.0, az);
        return az;
    }

    /**
     * This method computes the effective solar zenith angle.
     * <p>
     * The effective solar zenith angle is compute as a function of the
     * solar zenith angle and the solar zenith angle at day night transition.
     * </p>
     * @param month current month of the year
     * @param hours universal time (hours)
     * @param latitude in radians
     * @param longitude in radians
     * @return the effective solar zenith angle, radians
     */
    private double computeEffectiveSolarAngle(final int month,
                                              final double hours,
                                              final double latitude,
                                              final double longitude) {
        // Local time (Eq.4)
        final double lt = hours + longitude / FastMath.toRadians(15.0);
        // Day of year at the middle of the month (Eq. 20)
        final double dy = 30.5 * month - 15.0;
        // Time (Eq. 21)
        final double t = dy + (18 - hours) / 24;
        // Arguments am and al (Eq. 22 and 23)
        final double am = FastMath.toRadians(0.9856 * t - 3.289);
        final double al = am + FastMath.toRadians(1.916 * FastMath.sin(am) + 0.020 * FastMath.sin(2.0 * am) + 282.634);
        // Sine and cosine of solar declination (Eq. 24 and 25)
        final double sDec = 0.39782 * FastMath.sin(al);
        final double cDec = FastMath.sqrt(1. - sDec * sDec);
        // Solar zenith angle, deg (Eq. 26 and 27)
        final SinCos scLat   = FastMath.sinCos(latitude);
        final double coef    = (FastMath.PI / 12) * (12 - lt);
        final double cZenith = scLat.sin() * sDec + scLat.cos() * cDec * FastMath.cos(coef);
        final double angle   = FastMath.atan2(FastMath.sqrt(1.0 - cZenith * cZenith), cZenith);
        final double x       = FastMath.toDegrees(angle);
        // Effective solar zenith angle (Eq. 28)
        final double xeff = join(90.0 - 0.24 * clipExp(20.0 - 0.2 * x), x, 12.0, x - X0);
        return FastMath.toRadians(xeff);
    }

    /**
     * This method computes the E layer critical frequency at a given location.
     * @param month current month
     * @param az ffective ionisation level
     * @param xeff effective solar zenith angle in radians
     * @param latitude latitude in radians
     * @return the E layer critical frequency at a given location in MHz
     */
    private double computefoE(final int month, final double az,
                              final double xeff, final double latitude) {
        // The latitude has to be converted in degrees
        final double lat = FastMath.toDegrees(latitude);
        // Square root of the effective ionisation level
        final double sqAz = FastMath.sqrt(az);
        // seas parameter (Eq. 30 to 32)
        final int seas;
        if (month == 1 || month == 2 || month == 11 || month == 12) {
            seas = -1;
        } else if (month == 3 || month == 4 || month == 9 || month == 10) {
            seas = 0;
        } else {
            seas = 1;
        }
        // Latitudinal dependence (Eq. 33 and 34)
        final double ee = clipExp(0.3 * lat);
        final double seasp = seas * ((ee - 1.0) / (ee + 1.0));
        // Critical frequency (Eq. 35)
        final double coef = 1.112 - 0.019 * seasp;
        final double foE = FastMath.sqrt(coef * coef * sqAz * FastMath.pow(FastMath.cos(xeff), 0.6) + 0.49);
        return foE;
    }

    /**
     * Computes the F2 layer height of maximum electron density.
     * @param foE E layer layer critical frequency in MHz
     * @param foF2 F2 layer layer critical frequency in MHz
     * @param mF2 maximum usable frequency factor
     * @return hmF2 in km
     */
    private double computehmF2(final double foE, final double foF2, final double mF2) {
        // Ratio
        final double fo = foF2 / foE;
        final double ratio = join(fo, 1.75, 20.0, fo - 1.75);

        // deltaM parameter
        double deltaM = -0.012;
        if (foE >= 1e-30) {
            deltaM += 0.253 / (ratio - 1.215);
        }

        // hmF2 Eq. 80
        final double mF2Sq = mF2 * mF2;
        final double temp  = FastMath.sqrt((0.0196 * mF2Sq + 1) / (1.2967 * mF2Sq - 1.0));
        final double height  = ((1490.0 * mF2 * temp) / (mF2 + deltaM)) - 176.0;
        return height;
    }

    /**
     * Computes cf2 coefficients.
     * @param af2 interpolated coefficients for foF2
     * @param t time argument
     * @return the cf2 coefficients array
     */
    private double[] computeCF2(final double[][] af2, final double t) {
        // Eq. 50
        final double[] cf2 = new double[76];
        for (int i = 0; i < cf2.length; i++) {
            double sum = 0.0;
            for (int k = 0; k < 6; k++) {
                final SinCos sc = FastMath.sinCos((k + 1) * t);
                sum += af2[i][2 * k + 1] * sc.sin() + af2[i][2 * (k + 1)] * sc.cos();
            }
            cf2[i] = af2[i][0] + sum;
        }
        return cf2;
    }

    /**
     * Computes Cm3 coefficients.
     * @param am3 interpolated coefficients for foF2
     * @param t time argument
     * @return the Cm3 coefficients array
     */
    private double[] computeCm3(final double[][] am3, final double t) {
        // Eq. 51
        final double[] cm3 = new double[49];
        for (int i = 0; i < cm3.length; i++) {
            double sum = 0.0;
            for (int k = 0; k < 4; k++) {
                final SinCos sc = FastMath.sinCos((k + 1) * t);
                sum += am3[i][2 * k + 1] * sc.sin() + am3[i][2 * (k + 1)] * sc.cos();
            }
            cm3[i] = am3[i][0] + sum;
        }
        return cm3;
    }

    /**
     * This method computes the F2 layer critical frequency.
     * @param modip modified DIP latitude, in degrees
     * @param cf2 Fourier time series for foF2
     * @param latitude latitude in radians
     * @param longitude longitude in radians
     * @return the F2 layer critical frequency, MHz
     */
    private double computefoF2(final double modip, final double[] cf2,
                               final double latitude, final double longitude) {

        // Legendre grades (Eq. 63)
        final int[] q = new int[] {
            12, 12, 9, 5, 2, 1, 1, 1, 1
        };

        // Array for geographic terms
        final double[] g = new double[cf2.length];
        g[0] = 1.0;

        // MODIP coefficients Eq. 57
        final double sinMODIP = FastMath.sin(FastMath.toRadians(modip));
        final double[] m = new double[12];
        m[0] = 1.0;
        for (int i = 1; i < q[0]; i++) {
            m[i] = sinMODIP * m[i - 1];
            g[i] = m[i];
        }

        // Latitude coefficients (Eq. 58)
        final double cosLat = FastMath.cos(latitude);
        final double[] p = new double[8];
        p[0] = cosLat;
        for (int n = 2; n < 9; n++) {
            p[n - 1] = cosLat * p[n - 2];
        }

        // latitude and longitude terms
        int index = 12;
        for (int i = 1; i < q.length; i++) {
            final SinCos sc = FastMath.sinCos(i * longitude);
            for (int j = 0; j < q[i]; j++) {
                g[index++] = m[j] * p[i - 1] * sc.cos();
                g[index++] = m[j] * p[i - 1] * sc.sin();
            }
        }

        // Compute foF2 by linear combination
        final double frequency = MathArrays.linearCombination(cf2, g);
        return frequency;
    }

    /**
     * This method computes the Maximum Usable Frequency factor.
     * @param modip modified DIP latitude, in degrees
     * @param cm3 Fourier time series for M(3000)F2
     * @param latitude latitude in radians
     * @param longitude longitude in radians
     * @return the Maximum Usable Frequency factor
     */
    private double computeMF2(final double modip, final double[] cm3,
                              final double latitude, final double longitude) {

        // Legendre grades (Eq. 71)
        final int[] r = new int[] {
            7, 8, 6, 3, 2, 1, 1
        };

        // Array for geographic terms
        final double[] g = new double[cm3.length];
        g[0] = 1.0;

        // MODIP coefficients Eq. 57
        final double sinMODIP = FastMath.sin(FastMath.toRadians(modip));
        final double[] m = new double[12];
        m[0] = 1.0;
        for (int i = 1; i < 12; i++) {
            m[i] = sinMODIP * m[i - 1];
            if (i < 7) {
                g[i] = m[i];
            }
        }

        // Latitude coefficients (Eq. 58)
        final double cosLat = FastMath.cos(latitude);
        final double[] p = new double[8];
        p[0] = cosLat;
        for (int n = 2; n < 9; n++) {
            p[n - 1] = cosLat * p[n - 2];
        }

        // latitude and longitude terms
        int index = 7;
        for (int i = 1; i < r.length; i++) {
            final SinCos sc = FastMath.sinCos(i * longitude);
            for (int j = 0; j < r[i]; j++) {
                g[index++] = m[j] * p[i - 1] * sc.cos();
                g[index++] = m[j] * p[i - 1] * sc.sin();
            }
        }

        // Compute m3000 by linear combination
        final double m3000 = MathArrays.linearCombination(g, cm3);
        return m3000;
    }

    /**
     * This method computes the F1 layer critical frequency.
     * <p>
     * This computation performs the algorithm exposed in Annex F
     * of the reference document.
     * </p>
     * @param foE the E layer critical frequency, MHz
     * @return the F1 layer critical frequency, MHz
     * @param foF2 the F2 layer critical frequency, MHz
     */
    private double computefoF1(final double foE, final double foF2) {
        final double temp  = join(1.4 * foE, 0.0, 1000.0, foE - 2.0);
        final double temp2 = join(0.0, temp, 1000.0, foE - temp);
        final double value = join(temp2, 0.85 * temp2, 60.0, 0.85 * foF2 - temp2);
        if (value < 1.0E-6) {
            return 0.0;
        } else {
            return value;
        }
    }

    /**
     * This method allows the computation of the F2, F1 and E layer amplitudes.
     * <p>
     * The resulting element is an array having the following form:
     * <ul>
     * <li>double[0] = A1 → F2 layer amplitude
     * <li>double[1] = A2 → F1 layer amplitude
     * <li>double[2] = A3 → E  layer amplitude
     * </ul>
     * </p>
     * @param nmE E layer maximum density in 10^11 m-3
     * @param nmF1 F1 layer maximum density in 10^11 m-3
     * @param foF1 F1 layer critical frequency in MHz
     * @return a three components array containing the layer amplitudes
     */
    private double[] computeLayerAmplitudes(final double nmE, final double nmF1, final double foF1) {
        // Initialize array
        final double[] amplitude = new double[3];

        // F2 layer amplitude (Eq. 90)
        final double a1 = 4.0 * nmF2;
        amplitude[0]   = a1;

        // F1 and E layer amplitudes (Eq. 91 to 98)
        if (foF1 < 0.5) {
            amplitude[1] = 0.0;
            amplitude[2] = 4.0 * (nmE - epst(a1, hmF2, b2Bot, hmE));
        } else {
            double a2a = 0.0;
            double a3a = 4.0 * nmE;
            for (int i = 0; i < 5; i++) {
                a2a = 4.0 * (nmF1 - epst(a1, hmF2, b2Bot, hmF1) - epst(a3a, hmE, beTop, hmF1));
                a2a = join(a2a, 0.8 * nmF1, 1.0, a2a - 0.8 * nmF1);
                a3a = 4.0 * (nmE - epst(a2a, hmF1, b1Bot, hmE) - epst(a1, hmF2, b2Bot, hmE));
            }
            amplitude[1] = a2a;
            amplitude[2] = join(a3a, 0.05, 60.0, a3a - 0.005);
        }

        return amplitude;
    }

    /**
     * This method computes the topside thickness parameter H0.
     *
     * @param month current month
     * @param azr effective sunspot number
     * @return H0 in km
     */
    private double computeH0(final int month, final double azr) {

        // Auxiliary parameter ka (Eq. 99 and 100)
        final double ka;
        if (month > 3 && month < 10) {
            // month = 4,5,6,7,8,9
            ka = 6.705 - 0.014 * azr - 0.008 * hmF2;
        } else {
            // month = 1,2,3,10,11,12
            final double ratio = hmF2 / b2Bot;
            ka = -7.77 + 0.097 * ratio * ratio + 0.153 * nmF2;
        }

        // Auxiliary parameter kb (Eq. 101 and 102)
        double kb = join(ka, 2.0, 1.0, ka - 2.0);
        kb = join(8.0, kb, 1.0, kb - 8.0);

        // Auxiliary parameter Ha (Eq. 103)
        final double hA = kb * b2Bot;

        // Auxiliary parameters x and v (Eq. 104 and 105)
        final double x = 0.01 * (hA - 150.0);
        final double v = (0.041163 * x - 0.183981) * x + 1.424472;

        // Topside thickness parameter (Eq. 106)
        final double h = hA / v;
        return h;
    }

    /**
     * A clipped exponential function.
     * <p>
     * This function, describe in section F.2.12.2 of the reference document, is
     * recommanded for the computation of exponential values.
     * </p>
     * @param power power for exponential function
     * @return clipped exponential value
     */
    private double clipExp(final double power) {
        if (power > 80.0) {
            return 5.5406E34;
        } else if (power < -80) {
            return 1.8049E-35;
        } else {
            return FastMath.exp(power);
        }
    }

    /**
     * This method provides a third order interpolation function
     * as recommended in the reference document (Ref Eq. 128 to Eq. 138)
     *
     * @param z1 z1 coefficient
     * @param z2 z2 coefficient
     * @param z3 z3 coefficient
     * @param z4 z4 coefficient
     * @param x position
     * @return a third order interpolation
     */
    private double interpolate(final double z1, final double z2,
                               final double z3, final double z4,
                               final double x) {

        if (FastMath.abs(2.0 * x) < 1e-10) {
            return z2;
        }

        final double delta = 2.0 * x - 1.0;
        final double g1 = z3 + z2;
        final double g2 = z3 - z2;
        final double g3 = z4 + z1;
        final double g4 = (z4 - z1) / 3.0;
        final double a0 = 9.0 * g1 - g3;
        final double a1 = 9.0 * g2 - g4;
        final double a2 = g3 - g1;
        final double a3 = g4 - g2;
        final double zx = 0.0625 * (a0 + delta * (a1 + delta * (a2 + delta * a3)));

        return zx;
    }

    /**
     * Allows smooth joining of functions f1 and f2
     * (i.e. continuous first derivatives) at origin.
     * <p>
     * This function, describe in section F.2.12.1 of the reference document, is
     * recommanded for computational efficiency.
     * </p>
     * @param dF1 first function
     * @param dF2 second function
     * @param dA width of transition region
     * @param dX x value
     * @return the computed value
     */
    private double join(final double dF1, final double dF2,
                        final double dA, final double dX) {
        final double ee = clipExp(dA * dX);
        return (dF1 * ee + dF2) / (ee + 1.0);
    }

    /**
     * The Epstein function.
     * <p>
     * This function, describe in section 2.5.1 of the reference document, is used
     * as a basis analytical function in NeQuick for the construction of the ionospheric layers.
     * </p>
     * @param x x parameter
     * @param y y parameter
     * @param z z parameter
     * @param w w parameter
     * @return value of the epstein function
     */
    private double epst(final double x, final double y,
                        final double z, final double w) {
        final double ex  = clipExp((w - y) / z);
        final double opex = 1.0 + ex;
        final double epst = x * ex / (opex * opex);
        return epst;
    }

}