/*****************************************************************************\
*                                                                           *
*  The following code was derived from public domain ITM code available     *
*  at ftp://flattop.its.bldrdoc.gov/itm/ITMDLL.cpp that was released on     *
*  June 26, 2007.  It was modified to remove Microsoft Windows "dll-isms",  *
*  redundant and unnecessary #includes, redundant and unnecessary { }'s,    *
*  to initialize uninitialized variables, type cast some variables,         *
*  re-format the code for easier reading, and to replace pow() function     *
*  calls with explicit multiplications wherever possible to increase        *
*  execution speed and improve computational accuracy.                      *
*                                                                           *
*****************************************************************************
*                                                                           *
*  Added comments that refer to itm.pdf and itmalg.pdf in a way to easy     *
*  searching.                                                               *
*                                                                           *
* // :0: Blah, page 0     This is the implementation of the code from       *
*                         itm.pdf, section "0". The description is          *
*                         found on page 0.                                  *
* [Alg 0.0]               please refer to algorithm 0.0 from itm_alg.pdf    *
*                                                                           *
* Holger Schurig, DH3HS                                                     *
*                                                                           *
*****************************************************************************
*                                                                           *
* 2019.07.08 - Fix shadowing of static variable in adiff() function         *
* Ferran Obón Santacana                                                     *
*                                                                           *
\*****************************************************************************/

// SPDX-FileCopyrightText: 2011 Adrian Musceac YO8RZZ
// SPDX-License-Identifier: GPL-2.0-or-later

// *************************************
// C++ routines for this program are taken from
// a translation of the FORTRAN code written by
// U.S. Department of Commerce NTIA/ITS
// Institute for Telecommunication Sciences
// *****************
// Irregular Terrain Model (ITM) (Longley-Rice)
// *************************************



#include <math.h>
#include <complex>
#include <assert.h>
#include <stdio.h>


using namespace std;

