atcdialog.cxx
trafficcontrol.cxx
CommStation.cxx
- itm.cpp
)
flightgear_component(ATC "${SOURCES}")
+++ /dev/null
-/*****************************************************************************\
- * *
- * 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 *
-\*****************************************************************************/
-
-
-// *************************************
-// 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>
-
-
-const float THIRD = (1.0/3.0);
-const float f_0 = 47.7; // 47.7 MHz from [Alg 1.1], to convert frequency into wavenumber and vica versa
-
-
-using namespace std;
-
-
-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);
-
-
-// :13: single-knife attenuation, page 6
-/*
- * The attenuation due to a sigle 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];
- double 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 chaper 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 alos 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 attentuation, 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
- // Calulating 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.
- */
-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;
-}
-
-
-// :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.
- */
-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]
-
- // Frequecy -> 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();
-}
-
-
-// :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
-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];
- // fall throught
-
- 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)
- no_situation_variability -= 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);
- }
-
- // fall throught
-
- 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);
- // fall throught
-
- 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);
- }
- // fall throught
-
- 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;
-}
-
-
-// :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;
- assert((s = new double[n+2]) != 0);
- 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
-*/
-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);
-}
-
-
-//********************************************************
-//* 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 &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 dkm, xkm;
- 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];
- 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, propv);
- fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
- 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;
- 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
#endif
#include <algorithm>
-#include <math.h>
-#include <stdlib.h>
-#include <deque>
#include <osg/Geode>
#include <osg/Geometry>
#include <Airports/groundnetwork.hxx>
#include <Airports/dynamics.hxx>
#include <Airports/simple.hxx>
-#define WITH_POINT_TO_POINT
-#include "itm.cpp"
using std::sort;
|| (onBoardRadioFreqI1 == stationFreq)) {
if (rec->allowTransmissions()) {
- double snr = calculate_attenuation(rec, parent, ground_to_air);
- if (snr <= 0)
- return;
- if (snr > 0 && snr < 12) {
- //for low SNR values implement a way to make the conversation
- //hard to understand but audible
- //how this works in the real world, is the receiver AGC fails to capture the slope
- //and the signal, due to being amplitude modulated, decreases volume after demodulation
- //the implementation below is more akin to what would happen on a FM transmission
- //therefore the correct way would be to work on the volume
- string hash_noise = " ";
- int reps = fabs((int)snr - 11);
- int t_size = text.size();
- for (int n=1;n<=reps * 2;n++) {
- int pos = rand() % t_size -1;
- text.replace(pos,1, hash_noise);
- }
-
- }
+
fgSetString("/sim/messages/atc", text.c_str());
}
}
}
}
-double FGATCController::calculate_attenuation(FGTrafficRecord * rec, FGAirportDynamics *parent,
- int ground_to_air) {
-
- /// Implement radio attenuation
- /// based on the Longley-Rice propagation model
-
- FGScenery * scenery = globals->get_scenery();
- // player aircraft position
- double own_lat = fgGetDouble("/position/latitude-deg");
- double own_lon = fgGetDouble("/position/longitude-deg");
- double own_alt_ft = fgGetDouble("/position/altitude-ft");
- double own_alt= own_alt_ft * SG_FEET_TO_METER;
-
- //cerr << "ITM:: pilot Lat: " << own_lat << ", Lon: " << own_lon << ", Alt: " << own_alt << endl;
-
- SGGeod own_pos = SGGeod::fromDegM( own_lon, own_lat, own_alt );
- SGGeod max_own_pos = SGGeod::fromDegM( own_lon, own_lat, SG_MAX_ELEVATION_M );
- SGGeoc center = SGGeoc::fromGeod( max_own_pos );
- SGGeoc own_pos_c = SGGeoc::fromGeod( own_pos );
-
- // position of sender radio antenna (HAAT)
- // sender can be aircraft or ground station
- double ATC_HAAT = 30.0;
- double Aircraft_HAAT = 5.0;
- double sender_alt_ft,sender_alt;
- double transmitter_height=0.0;
- double receiver_height=0.0;
- SGGeod sender_pos;
- SGGeod max_sender_pos;
- if(ground_to_air) {
- sender_alt_ft = parent->getElevation();
- sender_alt = sender_alt_ft * SG_FEET_TO_METER;
- sender_pos= SGGeod::fromDegM( parent->getLongitude(),
- parent->getLatitude(), sender_alt );
- max_sender_pos = SGGeod::fromDegM( parent->getLongitude(),
- parent->getLatitude(), SG_MAX_ELEVATION_M );
- }
- else {
- sender_alt_ft = rec->getAltitude();
- sender_alt = sender_alt_ft * SG_FEET_TO_METER;
- sender_pos= SGGeod::fromDegM( rec->getLongitude(),
- rec->getLatitude(), sender_alt );
- max_sender_pos = SGGeod::fromDegM( rec->getLongitude(),
- rec->getLatitude(), SG_MAX_ELEVATION_M );
- }
- SGGeoc sender_pos_c = SGGeoc::fromGeod( sender_pos );
- //cerr << "ITM:: sender Lat: " << parent->getLatitude() << ", Lon: " << parent->getLongitude() << ", Alt: " << sender_alt << endl;
-
- double point_distance= 90.0; // regular SRTM is 90 meters
- double course = SGGeodesy::courseRad(own_pos_c, sender_pos_c);
- double distance_m = SGGeodesy::distanceM(own_pos, sender_pos);
- double probe_distance = 0.0;
- // If distance larger than this value (300 km), assume reception imposssible
- // technically 300 km is no problem if LOS conditions exist,
- // but we do this to spare resources
- if (distance_m > 300000)
- return -1.0;
-
-
- double max_points = distance_m / point_distance;
- deque<double> _elevations;
-
-
- double elevation_under_pilot = 0.0;
- if (scenery->get_elevation_m( max_own_pos, elevation_under_pilot, NULL )) {
- receiver_height = own_alt - elevation_under_pilot + 3; //assume antenna located 3 meters above ground
- }
-
-
- double elevation_under_sender = 0.0;
- if (scenery->get_elevation_m( max_sender_pos, elevation_under_sender, NULL )) {
- transmitter_height = sender_alt - elevation_under_sender;
- }
- else {
- transmitter_height = sender_alt;
- }
- if(ground_to_air)
- transmitter_height += ATC_HAAT;
- else
- transmitter_height += Aircraft_HAAT;
-
- cerr << "ITM:: RX-height: " << receiver_height << ", TX-height: " << transmitter_height << ", Distance: " << distance_m << endl;
-
-
- unsigned int e_size = (deque<unsigned>::size_type)max_points;
-
- while (_elevations.size() <= e_size) {
- probe_distance += point_distance;
- SGGeod probe = SGGeod::fromGeoc(center.advanceRadM( course, probe_distance ));
-
- double elevation_m = 0.0;
-
- if (scenery->get_elevation_m( probe, elevation_m, NULL )) {
- _elevations.push_front(elevation_m);
- //cerr << "ITM:: Probe elev: " << elevation_m << endl;
- }
- else {
- _elevations.push_front(0.0);
- }
- }
- _elevations.push_back(elevation_under_pilot);
- _elevations.push_front(elevation_under_sender);
- double max_alt_between=0.0;
- for( deque<double>::size_type i = 0; i < _elevations.size(); i++ ) {
- if (_elevations[i] > max_alt_between) {
- max_alt_between = _elevations[i];
- }
- }
-
- double num_points= (double)_elevations.size();
- //cerr << "ITM:: Max alt between: " << max_alt_between << ", num points:" << num_points << endl;
- _elevations.push_front(point_distance);
- _elevations.push_front(num_points -1);
- int size = _elevations.size();
- double itm_elev[size];
- for(int i=0;i<size;i++) {
- itm_elev[i]=_elevations[i];
- //cerr << "ITM:: itm_elev: " << _elevations[i] << endl;
- }
-
- ////////////// ITM default parameters //////////////
- // later perhaps take them from tile materials?
- double eps_dielect=15.0;
- double sgm_conductivity = 0.005;
- double eno = 301.0;
- double frq_mhz = 125.0; // middle of bandplan
- int radio_climate = 5; // continental temperate
- int pol=1; // assuming vertical polarization although this is more complex in reality
- double conf = 0.90; // my own tests in Radiomobile have worked best with these values
- double rel = 0.80; // ^^
- double dbloss;
- char strmode[150];
- int errnum;
-
- /////////// radio parameters ///////////
- double receiver_sensitivity = -110.0; // typical AM receiver sensitivity seems to be 0.8 microVolt at 12dB SINAD
- // AM transmitter power in dBm.
- // Note this value is calculated from the typical final transistor stage output
- // !!! small aircraft have portable transmitters which operate at 36 dBm output (4 Watts)
- // later store this value in aircraft description
- // ATC comms usually operate high power equipment, thus making the link asymetrical; this is ignored for now
- double transmitter_power = 43.0;
- double antenna_gain = 2.0; //real-life gain for conventional monopole/dipole antenna
- if(ground_to_air)
- transmitter_power = 49.0;
- else
- transmitter_power = 43.0;
- if(ground_to_air)
- antenna_gain = 5.0; //pilot plane's antenna gain + ground station antenna gain
- else
- antenna_gain = 2.0; //pilot plane's antenna gain + AI aircraft antenna gain
- double link_budget = transmitter_power - receiver_sensitivity + antenna_gain;
-
-
- // first Fresnel zone radius
- // frequency in the middle of the bandplan, more accuracy is not necessary
- double fz_clr= 8.657 * sqrt(distance_m / 0.125);
-
- // TODO: If we clear the first Fresnel zone, we are into line of sight teritory
-
- // else we need to calculate point to point link loss
-
- point_to_point(itm_elev, transmitter_height, receiver_height,
- eps_dielect, sgm_conductivity, eno, frq_mhz, radio_climate,
- pol, conf, rel, dbloss, strmode, errnum);
-
- cerr << "ITM:: Link budget: " << link_budget << ", Attenuation: " << dbloss << " dBm, " << strmode << ", Error: " << errnum << endl;
-
- //if (errnum !=0 && errnum !=1)
- // return -1;
- double snr = link_budget - dbloss;
- return snr;
-
-}
string FGATCController::formatATCFrequency3_2(int freq)
{
string getGateName(FGAIAircraft *aircraft);
virtual void render(bool) = 0;
virtual string getName() = 0;
- double calculate_attenuation(FGTrafficRecord * rec, FGAirportDynamics *parent, int ground_to_air);
private:
include(FlightGearComponent)
set(SOURCES
+ itm.cpp
adf.cxx
agradar.cxx
airspeed_indicator.cxx
--- /dev/null
+// commradio.cxx -- implementation of FGCommRadio
+//
+// Written by Adrian Musceac, started August 2011.
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of the
+// License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+
+
+
+#ifdef HAVE_CONFIG_H
+# include <config.h>
+#endif
+
+#include <math.h>
+#include <stdlib.h>
+#include <deque>
+
+#include <Scenery/scenery.hxx>
+
+#define WITH_POINT_TO_POINT 1
+#include "itm.cpp"
+
+
+FGCommRadio::FGCommRadio(SGPropertyNode *node) {
+
+ /////////// radio parameters ///////////
+ _receiver_sensitivity = -110.0; // typical AM receiver sensitivity seems to be 0.8 microVolt at 12dB SINAD
+ // AM transmitter power in dBm.
+ // Note this value is calculated from the typical final transistor stage output
+ // !!! small aircraft have portable transmitters which operate at 36 dBm output (4 Watts)
+ // later store this value in aircraft description
+ // ATC comms usually operate high power equipment, thus making the link asymetrical; this is ignored for now
+ _transmitter_power = 43.0;
+ //pilot plane's antenna gain + AI aircraft antenna gain
+ //real-life gain for conventional monopole/dipole antenna
+ _antenna_gain = 2.0;
+ _propagation_model = 2; // in the future choose between models via option: realistic radio on/off
+
+}
+
+FGCommRadio::~FGCommRadio()
+{
+}
+
+void FGCommRadio::init ()
+{
+}
+
+
+void FGCommRadio::bind ()
+{
+}
+
+
+void FGCommRadio::update ()
+{
+}
+
+double FGCommRadio::getFrequency(int radio) {
+ double freq = 118.0;
+ switch (radio) {
+ case 1:
+ freq = fgGetDouble("/instrumentation/comm[0]/frequencies/selected-mhz");
+ break;
+ case 2:
+ freq = fgGetDouble("/instrumentation/comm[1]/frequencies/selected-mhz");
+ break;
+ default:
+ freq = fgGetDouble("/instrumentation/comm[0]/frequencies/selected-mhz");
+
+ }
+ return freq;
+}
+
+
+
+
+void FGCommRadio::receive(SGGeod tx_pos, double freq, string text,
+ int ground_to_air) {
+
+ comm1 = getFrequency(1);
+ comm2 = getFrequency(2);
+ if ( (freq != comm1) && (freq != comm2) ) {
+ return;
+ }
+ else {
+ double signal = ITM_calculate_attenuation(tx_pos, freq, ground_to_air);
+ if (signal <= 0)
+ return;
+ if ((signal > 0) && (signal < 12)) {
+ //for low SNR values implement a way to make the conversation
+ //hard to understand but audible
+ //how this works in the real world, is the receiver AGC fails to capture the slope
+ //and the signal, due to being amplitude modulated, decreases volume after demodulation
+ //the implementation below is more akin to what would happen on a FM transmission
+ //therefore the correct way would be to work on the volume
+ string hash_noise = " ";
+ int reps = fabs((int)signal - 11);
+ int t_size = text.size();
+ for (int n=1;n<=reps * 2;n++) {
+ int pos = rand() % t_size -1;
+ text.replace(pos,1, hash_noise);
+ }
+
+ }
+ fgSetString("/sim/messages/atc", text.c_str());
+ }
+
+}
+
+double FGCommRadio::ITM_calculate_attenuation(SGGeod pos, double freq,
+ int ground_to_air) {
+
+ /// Implement radio attenuation
+ /// based on the Longley-Rice propagation model
+
+ ////////////// ITM default parameters //////////////
+ // in the future perhaps take them from tile materials?
+ double eps_dielect=15.0;
+ double sgm_conductivity = 0.005;
+ double eno = 301.0;
+ double frq_mhz;
+ if( (freq < 118.0) || (freq > 137.0) )
+ frq_mhz = 125.0; // sane value, middle of bandplan
+ else
+ frq_mhz = freq;
+ int radio_climate = 5; // continental temperate
+ int pol=1; // assuming vertical polarization although this is more complex in reality
+ double conf = 0.90; // 90% of situations and time, take into account speed
+ double rel = 0.90; // ^^
+ double dbloss;
+ char strmode[150];
+ int errnum;
+
+ double tx_pow,ant_gain;
+ double signal = 0.0;
+
+ if(ground_to_air)
+ tx_pow = _transmitter_power + 6.0;
+
+ if(ground_to_air)
+ ant_gain = _antenna_gain + 3.0; //pilot plane's antenna gain + ground station antenna gain
+
+ double link_budget = tx_pow - _receiver_sensitivity + ant_gain;
+
+ FGScenery * scenery = globals->get_scenery();
+
+ double own_lat = fgGetDouble("/position/latitude-deg");
+ double own_lon = fgGetDouble("/position/longitude-deg");
+ double own_alt_ft = fgGetDouble("/position/altitude-ft");
+ double own_alt= own_alt_ft * SG_FEET_TO_METER;
+
+
+ //cerr << "ITM:: pilot Lat: " << own_lat << ", Lon: " << own_lon << ", Alt: " << own_alt << endl;
+
+ SGGeod own_pos = SGGeod::fromDegM( own_lon, own_lat, own_alt );
+ SGGeod max_own_pos = SGGeod::fromDegM( own_lon, own_lat, SG_MAX_ELEVATION_M );
+ SGGeoc center = SGGeoc::fromGeod( max_own_pos );
+ SGGeoc own_pos_c = SGGeoc::fromGeod( own_pos );
+
+ // position of sender radio antenna (HAAT)
+ // sender can be aircraft or ground station
+ double ATC_HAAT = 30.0;
+ double Aircraft_HAAT = 5.0;
+ double sender_alt_ft,sender_alt;
+ double transmitter_height=0.0;
+ double receiver_height=0.0;
+ SGGeod sender_pos = pos;
+
+ sender_alt_ft = sender_pos.getElevationFt();
+ sender_alt = sender_alt_ft * SG_FEET_TO_METER;
+ SGGeod max_sender_pos = SGGeod::fromGeodM( pos, SG_MAX_ELEVATION_M );
+ SGGeoc sender_pos_c = SGGeoc::fromGeod( sender_pos );
+ //cerr << "ITM:: sender Lat: " << parent->getLatitude() << ", Lon: " << parent->getLongitude() << ", Alt: " << sender_alt << endl;
+
+ double point_distance= 90.0; // regular SRTM is 90 meters
+ double course = SGGeodesy::courseRad(own_pos_c, sender_pos_c);
+ double distance_m = SGGeodesy::distanceM(own_pos, sender_pos);
+ double probe_distance = 0.0;
+ // If distance larger than this value (300 km), assume reception imposssible to preserve resources
+ if (distance_m > 300000)
+ return -1.0;
+ // If above 9000, consider LOS mode and calculate free-space att
+ if (own_alt > 9000) {
+ dbloss = 20 * log10(distance_m) +20 * log10(frq_mhz) -27.55;
+ signal = link_budget - dbloss;
+ return signal;
+ }
+
+
+ double max_points = distance_m / point_distance;
+ deque<double> _elevations;
+
+ double elevation_under_pilot = 0.0;
+ if (scenery->get_elevation_m( max_own_pos, elevation_under_pilot, NULL )) {
+ receiver_height = own_alt - elevation_under_pilot + 3; //assume antenna located 3 meters above ground
+ }
+
+ double elevation_under_sender = 0.0;
+ if (scenery->get_elevation_m( max_sender_pos, elevation_under_sender, NULL )) {
+ transmitter_height = sender_alt - elevation_under_sender;
+ }
+ else {
+ transmitter_height = sender_alt;
+ }
+
+ if(ground_to_air)
+ transmitter_height += ATC_HAAT;
+ else
+ transmitter_height += Aircraft_HAAT;
+
+ cerr << "ITM:: RX-height: " << receiver_height << ", TX-height: " << transmitter_height << ", Distance: " << distance_m << endl;
+
+ unsigned int e_size = (deque<unsigned>::size_type)max_points;
+
+ while (_elevations.size() <= e_size) {
+ probe_distance += point_distance;
+ SGGeod probe = SGGeod::fromGeoc(center.advanceRadM( course, probe_distance ));
+
+ double elevation_m = 0.0;
+
+ if (scenery->get_elevation_m( probe, elevation_m, NULL )) {
+ _elevations.push_front(elevation_m);
+ //cerr << "ITM:: Probe elev: " << elevation_m << endl;
+ }
+ else {
+ _elevations.push_front(0.0);
+ }
+ }
+ _elevations.push_back(elevation_under_pilot);
+ _elevations.push_front(elevation_under_sender);
+
+ double max_alt_between=0.0;
+ for( deque<double>::size_type i = 0; i < _elevations.size(); i++ ) {
+ if (_elevations[i] > max_alt_between) {
+ max_alt_between = _elevations[i];
+ }
+ }
+
+ double num_points= (double)_elevations.size();
+ //cerr << "ITM:: Max alt between: " << max_alt_between << ", num points:" << num_points << endl;
+ _elevations.push_front(point_distance);
+ _elevations.push_front(num_points -1);
+ int size = _elevations.size();
+ double itm_elev[size];
+ for(int i=0;i<size;i++) {
+ itm_elev[i]=_elevations[i];
+ //cerr << "ITM:: itm_elev: " << _elevations[i] << endl;
+ }
+
+
+ // first Fresnel zone radius
+ // frequency in the middle of the bandplan, more accuracy is not necessary
+ double fz_clr= 8.657 * sqrt(distance_m / 0.125);
+
+ // TODO: If we clear the first Fresnel zone, we are into line of sight teritory
+
+ // else we need to calculate point to point link loss
+
+ point_to_point(itm_elev, transmitter_height, receiver_height,
+ eps_dielect, sgm_conductivity, eno, frq_mhz, radio_climate,
+ pol, conf, rel, dbloss, strmode, errnum);
+
+ cerr << "ITM:: Link budget: " << link_budget << ", Attenuation: " << dbloss << " dBm, " << strmode << ", Error: " << errnum << endl;
+
+ //if (errnum == 4)
+ // return -1;
+ signal = link_budget - dbloss;
+ return signal;
+
+}
--- /dev/null
+// commradio.hxx -- class to manage a comm radio instance
+//
+// Written by Adrian Musceac, started August 2011.
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of the
+// License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+
+
+#ifndef __cplusplus
+# error This library requires C++
+#endif
+
+#include <simgear/compiler.h>
+#include <simgear/structure/subsystem_mgr.hxx>
+
+#include <Main/fg_props.hxx>
+
+#include <simgear/math/sg_geodesy.hxx>
+#include <simgear/debug/logstream.hxx>
+#include <string>
+
+
+using std::string;
+
+
+class FGCommRadio : public SGSubsystem, public SGPropertyChangeListener
+{
+private:
+ bool isOperable() const
+ { return _operable; }
+ bool _operable; ///< is the unit serviceable, on, powered, etc
+
+ double _receiver_sensitivity;
+ double _transmitter_power;
+ double _antenna_gain;
+
+ int _propagation_model; /// 0 none, 1 round Earth, 2 ITM
+ double ITM_calculate_attenuation(SGGeod tx_pos, double freq, int ground_to_air)
+
+public:
+
+ FGCommRadio(SGPropertyNode *node);
+ ~FGCommRadio();
+
+ void init ();
+ void bind ();
+ void unbind ();
+ void update (double dt);
+
+ void setFrequency(double freq, int radio);
+ double getFrequency(int radio);
+ void setTxPower(double txpower) { _transmitter_power = txpower; };
+ void receive_text(SGGeod tx_pos, double freq, string text, int ground_to_air);
+ void setPropagationModel(int model) { _propagation_model = model; };
+
+};
+
+
--- /dev/null
+/*****************************************************************************\
+ * *
+ * 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 *
+\*****************************************************************************/
+
+
+// *************************************
+// 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>
+
+
+const float THIRD = (1.0/3.0);
+const float f_0 = 47.7; // 47.7 MHz from [Alg 1.1], to convert frequency into wavenumber and vica versa
+
+
+using namespace std;
+
+
+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);
+
+
+// :13: single-knife attenuation, page 6
+/*
+ * The attenuation due to a sigle 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];
+ double 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 chaper 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 alos 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 attentuation, 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
+ // Calulating 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.
+ */
+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;
+}
+
+
+// :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.
+ */
+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]
+
+ // Frequecy -> 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();
+}
+
+
+// :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
+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];
+ // fall throught
+
+ 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)
+ no_situation_variability -= 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);
+ }
+
+ // fall throught
+
+ 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);
+ // fall throught
+
+ 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);
+ }
+ // fall throught
+
+ 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;
+}
+
+
+// :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;
+ assert((s = new double[n+2]) != 0);
+ 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
+*/
+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);
+}
+
+
+//********************************************************
+//* 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 &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 dkm, xkm;
+ 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];
+ 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, propv);
+ fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
+ 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;
+ 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
projects / flightgear.git / commitdiff