namespace ITM {

const double THIRD = (1.0 / 3.0);
const double f_0 = 47.7; // 47.7 MHz from [Alg 1.1], to convert frequency into wavenumber and vica versa


struct prop_type {
    // General input
    double d;        // distance between the two terminals
    double h_g[2];   // antenna structural heights (above ground)
    double k;        // wave number (= radio frequency)
    double delta_h;  // terrain irregularity parameter
    double N_s;      // minimum monthly surface refractivity, n-Units
    double gamma_e;  // earth's effective curvature
    double Z_g_real; // surface transfer impedance of the ground
    double Z_g_imag;
    // Additional input for point-to-point mode
    double h_e[2];      // antenna effective heights
    double d_Lj[2];     // horizon distances
    double theta_ej[2]; // horizontal elevation angles
    int mdp;            // controlling mode: -1: point to point, 1 start of area, 0 area continuation
    // Output
    int kwx;      // error indicator
    double A_ref; // reference attenuation
    // used to be propa_type, computed in lrprop()
    double dx;       // scatter distance
    double ael;      // line-of-sight coefficients
    double ak1;      // dito
    double ak2;      // dito
    double aed;      // diffraction coefficients
    double emd;      // dito
    double aes;      // scatter coefficients
    double ems;      // dito
    double d_Lsj[2]; // smooth earth horizon distances
    double d_Ls;     // d_Lsj[] accumulated
    double d_L;      // d_Lj[] accumulated
    double theta_e;  // theta_ej[] accumulated, total bending angle
};


static int mymin(const int& i, const int& j)
{
    if (i < j)
        return i;
    else
        return j;
}


static int mymax(const int& i, const int& j)
{
    if (i > j)
        return i;
    else
        return j;
}


static double mymin(const double& a, const double& b)
{
    if (a < b)
        return a;
    else
        return b;
}


static double mymax(const double& a, const double& b)
{
    if (a > b)
        return a;
    else
        return b;
}


static double FORTRAN_DIM(const double& x, const double& y)
{
    // This performs the FORTRAN DIM function.
    // result is x-y if x is greater than y; otherwise result is 0.0

    if (x > y)
        return x - y;
    else
        return 0.0;
}

//#define set_warn(txt, err) printf("%s:%d %s\n", __func__, __LINE__, txt);
#define set_warn(txt, err)

// :13: single-knife attenuation, page 6
/*
	* The attenuation due to a single knife edge -- this is an approximation of
	* a Fresnel integral as a function of v^2. The non-approximated integral
	* is documented as [Alg 4.21]
	*
	* Now, what is "single knife attenuation" ? Googling found some paper
	* http://www.its.bldrdoc.gov/pub/ntia-rpt/81-86/81-86.pdf, which actually
	* talks about multi-knife attenuation calculation. However, it says there
	* that single-knife attenuation models attenuation over the edge of one
	* isolated hill.
	*
	* Note that the arguments to this function aren't v, but v^2
	*/
static double Fn(const double& v_square)
{
    double a;

    // [Alg 6.1]
    if (v_square <= 5.76) // this is the 2.40 from the text, but squared
        a = 6.02 + 9.11 * sqrt(v_square) - 1.27 * v_square;
    else
        a = 12.953 + 4.343 * log(v_square);

    return a;
}


// :14: page 6
/*
	* The heigh-gain over a smooth spherical earth -- to be used in the "three
	* radii" mode. The approximation is that given in in [Alg 6.4ff].
	*/
static double F(const double& x, const double& K)
{
    double w, fhtv;
    if (x <= 200.0) {
        // F = F_2(x, L), which is defined in [Alg 6.6]

        w = -log(K);

        // XXX the text says "or x * w^3 > 450"
        if (K < 1e-5 || (x * w * w * w) > 5495.0) {
            // F_2(x, k) = F_1(x), which is defined in [Alg 6.5]
            // XXX but this isn't the same as in itm_alg.pdf
            fhtv = -117.0;
            if (x > 1.0)
                fhtv = 17.372 * log(x) + fhtv;
        } else {
            // [Alg 6.6], lower part
            fhtv = 2.5e-5 * x * x / K - 8.686 * w - 15.0;
        }
    } else {
        // [Alg 6.3] and [Alg 6.4], lower part, which is G(x)
        fhtv = 0.05751 * x - 4.343 * log(x);

        // [Alg 6.4], middle part, but again XXX this doesn't match totally
        if (x < 2000.0) {
            w = 0.0134 * x * exp(-0.005 * x);
            fhtv = (1.0 - w) * fhtv + w * (17.372 * log(x) - 117.0);
        }
    }

    return fhtv;
}


// :25: Tropospheric scatter frequency gain, [Alg 6.10ff], page 12
static double H_0(double r, double et)
{
    // constants from [Alg 6.13]
    const double a[5] = {25.0, 80.0, 177.0, 395.0, 705.0};
    const double b[5] = {24.0, 45.0, 68.0, 80.0, 105.0};
    double q, x;
    int it;
    double h0fv;

    it = (int)et;

    if (it <= 0) {
        it = 1;
        q = 0.0;
    } else if (it >= 4) {
        it = 4;
        q = 0.0;
    } else {
        q = et - it;
    }

    x = (1.0 / r);
    x *= x;
    // [Alg 6.13], calculates something like H_01(r,j), but not really XXX
    h0fv = 4.343 * log((1.0 + a[it - 1] * x + b[it - 1]) * x);

    // XXX not sure what this means
    if (q != 0.0)
        h0fv = (1.0 - q) * h0fv + q * 4.343 * log((a[it] * x + b[it]) * x + 1.0);

    return h0fv;
}


// :25: This is the F(\Theta d) function for scatter fields, page 12
static double F_0(double td)
{
    // [Alg 6.9]
    if (td <= 10e3) // below 10 km
        return 133.4 + 104.6 * td + 71.8 * log(td);
    else if (td <= 70e3) // between 10 km and 70 km
        return 0.332e-3 + 0.212e-3 * td + 0.157e-3 * log(td);
    else // above 70 km
        return -4.343 + -1.086 * td + 2.171 * log(td);
}


// :10: Diffraction attenuation, page 4
static double adiff(double s, prop_type& prop)
{
    /*
		* The function adiff finds the "Diffraction attenuation" at the
		* distance s. It uses a convex combination of smooth earth
		* diffraction and knife-edge diffraction.
		*/
    static double wd1, xd1, A_fo, qk, aht, xht;
    const double A = 151.03;      // dimensionles constant from [Alg 4.20]
    const double D = 50e3;        // 50 km from [Alg 3.9], scale distance for \delta_h(s)
    const double H = 16;          // 16 m  from [Alg 3.10]
    const double ALPHA = 4.77e-4; // from [Alg 4.10]

    if (s == 0.0) {
        complex<double> prop_zgnd(prop.Z_g_real, prop.Z_g_imag);

        // :11: Prepare initial diffraction constants, page 5
        double q = prop.h_g[0] * prop.h_g[1];
        qk = prop.h_e[0] * prop.h_e[1] - q;

        if (prop.mdp < 0.0)
            q += 10.0; // "C" from [Alg 4.9]

        // wd1 and xd1 are parts of Q in [Alg 4.10], but I cannot find this there
        wd1 = sqrt(1.0 + qk / q);
        xd1 = prop.d_L + prop.theta_e / prop.gamma_e; // [Alg 4.9] upper right
        // \delta_h(s), [Alg 3.9]
        q = (1.0 - 0.8 * exp(-prop.d_Ls / D)) * prop.delta_h;
        // \sigma_h(s), [Alg 3.10]
        q *= 0.78 * exp(-pow(q / H, 0.25));

        // A_fo is the "clutter factor"
        A_fo = mymin(15.0, 2.171 * log(1.0 + ALPHA * prop.h_g[0] * prop.h_g[1] * prop.k * q)); // [Alg 4.10]
        // qk is part of the K_j calculation from [Alg 4.17]
        qk = 1.0 / abs(prop_zgnd);
        aht = 20.0; // 20 dB approximation for C_1(K) from [Alg 6.7], see also [Alg 4.25]
        xht = 0.0;

        for (int j = 0; j < 2; ++j) {
            double gamma_j_recip, alpha, K;
            // [Alg 4.15], a is reciproke of gamma_j
            gamma_j_recip = 0.5 * (prop.d_Lj[j] * prop.d_Lj[j]) / prop.h_e[j];
            // [Alg 4.16]
            alpha = pow(gamma_j_recip * prop.k, THIRD);
            // [Alg 4.17]
            K = qk / alpha;
            // [Alg 4.18 and 6.2]
            q = A * (1.607 - K) * alpha * prop.d_Lj[j] / gamma_j_recip;
            // [Alg 4.19, high gain part]
            xht += q;
            // [Alg 4.20] ?,    F(x, k) is in [Alg 6.4]
            aht += F(q, K);
        }
        return 0.0;
    }

    double gamma_o_recip, q, K, ds, theta, alpha, A_r, w;

    // :12: Compute diffraction attenuation, page 5
    // [Alg 4.12]
    theta = prop.theta_e + s * prop.gamma_e;
    // XXX this is not [Alg 4.13]
    ds = s - prop.d_L;
    q = 0.0795775 * prop.k * ds * theta * theta;
    // [Alg 4.14], double knife edge attenuation
    // Note that the arguments to Fn() are not v, but v^2
    double A_k = Fn(q * prop.d_Lj[0] / (ds + prop.d_Lj[0])) +
                 Fn(q * prop.d_Lj[1] / (ds + prop.d_Lj[1]));
    // kinda [Alg 4.15], just so that gamma_o is 1/a
    gamma_o_recip = ds / theta;
    // [Alg 4.16]
    alpha = pow(gamma_o_recip * prop.k, THIRD);
    // [Alg 4.17], note that qk is "1.0 / abs(prop_zgnd)" from above
    K = qk / alpha;
    // [Alg 4.19], q is now X_o
    q = A * (1.607 - K) * alpha * theta + xht;
    // looks a bit like [Alg 4.20], rounded earth attenuation, or??
    // note that G(x) should be "0.05751 * x - 10 * log(q)"
    A_r = 0.05751 * q - 4.343 * log(q) - aht;
    // I'm very unsure if this has anything to do with [Alg 4.9] or not
    q = (wd1 + xd1 / s) * mymin(((1.0 - 0.8 * exp(-s / 50e3)) * prop.delta_h * prop.k), 6283.2);
    // XXX this is NOT the same as the weighting factor from [Alg 4.9]
    w = 25.1 / (25.1 + sqrt(q));
    // [Alg 4.11]
    return (1.0 - w) * A_k + w * A_r + A_fo;
}


/*
	* :22: Scatter attenuation, page 9
	*
	* The function ascat finds the "scatter attenuation" at the distance d. It
	* uses an approximation to the methods of NBS TN101 with check for
	* inadmissable situations. For proper operation, the larger distance (d =
	* d6) must be the first called. A call with d = 0 sets up initial
	* constants.
	*
	* One needs to get TN101, especially chapter 9, to understand this function.
	*/
static double A_scat(double s, prop_type& prop)
{
    static double ad, rr, etq, h0s;

    if (s == 0.0) {
        // :23: Prepare initial scatter constants, page 10
        ad = prop.d_Lj[0] - prop.d_Lj[1];
        rr = prop.h_e[1] / prop.h_e[0];

        if (ad < 0.0) {
            ad = -ad;
            rr = 1.0 / rr;
        }

        etq = (5.67e-6 * prop.N_s - 2.32e-3) * prop.N_s + 0.031; // part of [Alg 4.67]
        h0s = -15.0;
        return 0.0;
    }

    double h0, r1, r2, z0, ss, et, ett, theta_tick, q, temp;

    // :24: Compute scatter attenuation, page 11
    if (h0s > 15.0)
        h0 = h0s;
    else {
        // [Alg 4.61]
        theta_tick = prop.theta_ej[0] + prop.theta_ej[1] + prop.gamma_e * s;
        // [Alg 4.62]
        r2 = 2.0 * prop.k * theta_tick;
        r1 = r2 * prop.h_e[0];
        r2 *= prop.h_e[1];

        if (r1 < 0.2 && r2 < 0.2)
            // The function is undefined
            return 1001.0;

        // XXX not like [Alg 4.65]
        ss = (s - ad) / (s + ad);
        q = rr / ss;
        ss = mymax(0.1, ss);
        q = mymin(mymax(0.1, q), 10.0);
        // XXX not like [Alg 4.66]
        z0 = (s - ad) * (s + ad) * theta_tick * 0.25 / s;
        // [Alg 4.67]
        temp = mymin(1.7, z0 / 8.0e3);
        temp = temp * temp * temp * temp * temp * temp;
        et = (etq * exp(-temp) + 1.0) * z0 / 1.7556e3;

        ett = mymax(et, 1.0);
        h0 = (H_0(r1, ett) + H_0(r2, ett)) * 0.5;                     // [Alg 6.12]
        h0 += mymin(h0, (1.38 - log(ett)) * log(ss) * log(q) * 0.49); // [Alg 6.10 and 6.11]
        h0 = FORTRAN_DIM(h0, 0.0);

        if (et < 1.0)
            // [Alg 6.14]
            h0 = et * h0 + (1.0 - et) * 4.343 * log(pow((1.0 + 1.4142 / r1) * (1.0 + 1.4142 / r2), 2.0) * (r1 + r2) / (r1 + r2 + 2.8284));

        if (h0 > 15.0 && h0s >= 0.0)
            h0 = h0s;
    }

    h0s = h0;
    double theta = prop.theta_e + s * prop.gamma_e; // [Alg 4.60]
    // [Alg 4.63 and 6.8]
    const double D_0 = 40e3; // 40 km from [Alg 6.8]
    const double H = 47.7;   // 47.7 m from [Alg 4.63]
    return 4.343 * log(prop.k * H * theta * theta * theta * theta) + F_0(theta * s) - 0.1 * (prop.N_s - 301.0) * exp(-theta * s / D_0) + h0;
}


static double abq_alos(complex<double> r)
{
    return r.real() * r.real() + r.imag() * r.imag();
}


/*
	* :17: line-of-sight attenuation, page 8
	*
	* The function A_los finds the "line-of-sight attenuation" at the distance
	* d. It uses a convex combination of plane earth fields and diffracted
	* fields. A call with d=0 sets up initial constants.
	*/
static double A_los(double d, prop_type& prop)
{
    static double wls;

    if (d == 0.0) {
        // :18: prepare initial line-of-sight constants, page 8
        const double D1 = 47.7; // 47.7 m from [Alg 4.43]
        const double D2 = 10e3; // 10 km  from [Alg 4.43]
        const double D1R = 1 / D1;
        // weighting factor w
        wls = D1R / (D1R + prop.k * prop.delta_h / mymax(D2, prop.d_Ls)); // [Alg 4.43]
        return 0;
    }

    complex<double> prop_zgnd(prop.Z_g_real, prop.Z_g_imag);
    complex<double> r;
    double s, sps, q;

    // :19: compute line of sight attenuation, page 8
    const double D = 50e3;                        // 50 km from [Alg 3.9]
    const double H = 16.0;                        // 16 m  from [Alg 3.10]
    q = (1.0 - 0.8 * exp(-d / D)) * prop.delta_h; // \Delta h(d), [Alg 3.9]
    s = 0.78 * q * exp(-pow(q / H, 0.25));        // \sigma_h(d), [Alg 3.10]
    q = prop.h_e[0] + prop.h_e[1];
    sps = q / sqrt(d * d + q * q);                                                   // sin(\psi), [Alg 4.46]
    r = (sps - prop_zgnd) / (sps + prop_zgnd) * exp(-mymin(10.0, prop.k * s * sps)); // [Alg 4.47]
    q = abq_alos(r);

    if (q < 0.25 || q < sps) // [Alg 4.48]
        r = r * sqrt(sps / q);

    double alosv = prop.emd * d + prop.aed;           // [Alg 4.45]
    q = prop.k * prop.h_e[0] * prop.h_e[1] * 2.0 / d; // [Alg 4.49]

    // M_PI is pi, M_PI_2 is pi/2
    if (q > M_PI_2) // [Alg 4.50]
        q = M_PI - (M_PI_2 * M_PI_2) / q;

    // [Alg 4.51 and 4.44]
    return (-4.343 * log(abq_alos(complex<double>(cos(q), -sin(q)) + r)) - alosv) * wls + alosv;
}


// :5: LRprop, page 2
/*
	* The value mdp controls some of the program flow. When it equals -1 we are
	* in point-to-point mode, when 1 we are beginning the area mode, and when 0
	* we are continuing the area mode. The assumption is that when one uses the
	* area mode, one will want a sequence of results for varying distances.
	*/
static void lrprop(double d, prop_type& prop)
{
    static bool wlos, wscat;
    static double dmin, xae;
    complex<double> prop_zgnd(prop.Z_g_real, prop.Z_g_imag);
    double a0, a1, a2, a3, a4, a5, a6;
    double d0, d1, d2, d3, d4, d5, d6;
    bool wq;
    double q;
    int j;


    if (prop.mdp != 0) {
        // :6: Do secondary parameters, page 3
        // [Alg 3.5]
        for (j = 0; j < 2; j++)
            prop.d_Lsj[j] = sqrt(2.0 * prop.h_e[j] / prop.gamma_e);

        // [Alg 3.6]
        prop.d_Ls = prop.d_Lsj[0] + prop.d_Lsj[1];
        // [Alg 3.7]
        prop.d_L = prop.d_Lj[0] + prop.d_Lj[1];
        // [Alg 3.8]
        prop.theta_e = mymax(prop.theta_ej[0] + prop.theta_ej[1], -prop.d_L * prop.gamma_e);
        wlos = false;
        wscat = false;

        // :7: Check parameters range, page 3

        // kwx is some kind of error indicator. Setting kwx to 1
        // means "results are slightly bad", while setting it to 4
        // means "my calculations will be bogus"

        // Frequency must be between 40 MHz and 10 GHz
        // Wave number (wn) = 2 * M_PI / lambda, or f/f0, where f0 = 47.7 MHz*m;
        // 0.838 => 40 MHz, 210 => 10 GHz
        if (prop.k < 0.838 || prop.k > 210.0) {
            set_warn("Frequency not optimal", 1);
            prop.kwx = mymax(prop.kwx, 1);
        }

        // Surface refractivity, see [Alg 1.2]
        if (prop.N_s < 250.0 || prop.N_s > 400.0) {
            set_warn("Surface refractivity out-of-bounds", 4);
            prop.kwx = 4;
        } else
            // Earth's effective curvature, see [Alg 1.3]
            if (prop.gamma_e < 75e-9 || prop.gamma_e > 250e-9) {
                set_warn("Earth's curvature out-of-bounds", 4);
                prop.kwx = 4;
            } else
                // Surface transfer impedance, see [Alg 1.4]
                if (prop_zgnd.real() <= abs(prop_zgnd.imag())) {
                    set_warn("Surface transfer impedance out-of-bounds", 4);
                    prop.kwx = 4;
                } else
                    // Calculating outside of 20 MHz to 40 GHz is really bad
                    if (prop.k < 0.419 || prop.k > 420.0) {
                        set_warn("Frequency out-of-bounds", 4);
                        prop.kwx = 4;
                    } else
                        for (j = 0; j < 2; j++) {
                            // Antenna structural height should be between 1 and 1000 m
                            if (prop.h_g[j] < 1.0 || prop.h_g[j] > 1000.0) {
                                set_warn("Antenna height not optimal", 1);
                                prop.kwx = mymax(prop.kwx, 1);
                            }

                            // Horizon elevation angle
                            if (abs(prop.theta_ej[j]) > 200e-3) {
                                set_warn("Horizon elevation weird", 3);
                                prop.kwx = mymax(prop.kwx, 3);
                            }

                            // Horizon distance dl,
                            // Smooth earth horizon distance dls
                            if (prop.d_Lj[j] < 0.1 * prop.d_Lsj[j] ||
                                prop.d_Lj[j] > 3.0 * prop.d_Lsj[j]) {
                                set_warn("Horizon distance weird", 3);
                                prop.kwx = mymax(prop.kwx, 3);
                            }
                            // Antenna structural height must between  0.5 and 3000 m
                            if (prop.h_g[j] < 0.5 || prop.h_g[j] > 3000.0) {
                                set_warn("Antenna heights out-of-bounds", 4);
                                prop.kwx = 4;
                            }
                        }

        dmin = abs(prop.h_e[0] - prop.h_e[1]) / 200e-3;

        // :9: Diffraction coefficients, page 4
        /*
				* This is the region beyond the smooth-earth horizon at
				* D_Lsa and short of where isotropic scatter takes over. It
				* is a key to central region and associated coefficients
				* must always be computed.
				*/
        q = adiff(0.0, prop);
        xae = pow(prop.k * prop.gamma_e * prop.gamma_e, -THIRD); // [Alg 4.2]
        d3 = mymax(prop.d_Ls, 1.3787 * xae + prop.d_L);          // [Alg 4.3]
        d4 = d3 + 2.7574 * xae;                                  // [Alg 4.4]
        a3 = adiff(d3, prop);                                    // [Alg 4.5]
        a4 = adiff(d4, prop);                                    // [Alg 4.7]

        prop.emd = (a4 - a3) / (d4 - d3); // [Alg 4.8]
        prop.aed = a3 - prop.emd * d3;
    }

    if (prop.mdp >= 0) {
        prop.mdp = 0;
        prop.d = d;
    }

    if (prop.d > 0.0) {
        // :8: Check distance, page 3

        // Distance above 1000 km is guessing
        if (prop.d > 1000e3) {
            set_warn("Distance not optimal", 1);
            prop.kwx = mymax(prop.kwx, 1);
        }

        // Distance too small, use some indoor algorithm :-)
        if (prop.d < dmin) {
            set_warn("Distance too small", 3);
            prop.kwx = mymax(prop.kwx, 3);
        }

        // Distance above 2000 km is really bad, don't do that
        if (prop.d < 1e3 || prop.d > 2000e3) {
            set_warn("Distance out-of-bounds", 4);
            prop.kwx = 4;
        }
    }

    if (prop.d < prop.d_Ls) {
        // :15: Line-of-sight calculations, page 7
        if (!wlos) {
            // :16: Line-of-sight coefficients, page 7
            q = A_los(0.0, prop);
            d2 = prop.d_Ls;
            a2 = prop.aed + d2 * prop.emd;
            d0 = 1.908 * prop.k * prop.h_e[0] * prop.h_e[1]; // [Alg 4.38]

            if (prop.aed >= 0.0) {
                d0 = mymin(d0, 0.5 * prop.d_L);   // [Alg 4.28]
                d1 = d0 + 0.25 * (prop.d_L - d0); // [Alg 4.29]
            } else {
                d1 = mymax(-prop.aed / prop.emd, 0.25 * prop.d_L); // [Alg 4.39]
            }

            a1 = A_los(d1, prop); // [Alg 4.31]
            wq = false;

            if (d0 < d1) {
                a0 = A_los(d0, prop); // [Alg 4.30]
                q = log(d2 / d0);
                prop.ak2 = mymax(0.0, ((d2 - d0) * (a1 - a0) - (d1 - d0) * (a2 - a0)) / ((d2 - d0) * log(d1 / d0) - (d1 - d0) * q)); // [Alg 4.32]
                wq = prop.aed >= 0.0 || prop.ak2 > 0.0;

                if (wq) {
                    prop.ak1 = (a2 - a0 - prop.ak2 * q) / (d2 - d0); // [Alg 4.33]

                    if (prop.ak1 < 0.0) {
                        prop.ak1 = 0.0;                     // [Alg 4.36]
                        prop.ak2 = FORTRAN_DIM(a2, a0) / q; // [Alg 4.35]

                        if (prop.ak2 == 0.0) // [Alg 4.37]
                            prop.ak1 = prop.emd;
                    }
                }
            }

            if (!wq) {
                prop.ak1 = FORTRAN_DIM(a2, a1) / (d2 - d1); // [Alg 4.40]
                prop.ak2 = 0.0;                             // [Alg 4.41]

                if (prop.ak1 == 0.0) // [Alg 4.37]
                    prop.ak1 = prop.emd;
            }

            prop.ael = a2 - prop.ak1 * d2 - prop.ak2 * log(d2); // [Alg 4.42]
            wlos = true;
        }

        if (prop.d > 0.0)
            // [Alg 4.1]
            /*
				* The reference attenuation is determined as a function of the distance
				* d from 3 piecewise formulatios. This is part one.
				*/
            prop.A_ref = prop.ael + prop.ak1 * prop.d + prop.ak2 * log(prop.d);
    }

    if (prop.d <= 0.0 || prop.d >= prop.d_Ls) {
        // :20: Troposcatter calculations, page 9
        if (!wscat) {
            // :21: Troposcatter coefficients
            const double DS = 200e3; // 200 km from [Alg 4.53]
            const double HS = 47.7;  // 47.7 m from [Alg 4.59]
            q = A_scat(0.0, prop);
            d5 = prop.d_L + DS;    // [Alg 4.52]
            d6 = d5 + DS;          // [Alg 4.53]
            a6 = A_scat(d6, prop); // [Alg 4.54]
            a5 = A_scat(d5, prop); // [Alg 4.55]

            if (a5 < 1000.0) {
                prop.ems = (a6 - a5) / DS; // [Alg 4.57]
                prop.dx = mymax(prop.d_Ls, // [Alg 4.58]
                                mymax(prop.d_L + 0.3 * xae * log(HS * prop.k),
                                      (a5 - prop.aed - prop.ems * d5) / (prop.emd - prop.ems)));
                prop.aes = (prop.emd - prop.ems) * prop.dx + prop.aed; // [Alg 4.59]
            } else {
                prop.ems = prop.emd;
                prop.aes = prop.aed;
                prop.dx = 10.e6; // [Alg 4.56]
            }
            wscat = true;
        }

        // [Alg 4.1], part two and three.
        if (prop.d > prop.dx)
            prop.A_ref = prop.aes + prop.ems * prop.d; // scatter region
        else
            prop.A_ref = prop.aed + prop.emd * prop.d; // diffraction region
    }

    prop.A_ref = mymax(prop.A_ref, 0.0);
}


/*****************************************************************************/

// :27: Variablility parameters, page 13

struct propv_type {
    // Input:
    int lvar; // control switch for initialisation and/or recalculation
    // 0: no recalculations needed
    // 1: distance changed
    // 2: antenna heights changed
    // 3: frequency changed
    // 4: mdvar changed
    // 5: climate changed, or initialize everything
    int mdvar; // desired mode of variability
    int klim;  // climate indicator
    // Output
    double sgc; // standard deviation of situation variability (confidence)
};


#ifdef WITH_QRLA
// :40: Prepare model for area prediction mode, page 21
/*
			* This is used to prepare the model in the area prediction mode. Normally,
			* one first calls qlrps() and then qlra(). Before calling the latter, one
			* should have defined in prop_type the antenna heights @hg, the terrain
			* irregularity parameter @dh , and (probably through qlrps() the variables
			* @wn, @ens, @gme, and @zgnd. The input @kst will define siting criteria
			* for the terminals, @klimx the climate, and @mdvarx the mode of
			* variability. If @klimx <= 0 or @mdvarx < 0 the associated parameters
			* remain unchanged.
			*/
static void qlra(int kst[], int klimx, int mdvarx, prop_type& prop, propv_type& propv)
{
    double q;

    for (int j = 0; j < 2; ++j) {
        if (kst[j] <= 0)
            // [Alg 3.1]
            prop.h_e[j] = prop.h_g[j];
        else {
            q = 4.0;

            if (kst[j] != 1)
                q = 9.0;

            if (prop.h_g[j] < 5.0)
                q *= sin(0.3141593 * prop.h_g[j]);

            prop.h_e[j] = prop.h_g[j] + (1.0 + q) * exp(-mymin(20.0, 2.0 * prop.h_g[j] / mymax(1e-3, prop.delta_h)));
        }

        // [Alg 3.3], upper function. q is d_Ls_j
        const double H_3 = 5; // 5m from [Alg 3.3]
        q = sqrt(2.0 * prop.h_e[j] / prop.gamma_e);
        prop.d_Lj[j] = q * exp(-0.07 * sqrt(prop.delta_h / mymax(prop.h_e[j], H_3)));
        // [Alg 3.4]
        prop.theta_ej[j] = (0.65 * prop.delta_h * (q / prop.d_Lj[j] - 1.0) - 2.0 * prop.h_e[j]) / q;
    }

    prop.mdp = 1;
    propv.lvar = mymax(propv.lvar, 3);

    if (mdvarx >= 0) {
        propv.mdvar = mdvarx;
        propv.lvar = mymax(propv.lvar, 4);
    }

    if (klimx > 0) {
        propv.klim = klimx;
        propv.lvar = 5;
    }
}
#endif


// :51: Inverse of standard normal complementary probability
/*
			* The approximation is due to C. Hastings, Jr. ("Approximations for digital
			* computers," Princeton Univ. Press, 1955) and the maximum error should be
			* 4.5e-4.
			*/
#if 1
static double qerfi(double q)
{
    double x, t, v;
    const double c0 = 2.515516698;
    const double c1 = 0.802853;
    const double c2 = 0.010328;
    const double d1 = 1.432788;
    const double d2 = 0.189269;
    const double d3 = 0.001308;

    x = 0.5 - q;
    t = mymax(0.5 - fabs(x), 0.000001);
    t = sqrt(-2.0 * log(t));
    v = t - ((c2 * t + c1) * t + c0) / (((d3 * t + d2) * t + d1) * t + 1.0);

    if (x < 0.0)
        v = -v;

    return v;
}
#endif


// :41: preparatory routine, page 20
/*
			* This subroutine converts
			*   frequency @fmhz
			*   surface refractivity reduced to sea level @en0
			*   general system elevation @zsys
			*   polarization and ground constants @eps, @sgm
			* to
			*   wave number @wn,
			*   surface refractivity @ens
			*   effective earth curvature @gme
			*   surface impedance @zgnd
			* It may be used with either the area prediction or the point-to-point mode.
			*/
#if 1
static void qlrps(double fmhz, double zsys, double en0, int ipol, double eps, double sgm, prop_type& prop)

{
    const double Z_1 = 9.46e3; // 9.46 km       from [Alg 1.2]
    const double gma = 157e-9; // 157e-9 1/m    from [Alg 1.3]
    const double N_1 = 179.3;  // 179.3 N-units from [Alg 1.3]
    const double Z_0 = 376.62; // 376.62 Ohm    from [Alg 1.5]

    // Frequency -> Wave number
    prop.k = fmhz / f_0; // [Alg 1.1]

    // Surface refractivity
    prop.N_s = en0;
    if (zsys != 0.0)
        prop.N_s *= exp(-zsys / Z_1); // [Alg 1.2]

    // Earths effective curvature
    prop.gamma_e = gma * (1.0 - 0.04665 * exp(prop.N_s / N_1)); // [Alg 1.3]

    // Surface transfer impedance
    complex<double> zq, prop_zgnd(prop.Z_g_real, prop.Z_g_imag);
    zq = complex<double>(eps, Z_0 * sgm / prop.k); // [Alg 1.5]
    prop_zgnd = sqrt(zq - 1.0);

    // adjust surface transfer impedance for Polarization
    if (ipol != 0.0)
        prop_zgnd = prop_zgnd / zq; // [Alg 1.4]

    prop.Z_g_real = prop_zgnd.real();
    prop.Z_g_imag = prop_zgnd.imag();
}
#endif


// :30: Function curv, page 15
static double curve(double const& c1, double const& c2, double const& x1, double const& x2, double const& x3, double const& de)
{
    double temp1, temp2;

    temp1 = (de - x2) / x3;
    temp2 = de / x1;

    temp1 *= temp1;
    temp2 *= temp2;

    return (c1 + c2 / (1.0 + temp1)) * temp2 / (1.0 + temp2);
}


// :28: Area variablity, page 14
#if 1
static double avar(double zzt, double zzl, double zzc, prop_type& prop, propv_type& propv)
{
    static int kdv;
    static double dexa, de, vmd, vs0, sgl, sgtm, sgtp, sgtd, tgtd, gm, gp, cv1, cv2, yv1, yv2, yv3, csm1, csm2, ysm1, ysm2, ysm3, csp1, csp2, ysp1, ysp2, ysp3, csd1, zd, cfm1, cfm2, cfm3, cfp1, cfp2, cfp3;

    // :29: Climatic constants, page 15
    // Indexes are:
    //   0: equatorial
    //   1: continental suptropical
    //   2: maritime subtropical
    //   3: desert
    //   4: continental temperature
    //   5: maritime over land
    //   6: maritime over sea
    //                      equator  contsup maritsup   desert conttemp mariland  marisea
    const double bv1[7] = {-9.67, -0.62, 1.26, -9.21, -0.62, -0.39, 3.15};
    const double bv2[7] = {12.7, 9.19, 15.5, 9.05, 9.19, 2.86, 857.9};
    const double xv1[7] = {144.9e3, 228.9e3, 262.6e3, 84.1e3, 228.9e3, 141.7e3, 2222.e3};
    const double xv2[7] = {190.3e3, 205.2e3, 185.2e3, 101.1e3, 205.2e3, 315.9e3, 164.8e3};
    const double xv3[7] = {133.8e3, 143.6e3, 99.8e3, 98.6e3, 143.6e3, 167.4e3, 116.3e3};
    const double bsm1[7] = {2.13, 2.66, 6.11, 1.98, 2.68, 6.86, 8.51};
    const double bsm2[7] = {159.5, 7.67, 6.65, 13.11, 7.16, 10.38, 169.8};
    const double xsm1[7] = {762.2e3, 100.4e3, 138.2e3, 139.1e3, 93.7e3, 187.8e3, 609.8e3};
    const double xsm2[7] = {123.6e3, 172.5e3, 242.2e3, 132.7e3, 186.8e3, 169.6e3, 119.9e3};
    const double xsm3[7] = {94.5e3, 136.4e3, 178.6e3, 193.5e3, 133.5e3, 108.9e3, 106.6e3};
    const double bsp1[7] = {2.11, 6.87, 10.08, 3.68, 4.75, 8.58, 8.43};
    const double bsp2[7] = {102.3, 15.53, 9.60, 159.3, 8.12, 13.97, 8.19};
    const double xsp1[7] = {636.9e3, 138.7e3, 165.3e3, 464.4e3, 93.2e3, 216.0e3, 136.2e3};
    const double xsp2[7] = {134.8e3, 143.7e3, 225.7e3, 93.1e3, 135.9e3, 152.0e3, 188.5e3};
    const double xsp3[7] = {95.6e3, 98.6e3, 129.7e3, 94.2e3, 113.4e3, 122.7e3, 122.9e3};
    // bds1 -> is similar to C_D from table 5.1 at [Alg 5.8]
    // bzd1 -> is similar to z_D from table 5.1 at [Alg 5.8]
    const double bsd1[7] = {1.224, 0.801, 1.380, 1.000, 1.224, 1.518, 1.518};
    const double bzd1[7] = {1.282, 2.161, 1.282, 20.0, 1.282, 1.282, 1.282};
    const double bfm1[7] = {1.0, 1.0, 1.0, 1.0, 0.92, 1.0, 1.0};
    const double bfm2[7] = {0.0, 0.0, 0.0, 0.0, 0.25, 0.0, 0.0};
    const double bfm3[7] = {0.0, 0.0, 0.0, 0.0, 1.77, 0.0, 0.0};
    const double bfp1[7] = {1.0, 0.93, 1.0, 0.93, 0.93, 1.0, 1.0};
    const double bfp2[7] = {0.0, 0.31, 0.0, 0.19, 0.31, 0.0, 0.0};
    const double bfp3[7] = {0.0, 2.00, 0.0, 1.79, 2.00, 0.0, 0.0};
    const double rt = 7.8, rl = 24.0;
    static bool no_location_variability, no_situation_variability;
    double avarv, q, vs, zt, zl, zc;
    double sgt, yr;
    int temp_klim;

    if (propv.lvar > 0) {
        // :31: Setup variablity constants, page 16
        switch (propv.lvar) {
        default:
            // Initialization or climate change

            // if climate is wrong, use some "continental temperature" as default
            // and set error indicator
            if (propv.klim <= 0 || propv.klim > 7) {
                propv.klim = 5;
                prop.kwx = mymax(prop.kwx, 2);
                set_warn("Climate index set to continental", 2);
            }

            // convert climate number into index into the climate tables
            temp_klim = propv.klim - 1;

            // :32: Climatic coefficients, page 17
            cv1 = bv1[temp_klim];
            cv2 = bv2[temp_klim];
            yv1 = xv1[temp_klim];
            yv2 = xv2[temp_klim];
            yv3 = xv3[temp_klim];
            csm1 = bsm1[temp_klim];
            csm2 = bsm2[temp_klim];
            ysm1 = xsm1[temp_klim];
            ysm2 = xsm2[temp_klim];
            ysm3 = xsm3[temp_klim];
            csp1 = bsp1[temp_klim];
            csp2 = bsp2[temp_klim];
            ysp1 = xsp1[temp_klim];
            ysp2 = xsp2[temp_klim];
            ysp3 = xsp3[temp_klim];
            csd1 = bsd1[temp_klim];
            zd = bzd1[temp_klim];
            cfm1 = bfm1[temp_klim];
            cfm2 = bfm2[temp_klim];
            cfm3 = bfm3[temp_klim];
            cfp1 = bfp1[temp_klim];
            cfp2 = bfp2[temp_klim];
            cfp3 = bfp3[temp_klim];
            [[fallthrough]];

        case 4:
            // :33: Mode of variablity coefficients, page 17

            // This code means that propv.mdvar can be
            //  0 ..  3
            // 10 .. 13, then no_location_variability is set              (no location variability)
            // 20 .. 23, then no_situation_variability is set             (no situatian variability)
            // 30 .. 33, then no_location_variability and no_situation_variability is set
            kdv = propv.mdvar;
            no_situation_variability = kdv >= 20;
            if (no_situation_variability)
                kdv -= 20;

            no_location_variability = kdv >= 10;
            if (no_location_variability)
                kdv -= 10;

            if (kdv < 0 || kdv > 3) {
                kdv = 0;
                set_warn("kdv set to 0", 2);
                prop.kwx = mymax(prop.kwx, 2);
            }
            [[fallthrough]];

        case 3:
            // :34: Frequency coefficients, page 18
            q = log(0.133 * prop.k);
            gm = cfm1 + cfm2 / ((cfm3 * q * cfm3 * q) + 1.0);
            gp = cfp1 + cfp2 / ((cfp3 * q * cfp3 * q) + 1.0);
            [[fallthrough]];

        case 2:
            // :35: System coefficients, page 18
            // [Alg 5.3], effective distance
            {
                const double a_1 = 9000e3; // 9000 km from [[Alg 5.3]
                //XXX I don't have any idea how they made up the third summand,
                //XXX text says    a_1 * pow(k * D_1, -THIRD)
                //XXX with const double D_1 = 1266; // 1266 km
                dexa = sqrt(2 * a_1 * prop.h_e[0]) +
                       sqrt(2 * a_1 * prop.h_e[1]) +
                       pow((575.7e12 / prop.k), THIRD);
            }
            [[fallthrough]];

        case 1:
            // :36: Distance coefficients, page 18
            {
                // [Alg 5.4]
                const double D_0 = 130e3; // 130 km from [Alg 5.4]
                if (prop.d < dexa)
                    de = D_0 * prop.d / dexa;
                else
                    de = D_0 + prop.d - dexa;
            }
        }
        /*
					* Quantiles of time variability are computed using a
					* variation of the methods described in Section 10 and
					* Annex III.7 of NBS~TN101, and also in CCIR
					* Report {238-3}. Those methods speak of eight or nine
					* discrete radio climates, of which seven have been
					* documented with corresponding empirical curves. It is
					* these empirical curves to which we refer below. They are
					* all curves of quantiles of deviations versus the
					* effective distance @de.
					*/
        vmd = curve(cv1, cv2, yv1, yv2, yv3, de); // [Alg 5.5]
        // [Alg 5.7], the slopes or "pseudo-standard deviations":
        // sgtm -> \sigma T-
        // sgtp -> \sigma T+
        sgtm = curve(csm1, csm2, ysm1, ysm2, ysm3, de) * gm;
        sgtp = curve(csp1, csp2, ysp1, ysp2, ysp3, de) * gp;
        // [Alg 5.8], ducting, "sgtd" -> \sigma TD
        sgtd = sgtp * csd1;
        tgtd = (sgtp - sgtd) * zd;

        // Location variability
        if (no_location_variability) {
            sgl = 0.0;
        } else {
            // Alg [3.9]
            q = (1.0 - 0.8 * exp(-prop.d / 50e3)) * prop.delta_h * prop.k;
            // [Alg 5.9]
            sgl = 10.0 * q / (q + 13.0);
        }

        // Situation variability
        if (no_situation_variability) {
            vs0 = 0.0;
        } else {
            const double D = 100e3;           // 100 km
            vs0 = (5.0 + 3.0 * exp(-de / D)); // [Alg 5.10]
            vs0 *= vs0;
        }

        // Mark all constants as initialized
        propv.lvar = 0;
    }


    // :37: Correct normal deviates, page 19
    zt = zzt;
    zl = zzl;
    zc = zzc;
    // kdv is derived from prop.mdvar
    switch (kdv) {
    case 0:
        // Single message mode. Time, location and situation variability
        // are combined to form a confidence level.
        zt = zc;
        zl = zc;
        break;

    case 1:
        // Individual mode. Reliability is given by time
        // variability. Confidence. is a combination of location
        // and situation variability.
        zl = zc;
        break;

    case 2:
        // Mobile modes. Reliability is a combination of time and
        // location variability. Confidence. is given by the
        // situation variability.
        zl = zt;
        // case 3: Broadcast mode. like avar(zt, zl, zc).
        // Reliability is given by the two-fold statement of at
        // least qT of the time in qL of the locations. Confidence.
        // is given by the situation variability.
    }
    if (fabs(zt) > 3.1 || fabs(zl) > 3.1 || fabs(zc) > 3.1) {
        set_warn("Situations variables not optimal", 1);
        prop.kwx = mymax(prop.kwx, 1);
    }


    // :38: Resolve standard deviations, page 19
    if (zt < 0.0)
        sgt = sgtm;
    else if (zt <= zd)
        sgt = sgtp;
    else
        sgt = sgtd + tgtd / zt;
    // [Alg 5.11], situation variability
    vs = vs0 + (sgt * zt * sgt * zt) / (rt + zc * zc) + (sgl * zl * sgl * zl) / (rl + zc * zc);


    // :39: Resolve deviations, page 19
    if (kdv == 0) {
        yr = 0.0;
        propv.sgc = sqrt(sgt * sgt + sgl * sgl + vs);
    } else if (kdv == 1) {
        yr = sgt * zt;
        propv.sgc = sqrt(sgl * sgl + vs);
    } else if (kdv == 2) {
        yr = sqrt(sgt * sgt + sgl * sgl) * zt;
        propv.sgc = sqrt(vs);
    } else {
        yr = sgt * zt + sgl * zl;
        propv.sgc = sqrt(vs);
    }

    // [Alg 5.1], area variability
    avarv = prop.A_ref - vmd - yr - propv.sgc * zc;

    // [Alg 5.2]
    if (avarv < 0.0)
        avarv = avarv * (29.0 - avarv) / (29.0 - 10.0 * avarv);

    return avarv;
}
#endif

// :45: Find to horizons, page 24
/*
			* Here we use the terrain profile @pfl to find the two horizons. Output
			* consists of the horizon distances @dl and the horizon take-off angles
			* @the. If the path is line-of-sight, the routine sets both horizon
			* distances equal to @dist.
			*/
static void hzns(double pfl[], prop_type& prop)
{
    bool wq;
    int np;
    double xi, za, zb, qc, q, sb, sa;

    np = (int)pfl[0];
    xi = pfl[1];
    za = pfl[2] + prop.h_g[0];
    zb = pfl[np + 2] + prop.h_g[1];
    qc = 0.5 * prop.gamma_e;
    q = qc * prop.d;
    prop.theta_ej[1] = (zb - za) / prop.d;
    prop.theta_ej[0] = prop.theta_ej[1] - q;
    prop.theta_ej[1] = -prop.theta_ej[1] - q;
    prop.d_Lj[0] = prop.d;
    prop.d_Lj[1] = prop.d;

    if (np >= 2) {
        sa = 0.0;
        sb = prop.d;
        wq = true;

        for (int i = 1; i < np; i++) {
            sa += xi;
            sb -= xi;
            q = pfl[i + 2] - (qc * sa + prop.theta_ej[0]) * sa - za;

            if (q > 0.0) {
                prop.theta_ej[0] += q / sa;
                prop.d_Lj[0] = sa;
                wq = false;
            }

            if (!wq) {
                q = pfl[i + 2] - (qc * sb + prop.theta_ej[1]) * sb - zb;

                if (q > 0.0) {
                    prop.theta_ej[1] += q / sb;
                    prop.d_Lj[1] = sb;
                }
            }
        }
    }
}


// :53: Linear least square fit, page 28
static void zlsq1(double z[], const double& x1, const double& x2, double& z0, double& zn)

{
    double xn, xa, xb, x, a, b;
    int n, ja, jb;

    xn = z[0];
    xa = int(FORTRAN_DIM(x1 / z[1], 0.0));
    xb = xn - int(FORTRAN_DIM(xn, x2 / z[1]));

    if (xb <= xa) {
        xa = FORTRAN_DIM(xa, 1.0);
        xb = xn - FORTRAN_DIM(xn, xb + 1.0);
    }

    ja = (int)xa;
    jb = (int)xb;
    n = jb - ja;
    xa = xb - xa;
    x = -0.5 * xa;
    xb += x;
    a = 0.5 * (z[ja + 2] + z[jb + 2]);
    b = 0.5 * (z[ja + 2] - z[jb + 2]) * x;

    for (int i = 2; i <= n; ++i) {
        ++ja;
        x += 1.0;
        a += z[ja + 2];
        b += z[ja + 2] * x;
    }

    a /= xa;
    b = b * 12.0 / ((xa * xa + 2.0) * xa);
    z0 = a - b * xb;
    zn = a + b * (xn - xb);
}


// :52: Provide a quantile and reorders array @a, page 27
static double qtile(const int& nn, double a[], const int& ir)
{
    double q = 0.0, r;                 /* Initializations by KD2BD */
    int m, n, i, j, j1 = 0, i0 = 0, k; /* Initializations by KD2BD */
    bool done = false;
    bool goto10 = true;

    m = 0;
    n = nn;
    k = mymin(mymax(0, ir), n);

    while (!done) {
        if (goto10) {
            q = a[k];
            i0 = m;
            j1 = n;
        }

        i = i0;

        while (i <= n && a[i] >= q)
            i++;

        if (i > n)
            i = n;
        j = j1;

        while (j >= m && a[j] <= q)
            j--;

        if (j < m)
            j = m;

        if (i < j) {
            r = a[i];
            a[i] = a[j];
            a[j] = r;
            i0 = i + 1;
            j1 = j - 1;
            goto10 = false;
        } else if (i < k) {
            a[k] = a[i];
            a[i] = q;
            m = i + 1;
            goto10 = true;
        } else if (j > k) {
            a[k] = a[j];
            a[j] = q;
            n = j - 1;
            goto10 = true;
        } else {
            done = true;
        }
    }

    return q;
}


#ifdef WITH_QERF
// :50: Standard normal complementary probability, page 26
static double qerf(const double& z)
{
    const double b1 = 0.319381530, b2 = -0.356563782, b3 = 1.781477937;
    const double b4 = -1.821255987, b5 = 1.330274429;
    const double rp = 4.317008, rrt2pi = 0.398942280;
    double t, x, qerfv;

    x = z;
    t = fabs(x);

    if (t >= 10.0)
        qerfv = 0.0;
    else {
        t = rp / (t + rp);
        qerfv = exp(-0.5 * x * x) * rrt2pi * ((((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t;
    }

    if (x < 0.0)
        qerfv = 1.0 - qerfv;

    return qerfv;
}
#endif


// :48: Find interdecile range of elevations, page 25
/*
			* Using the terrain profile @pfl we find \Delta h, the interdecile range of
			* elevations between the two points @x1 and @x2.
			*/
static double dlthx(double pfl[], const double& x1, const double& x2)
{
    int np, ka, kb, n, k, j;
    double dlthxv, sn, xa, xb;
    double* s;

    np = (int)pfl[0];
    xa = x1 / pfl[1];
    xb = x2 / pfl[1];
    dlthxv = 0.0;

    if (xb - xa < 2.0) // exit out
        return dlthxv;

    ka = (int)(0.1 * (xb - xa + 8.0));
    ka = mymin(mymax(4, ka), 25);
    n = 10 * ka - 5;
    kb = n - ka + 1;
    sn = n - 1;
    s = new double[n + 2];
    s[0] = sn;
    s[1] = 1.0;
    xb = (xb - xa) / sn;
    k = (int)(xa + 1.0);
    xa -= (double)k;

    for (j = 0; j < n; j++) {
        // Reduce
        while (xa > 0.0 && k < np) {
            xa -= 1.0;
            ++k;
        }

        s[j + 2] = pfl[k + 2] + (pfl[k + 2] - pfl[k + 1]) * xa;
        xa = xa + xb;
    }

    zlsq1(s, 0.0, sn, xa, xb);
    xb = (xb - xa) / sn;

    for (j = 0; j < n; j++) {
        s[j + 2] -= xa;
        xa = xa + xb;
    }

    dlthxv = qtile(n - 1, s + 2, ka - 1) - qtile(n - 1, s + 2, kb - 1);
    dlthxv /= 1.0 - 0.8 * exp(-(x2 - x1) / 50.0e3);
    delete[] s;

    return dlthxv;
}


// :41: Prepare model for point-to-point operation, page 22
/*
			* This mode requires the terrain profile lying between the terminals. This
			* should be a sequence of surface elevations at points along the great
			* circle path joining the two points. It should start at the ground beneath
			* the first terminal and end at the ground beneath the second. In the
			* present routine it is assumed that the elevations are equispaced along
			* the path. They are stored in the array @pfl along with two defining
			* parameters.
			*
			* We will have
			*   pfl[0] = np, the number of points in the path
			*   pfl[1] = xi, the length of each increment
			*/
#if 1
static void qlrpfl(double pfl[], int klimx, int mdvarx, prop_type& prop, propv_type& propv)
{
    int np, j;
    double xl[2], q, za, zb;

    prop.d = pfl[0] * pfl[1];
    np = (int)pfl[0];

    // :44: determine horizons and dh from pfl, page 23
    hzns(pfl, prop);
    for (j = 0; j < 2; j++)
        xl[j] = mymin(15.0 * prop.h_g[j], 0.1 * prop.d_Lj[j]);

    xl[1] = prop.d - xl[1];
    prop.delta_h = dlthx(pfl, xl[0], xl[1]);

    if (prop.d_Lj[0] + prop.d_Lj[1] > 1.5 * prop.d) {
        // :45: Redo line-of-sight horizons, page 23

        /*
					* If the path is line-of-sight, we still need to know where
					* the horizons might have been, and so we turn to
					* techniques used in area prediction mode.
					*/
        zlsq1(pfl, xl[0], xl[1], za, zb);
        prop.h_e[0] = prop.h_g[0] + FORTRAN_DIM(pfl[2], za);
        prop.h_e[1] = prop.h_g[1] + FORTRAN_DIM(pfl[np + 2], zb);

        for (j = 0; j < 2; j++)
            prop.d_Lj[j] = sqrt(2.0 * prop.h_e[j] / prop.gamma_e) * exp(-0.07 * sqrt(prop.delta_h / mymax(prop.h_e[j], 5.0)));

        q = prop.d_Lj[0] + prop.d_Lj[1];

        if (q <= prop.d) {
            q = ((prop.d / q) * (prop.d / q));

            for (j = 0; j < 2; j++) {
                prop.h_e[j] *= q;
                prop.d_Lj[j] = sqrt(2.0 * prop.h_e[j] / prop.gamma_e) * exp(-0.07 * sqrt(prop.delta_h / mymax(prop.h_e[j], 5.0)));
            }
        }

        for (j = 0; j < 2; j++) {
            q = sqrt(2.0 * prop.h_e[j] / prop.gamma_e);
            prop.theta_ej[j] = (0.65 * prop.delta_h * (q / prop.d_Lj[j] - 1.0) - 2.0 * prop.h_e[j]) / q;
        }
    } else {
        // :46: Transhorizon effective heights, page 23
        zlsq1(pfl, xl[0], 0.9 * prop.d_Lj[0], za, q);
        zlsq1(pfl, prop.d - 0.9 * prop.d_Lj[1], xl[1], q, zb);
        prop.h_e[0] = prop.h_g[0] + FORTRAN_DIM(pfl[2], za);
        prop.h_e[1] = prop.h_g[1] + FORTRAN_DIM(pfl[np + 2], zb);
    }

    prop.mdp = -1;
    propv.lvar = mymax(propv.lvar, 3);

    if (mdvarx >= 0) {
        propv.mdvar = mdvarx;
        propv.lvar = mymax(propv.lvar, 4);
    }

    if (klimx > 0) {
        propv.klim = klimx;
        propv.lvar = 5;
    }

    lrprop(0.0, prop);
}
#endif


//********************************************************
//* Point-To-Point Mode Calculations                     *
//********************************************************


#ifdef WITH_POINT_TO_POINT
#include <string.h>
static void point_to_point(double elev[],
                           double tht_m,            // Transceiver above ground level
                           double rht_m,            // Receiver above groud level
                           double eps_dielect,      // Earth dielectric constant (rel. permittivity)
                           double sgm_conductivity, // Earth conductivity (Siemens/m)
                           double eno,              // Atmospheric bending const, n-Units
                           double frq_mhz,
                           int radio_climate, // 1..7
                           int pol,           // 0 horizontal, 1 vertical
                           double conf,       // 0.01 .. .99, Fractions of situations
                           double rel,        // 0.01 .. .99, Fractions of time
                           double& dbloss,
                           char* strmode,
                           int& p_mode,           // propagation mode selector
                           double (&horizons)[2], // horizon distances
                           int& errnum)
{
    // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
    //                4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
    //                7-Maritime Temperate, Over Sea
    // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
    // errnum: 0- No Error.
    //         1- Warning: Some parameters are nearly out of range.
    //                     Results should be used with caution.
    //         2- Note: Default parameters have been substituted for impossible ones.
    //         3- Warning: A combination of parameters is out of range.
    //                     Results are probably invalid.
    //         Other-  Warning: Some parameters are out of range.
    //                          Results are probably invalid.

    prop_type prop;
    propv_type propv;

    double zsys = 0;
    double zc, zr;
    double q = eno;
    long ja, jb, i, np;
    double fs;

    prop.h_g[0] = tht_m; // Tx height above ground level
    prop.h_g[1] = rht_m; // Rx height above ground level
    propv.klim = radio_climate;
    prop.kwx = 0;   // no error yet
    propv.lvar = 5; // initialize all constants
    prop.mdp = -1;  // point-to-point
    zc = qerfi(conf);
    zr = qerfi(rel);
    np = (long)elev[0];
#if 0
					double dkm = (elev[1] * elev[0]) / 1000.0;
					double xkm = elev[1] / 1000.0;
#endif

    ja = (long)(3.0 + 0.1 * elev[0]);
    jb = np - ja + 6;
    for (i = ja - 1; i < jb; ++i)
        zsys += elev[i];
    zsys /= (jb - ja + 1);

    propv.mdvar = 12;
    qlrps(frq_mhz, zsys, q, pol, eps_dielect, sgm_conductivity, prop);
    qlrpfl(elev, propv.klim, propv.mdvar, prop, propv);
    fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
    q = prop.d - prop.d_L;

    horizons[0] = 0.0;
    horizons[1] = 0.0;
    if (int(q) < 0.0) {
        strcpy(strmode, "Line-Of-Sight Mode");
        p_mode = 0;
    } else {
        if (int(q) == 0.0) {
            strcpy(strmode, "Single Horizon");
            horizons[0] = prop.d_Lj[0];
            p_mode = 1;
        }

        else {
            if (int(q) > 0.0) {
                strcpy(strmode, "Double Horizon");
                horizons[0] = prop.d_Lj[0];
                horizons[1] = prop.d_Lj[1];
                p_mode = 1;
            }
        }

        if (prop.d <= prop.d_Ls || prop.d <= prop.dx) {
            strcat(strmode, ", Diffraction Dominant");
            p_mode = 1;
        }

        else {
            if (prop.d > prop.dx) {
                strcat(strmode, ", Troposcatter Dominant");
                p_mode = 2;
            }
        }
    }

    dbloss = avar(zr, 0.0, zc, prop, propv) + fs;
    errnum = prop.kwx;
}
#endif


#ifdef WITH_POINT_TO_POINT_MDH
static void point_to_pointMDH(double elev[], double tht_m, double rht_m, double eps_dielect, double sgm_conductivity, double eno, double frq_mhz, int radio_climate, int pol, double timepct, double locpct, double confpct, double& dbloss, int& propmode, double& deltaH, int& errnum)
{
    // pol: 0-Horizontal, 1-Vertical
    // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
    //                4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
    //                7-Maritime Temperate, Over Sea
    // timepct, locpct, confpct: .01 to .99
    // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
    // propmode:  Value   Mode
    //             -1     mode is undefined
    //              0     Line of Sight
    //              5     Single Horizon, Diffraction
    //              6     Single Horizon, Troposcatter
    //              9     Double Horizon, Diffraction
    //             10     Double Horizon, Troposcatter
    // errnum: 0- No Error.
    //         1- Warning: Some parameters are nearly out of range.
    //                     Results should be used with caution.
    //         2- Note: Default parameters have been substituted for impossible ones.
    //         3- Warning: A combination of parameters is out of range.
    //                     Results are probably invalid.
    //         Other-  Warning: Some parameters are out of range.
    //                          Results are probably invalid.

    prop_type prop;
    propv_type propv;
    propa_type propa;
    double zsys = 0;
    double ztime, zloc, zconf;
    double q = eno;
    long ja, jb, i, np;
    double dkm, xkm;
    double fs;

    propmode = -1; // mode is undefined
    prop.h_g[0] = tht_m;
    prop.h_g[1] = rht_m;
    propv.klim = radio_climate;
    prop.kwx = 0;
    propv.lvar = 5;
    prop.mdp = -1;
    ztime = qerfi(timepct);
    zloc = qerfi(locpct);
    zconf = qerfi(confpct);

    np = (long)elev[0];
    dkm = (elev[1] * elev[0]) / 1000.0;
    xkm = elev[1] / 1000.0;

    ja = (long)(3.0 + 0.1 * elev[0]);
    jb = np - ja + 6;
    for (i = ja - 1; i < jb; ++i)
        zsys += elev[i];
    zsys /= (jb - ja + 1);

    propv.mdvar = 12;
    qlrps(frq_mhz, zsys, q, pol, eps_dielect, sgm_conductivity, prop);
    qlrpfl(elev, propv.klim, propv.mdvar, prop, propa, propv);
    fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
    deltaH = prop.delta_h;
    q = prop.d - prop.d_L;

    if (int(q) < 0.0) {
        propmode = 0; // Line-Of-Sight Mode
    } else {
        if (int(q) == 0.0)
            propmode = 4; // Single Horizon
        else if (int(q) > 0.0)
            propmode = 8; // Double Horizon

        if (prop.d <= prop.d_Ls || prop.d <= prop.dx)
            propmode += 1; // Diffraction Dominant

        else if (prop.d > prop.dx)
            propmode += 2; // Troposcatter Dominant
    }

    dbloss = avar(ztime, zloc, zconf, prop, propv) + fs; //avar(time,location,confidence)
    errnum = prop.kwx;
}
#endif


#ifdef WITH_POINT_TO_POINT_DH
static void point_to_pointDH(double elev[], double tht_m, double rht_m, double eps_dielect, double sgm_conductivity, double eno, double frq_mhz, int radio_climate, int pol, double conf, double rel, double& dbloss, double& deltaH, int& errnum)
{
    // pol: 0-Horizontal, 1-Vertical
    // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
    //                4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
    //                7-Maritime Temperate, Over Sea
    // conf, rel: .01 to .99
    // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
    // errnum: 0- No Error.
    //         1- Warning: Some parameters are nearly out of range.
    //                     Results should be used with caution.
    //         2- Note: Default parameters have been substituted for impossible ones.
    //         3- Warning: A combination of parameters is out of range.
    //                     Results are probably invalid.
    //         Other-  Warning: Some parameters are out of range.
    //                          Results are probably invalid.

    char strmode[100];
    prop_type prop;
    propv_type propv;
    propa_type propa;
    double zsys = 0;
    double zc, zr;
    double q = eno;
    long ja, jb, i, np;
    double dkm, xkm;
    double fs;

    prop.h_g[0] = tht_m;
    prop.h_g[1] = rht_m;
    propv.klim = radio_climate;
    prop.kwx = 0;
    propv.lvar = 5;
    prop.mdp = -1;
    zc = qerfi(conf);
    zr = qerfi(rel);
    np = (long)elev[0];
    dkm = (elev[1] * elev[0]) / 1000.0;
    xkm = elev[1] / 1000.0;

    ja = (long)(3.0 + 0.1 * elev[0]);
    jb = np - ja + 6;
    for (i = ja - 1; i < jb; ++i)
        zsys += elev[i];
    zsys /= (jb - ja + 1);

    propv.mdvar = 12;
    qlrps(frq_mhz, zsys, q, pol, eps_dielect, sgm_conductivity, prop);
    qlrpfl(elev, propv.klim, propv.mdvar, prop, propa, propv);
    fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
    deltaH = prop.delta_h;
    q = prop.d - prop.d_L;

    if (int(q) < 0.0)
        strcpy(strmode, "Line-Of-Sight Mode");
    else {
        if (int(q) == 0.0)
            strcpy(strmode, "Single Horizon");

        else if (int(q) > 0.0)
            strcpy(strmode, "Double Horizon");

        if (prop.d <= prop.d_Ls || prop.d <= prop.dx)
            strcat(strmode, ", Diffraction Dominant");

        else if (prop.d > prop.dx)
            strcat(strmode, ", Troposcatter Dominant");
    }

    dbloss = avar(zr, 0.0, zc, prop, propv) + fs; //avar(time,location,confidence)
    errnum = prop.kwx;
}
#endif


//********************************************************
//* Area Mode Calculations                               *
//********************************************************


#ifdef WITH_AREA
static void area(long ModVar, double deltaH, double tht_m, double rht_m, double dist_km, int TSiteCriteria, int RSiteCriteria, double eps_dielect, double sgm_conductivity, double eno, double frq_mhz, int radio_climate, int pol, double pctTime, double pctLoc, double pctConf, double& dbloss, int& errnum)
{
    // pol: 0-Horizontal, 1-Vertical
    // TSiteCriteria, RSiteCriteria:
    //		   0 - random, 1 - careful, 2 - very careful
    // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
    //                4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
    //                7-Maritime Temperate, Over Sea
    // ModVar: 0 - Single: pctConf is "Time/Situation/Location", pctTime, pctLoc not used
    //         1 - Individual: pctTime is "Situation/Location", pctConf is "Confidence", pctLoc not used
    //         2 - Mobile: pctTime is "Time/Locations (Reliability)", pctConf is "Confidence", pctLoc not used
    //         3 - Broadcast: pctTime is "Time", pctLoc is "Location", pctConf is "Confidence"
    // pctTime, pctLoc, pctConf: .01 to .99
    // errnum: 0- No Error.
    //         1- Warning: Some parameters are nearly out of range.
    //                     Results should be used with caution.
    //         2- Note: Default parameters have been substituted for impossible ones.
    //         3- Warning: A combination of parameters is out of range.
    //                     Results are probably invalid.
    //         Other-  Warning: Some parameters are out of range.
    //                          Results are probably invalid.
    // NOTE: strmode is not used at this time.

    prop_type prop;
    propv_type propv;
    propa_type propa;
    double zt, zl, zc, xlb;
    double fs;
    long ivar;
    double eps, sgm;
    long ipol;
    int kst[2];

    kst[0] = (int)TSiteCriteria;
    kst[1] = (int)RSiteCriteria;
    zt = qerfi(pctTime);
    zl = qerfi(pctLoc);
    zc = qerfi(pctConf);
    eps = eps_dielect;
    sgm = sgm_conductivity;
    prop.delta_h = deltaH;
    prop.h_g[0] = tht_m;
    prop.h_g[1] = rht_m;
    propv.klim = (long)radio_climate;
    prop.N_s = eno;
    prop.kwx = 0;
    ivar = (long)ModVar;
    ipol = (long)pol;
    qlrps(frq_mhz, 0.0, eno, ipol, eps, sgm, prop);
    qlra(kst, propv.klim, ivar, prop, propv);

    if (propv.lvar < 1)
        propv.lvar = 1;

    lrprop(dist_km * 1000.0, prop, propa);
    fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
    xlb = fs + avar(zt, zl, zc, prop, propv);
    dbloss = xlb;

    if (prop.kwx == 0)
        errnum = 0;
    else
        errnum = prop.kwx;
}
#endif

} // namespace ITM
