]> git.mxchange.org Git - flightgear.git/blob - src/Radio/itm.cpp
8f75162ce6263f8b7b6aef8d4cedf232a1881eb9
[flightgear.git] / src / Radio / itm.cpp
1 /*****************************************************************************\
2  *                                                                           *
3  *  The following code was derived from public domain ITM code available     *
4  *  at ftp://flattop.its.bldrdoc.gov/itm/ITMDLL.cpp that was released on     *
5  *  June 26, 2007.  It was modified to remove Microsoft Windows "dll-isms",  *
6  *  redundant and unnecessary #includes, redundant and unnecessary { }'s,    *
7  *  to initialize uninitialized variables, type cast some variables,         *
8  *  re-format the code for easier reading, and to replace pow() function     *
9  *  calls with explicit multiplications wherever possible to increase        *
10  *  execution speed and improve computational accuracy.                      *
11  *                                                                           *
12  *****************************************************************************
13  *                                                                           *
14  *  Added comments that refer to itm.pdf and itmalg.pdf in a way to easy     *
15  *  searching.                                                               *
16  *                                                                           *
17  * // :0: Blah, page 0     This is the implementation of the code from       *
18  *                         itm.pdf, section "0". The description is          *
19  *                         found on page 0.                                  *
20  * [Alg 0.0]               please refer to algorithm 0.0 from itm_alg.pdf    *
21  *                                                                           *
22  * Holger Schurig, DH3HS                                                     *
23 \*****************************************************************************/
24
25
26 // *************************************
27 // C++ routines for this program are taken from
28 // a translation of the FORTRAN code written by
29 // U.S. Department of Commerce NTIA/ITS
30 // Institute for Telecommunication Sciences
31 // *****************
32 // Irregular Terrain Model (ITM) (Longley-Rice)
33 // *************************************
34
35
36
37 #include <math.h>
38 #include <complex>
39 #include <assert.h>
40 #include <stdio.h>
41
42
43 using namespace std;
44
45 namespace ITM {
46         
47 const double THIRD = (1.0/3.0);
48 const double f_0 = 47.7; // 47.7 MHz from [Alg 1.1], to convert frequency into wavenumber and vica versa
49
50
51 struct prop_type {
52         // General input
53         double d;          // distance between the two terminals
54         double h_g[2];     // antenna structural heights (above ground)
55         double k;          // wave number (= radio frequency)
56         double delta_h;    // terrain irregularity parameter
57         double N_s;        // minimum monthly surface refractivity, n-Units
58         double gamma_e;    // earth's effective curvature
59         double Z_g_real;   // surface transfer impedance of the ground
60         double Z_g_imag;
61         // Additional input for point-to-point mode
62         double h_e[2];     // antenna effective heights
63         double d_Lj[2];    // horizon distances
64         double theta_ej[2];// horizontal elevation angles
65         int mdp;           // controlling mode: -1: point to point, 1 start of area, 0 area continuation
66         // Output
67         int kwx;           // error indicator
68         double A_ref;      // reference attenuation
69         // used to be propa_type, computed in lrprop()
70         double dx;         // scatter distance
71         double ael;        // line-of-sight coefficients
72         double ak1;        // dito
73         double ak2;        // dito
74         double aed;        // diffraction coefficients
75         double emd;        // dito
76         double aes;        // scatter coefficients
77         double ems;        // dito
78         double d_Lsj[2];   // smooth earth horizon distances
79         double d_Ls;       // d_Lsj[] accumulated
80         double d_L;        // d_Lj[] accumulated
81         double theta_e;    // theta_ej[] accumulated, total bending angle
82 };
83
84
85 static
86 int mymin(const int &i, const int &j)
87 {
88         if (i < j)
89                 return i;
90         else
91                 return j;
92 }
93
94
95 static
96 int mymax(const int &i, const int &j)
97 {
98         if (i > j)
99                 return i;
100         else
101                 return j;
102 }
103
104
105 static
106 double mymin(const double &a, const double &b)
107 {
108         if (a < b)
109                 return a;
110         else
111                 return b;
112 }
113
114
115 static
116 double mymax(const double &a, const double &b)
117 {
118         if (a > b)
119                 return a;
120         else
121                 return b;
122 }
123
124
125 static
126 double FORTRAN_DIM(const double &x, const double &y)
127 {
128         // This performs the FORTRAN DIM function.
129         // result is x-y if x is greater than y; otherwise result is 0.0
130
131         if (x > y)
132                 return x - y;
133         else
134                 return 0.0;
135 }
136
137 //#define set_warn(txt, err) printf("%s:%d %s\n", __func__, __LINE__, txt);
138 #define set_warn(txt, err)
139
140 // :13: single-knife attenuation, page 6
141 /*
142  * The attenuation due to a sigle knife edge -- this is an approximation of
143  * a Fresnel integral as a function of v^2. The non-approximated integral
144  * is documented as [Alg 4.21]
145  *
146  * Now, what is "single knife attenuation" ? Googling found some paper
147  * http://www.its.bldrdoc.gov/pub/ntia-rpt/81-86/81-86.pdf, which actually
148  * talks about multi-knife attenuation calculation. However, it says there
149  * that single-knife attenuation models attenuation over the edge of one
150  * isolated hill.
151  *
152  * Note that the arguments to this function aren't v, but v^2
153  */
154 static
155 double Fn(const double &v_square)
156 {
157         double a;
158
159         // [Alg 6.1]
160         if (v_square <= 5.76)  // this is the 2.40 from the text, but squared
161                 a = 6.02 + 9.11 * sqrt(v_square) - 1.27 * v_square;
162         else
163                 a = 12.953 + 4.343 * log(v_square);
164
165         return a;
166 }
167
168
169 // :14: page 6
170 /*
171  * The heigh-gain over a smooth spherical earth -- to be used in the "three
172  * radii" mode. The approximation is that given in in [Alg 6.4ff].
173  */
174 static
175 double F(const double& x, const double& K)
176 {
177         double w, fhtv;
178         if (x <= 200.0) {
179                 // F = F_2(x, L), which is defined in [Alg 6.6]
180
181                 w = -log(K);
182
183                 // XXX the text says "or x * w^3 > 450"
184                 if (K < 1e-5 || (x * w * w * w) > 5495.0) {
185                         // F_2(x, k) = F_1(x), which is defined in [Alg 6.5]
186                         // XXX but this isn't the same as in itm_alg.pdf
187                         fhtv = -117.0;
188                         if (x > 1.0)
189                                 fhtv = 17.372 * log(x) + fhtv;
190                 } else {
191                         // [Alg 6.6], lower part
192                         fhtv = 2.5e-5 * x * x / K - 8.686 * w - 15.0;
193                 }
194         } else {
195                 // [Alg 6.3] and [Alg 6.4], lower part, which is G(x)
196                 fhtv = 0.05751 * x - 4.343 * log(x);
197
198                 // [Alg 6.4], middle part, but again XXX this doesn't match totally
199                 if (x < 2000.0) {
200                         w = 0.0134 * x * exp(-0.005 * x);
201                         fhtv = (1.0 - w) * fhtv + w * (17.372 * log(x) - 117.0);
202                 }
203         }
204
205         return fhtv;
206 }
207
208
209 // :25: Tropospheric scatter frequency gain, [Alg 6.10ff], page 12
210 static
211 double H_0(double r, double et)
212 {
213         // constants from [Alg 6.13]
214         const double a[5] = {25.0, 80.0, 177.0, 395.0, 705.0};
215         const double b[5] = {24.0, 45.0,  68.0,  80.0, 105.0};
216         double q, x;
217         int it;
218         double h0fv;
219
220         it = (int)et;
221
222         if (it <= 0) {
223                 it = 1;
224                 q = 0.0;
225         } else
226         if (it >= 4) {
227                 it = 4;
228                 q = 0.0;
229         } else {
230                 q = et - it;
231         }
232
233         x = (1.0 / r);
234         x *= x;
235         // [Alg 6.13], calculates something like H_01(r,j), but not really XXX
236         h0fv = 4.343 * log((1.0 + a[it-1] * x + b[it-1]) * x);
237
238         // XXX not sure what this means
239         if (q != 0.0)
240                 h0fv = (1.0 - q) * h0fv + q * 4.343 * log((a[it] * x + b[it]) * x + 1.0);
241
242         return h0fv;
243 }
244
245
246 // :25: This is the F(\Theta d) function for scatter fields, page 12
247 static
248 double F_0(double td)
249 {
250         // [Alg 6.9]
251         if (td <= 10e3) // below 10 km
252                 return 133.4    + 104.6    * td + 71.8     * log(td);
253         else
254         if (td <= 70e3) // between 10 km and 70 km
255                 return 0.332e-3 + 0.212e-3 * td + 0.157e-3 * log(td);
256         else // above 70 km
257                 return-4.343    + -1.086   * td + 2.171    * log(td);
258 }
259
260
261 // :10: Diffraction attenuation, page 4
262 static
263 double  adiff(double s, prop_type &prop)
264 {
265         /*
266          * The function adiff finds the "Diffraction attenuation" at the
267          * distance s. It uses a convex combination of smooth earth
268          * diffraction and knife-edge diffraction.
269          */
270         static double wd1, xd1, A_fo, qk, aht, xht;
271         const double A = 151.03;      // dimensionles constant from [Alg 4.20]
272         const double D = 50e3;        // 50 km from [Alg 3.9], scale distance for \delta_h(s)
273         const double H = 16;          // 16 m  from [Alg 3.10]
274         const double ALPHA = 4.77e-4; // from [Alg 4.10]
275
276         if (s == 0.0) {
277                 complex<double> prop_zgnd(prop.Z_g_real, prop.Z_g_imag);
278
279                 // :11: Prepare initial diffraction constants, page 5
280                 double q = prop.h_g[0] * prop.h_g[1];
281                 double qk = prop.h_e[0] * prop.h_e[1] - q;
282
283                 if (prop.mdp < 0.0)
284                         q += 10.0; // "C" from [Alg 4.9]
285
286                 // wd1 and xd1 are parts of Q in [Alg 4.10], but I cannot find this there
287                 wd1 = sqrt(1.0 + qk / q);
288                 xd1 = prop.d_L + prop.theta_e / prop.gamma_e;  // [Alg 4.9] upper right
289                 // \delta_h(s), [Alg 3.9]
290                 q = (1.0 - 0.8 * exp(-prop.d_Ls / D)) * prop.delta_h;
291                 // \sigma_h(s), [Alg 3.10]
292                 q *= 0.78 * exp(-pow(q / H, 0.25));
293
294                 // A_fo is the "clutter factor"
295                 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]
296                 // qk is part of the K_j calculation from [Alg 4.17]
297                 qk = 1.0 / abs(prop_zgnd);
298                 aht = 20.0; // 20 dB approximation for C_1(K) from [Alg 6.7], see also [Alg 4.25]
299                 xht = 0.0;
300
301                 for (int j = 0; j < 2; ++j) {
302                         double gamma_j_recip, alpha, K;
303                         // [Alg 4.15], a is reciproke of gamma_j
304                         gamma_j_recip = 0.5 * (prop.d_Lj[j] * prop.d_Lj[j]) / prop.h_e[j];
305                         // [Alg 4.16]
306                         alpha = pow(gamma_j_recip * prop.k, THIRD);
307                         // [Alg 4.17]
308                         K = qk / alpha;
309                         // [Alg 4.18 and 6.2]
310                         q = A * (1.607 - K) * alpha * prop.d_Lj[j] / gamma_j_recip;
311                         // [Alg 4.19, high gain part]
312                         xht += q;
313                         // [Alg 4.20] ?,    F(x, k) is in [Alg 6.4]
314                         aht += F(q, K);
315                 }
316                 return 0.0;
317         }
318
319         double gamma_o_recip, q, K, ds, theta, alpha, A_r, w;
320
321         // :12: Compute diffraction attenuation, page 5
322         // [Alg 4.12]
323         theta = prop.theta_e + s * prop.gamma_e;
324         // XXX this is not [Alg 4.13]
325         ds = s - prop.d_L;
326         q = 0.0795775 * prop.k * ds * theta * theta;
327         // [Alg 4.14], double knife edge attenuation
328         // Note that the arguments to Fn() are not v, but v^2
329         double A_k = Fn(q * prop.d_Lj[0] / (ds + prop.d_Lj[0])) +
330                      Fn(q * prop.d_Lj[1] / (ds + prop.d_Lj[1]));
331         // kinda [Alg 4.15], just so that gamma_o is 1/a
332         gamma_o_recip = ds / theta;
333         // [Alg 4.16]
334         alpha = pow(gamma_o_recip * prop.k, THIRD);
335         // [Alg 4.17], note that qk is "1.0 / abs(prop_zgnd)" from above
336         K = qk / alpha;
337         // [Alg 4.19], q is now X_o
338         q = A * (1.607 - K) * alpha * theta + xht;
339         // looks a bit like [Alg 4.20], rounded earth attenuation, or??
340         // note that G(x) should be "0.05751 * x - 10 * log(q)"
341         A_r = 0.05751 * q - 4.343 * log(q) - aht;
342         // I'm very unsure if this has anything to do with [Alg 4.9] or not
343         q = (wd1 + xd1 / s) * mymin(((1.0 - 0.8 * exp(-s / 50e3)) * prop.delta_h * prop.k), 6283.2);
344         // XXX this is NOT the same as the weighting factor from [Alg 4.9]
345         w = 25.1 / (25.1 + sqrt(q));
346         // [Alg 4.11]
347         return (1.0 - w) * A_k + w * A_r + A_fo;
348 }
349
350
351 /*
352  * :22: Scatter attenuation, page 9
353  *
354  * The function ascat finds the "scatter attenuation" at the distance d. It
355  * uses an approximation to the methods of NBS TN101 with check for
356  * inadmissable situations. For proper operation, the larger distance (d =
357  * d6) must be the first called. A call with d = 0 sets up initial
358  * constants.
359  *
360  * One needs to get TN101, especially chaper 9, to understand this function.
361  */
362 static
363 double A_scat(double s, prop_type &prop)
364 {
365         static double ad, rr, etq, h0s;
366
367         if (s == 0.0) {
368                 // :23: Prepare initial scatter constants, page 10
369                 ad = prop.d_Lj[0] - prop.d_Lj[1];
370                 rr = prop.h_e[1] / prop.h_e[0];
371
372                 if (ad < 0.0) {
373                         ad = -ad;
374                         rr = 1.0 / rr;
375                 }
376
377                 etq = (5.67e-6 * prop.N_s - 2.32e-3) * prop.N_s + 0.031; // part of [Alg 4.67]
378                 h0s = -15.0;
379                 return 0.0;
380         }
381
382         double h0, r1, r2, z0, ss, et, ett, theta_tick, q, temp;
383
384         // :24: Compute scatter attenuation, page 11
385         if (h0s > 15.0)
386                 h0 = h0s;
387         else {
388                 // [Alg 4.61]
389                 theta_tick = prop.theta_ej[0] + prop.theta_ej[1] + prop.gamma_e * s;
390                 // [Alg 4.62]
391                 r2 = 2.0 * prop.k * theta_tick;
392                 r1 = r2 * prop.h_e[0];
393                 r2 *= prop.h_e[1];
394
395                 if (r1 < 0.2 && r2 < 0.2)
396                         // The function is undefined
397                         return 1001.0;
398
399                 // XXX not like [Alg 4.65]
400                 ss = (s - ad) / (s + ad);
401                 q = rr / ss;
402                 ss = mymax(0.1, ss);
403                 q = mymin(mymax(0.1, q), 10.0);
404                 // XXX not like [Alg 4.66]
405                 z0 = (s - ad) * (s + ad) * theta_tick * 0.25 / s;
406                 // [Alg 4.67]
407                 temp = mymin(1.7, z0 / 8.0e3);
408                 temp = temp * temp * temp * temp * temp * temp;
409                 et = (etq * exp(-temp) + 1.0) * z0 / 1.7556e3;
410
411                 ett = mymax(et, 1.0);
412                 h0 = (H_0(r1, ett) + H_0(r2, ett)) * 0.5;  // [Alg 6.12]
413                 h0 += mymin(h0, (1.38 - log(ett)) * log(ss) * log(q) * 0.49);  // [Alg 6.10 and 6.11]
414                 h0 = FORTRAN_DIM(h0, 0.0);
415
416                 if (et < 1.0)
417                         // [Alg 6.14]
418                         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));
419
420                 if (h0 > 15.0 && h0s >= 0.0)
421                         h0 = h0s;
422         }
423
424         h0s = h0;
425         double theta = prop.theta_e + s * prop.gamma_e;  // [Alg 4.60]
426         // [Alg 4.63 and 6.8]
427         const double D_0 = 40e3; // 40 km from [Alg 6.8]
428         const double H = 47.7;   // 47.7 m from [Alg 4.63]
429         return 4.343 * log(prop.k * H * theta * theta * theta * theta)
430                + F_0(theta * s)
431                - 0.1 * (prop.N_s - 301.0) * exp(-theta * s / D_0)
432                + h0;
433 }
434
435
436 static
437 double abq_alos(complex<double> r)
438 {
439         return r.real() * r.real() + r.imag() * r.imag();
440 }
441
442
443 /*
444  * :17: line-of-sight attenuation, page 8
445  *
446  * The function alos finds the "line-of-sight attenuation" at the distance
447  * d. It uses a convex combination of plane earth fields and diffracted
448  * fields. A call with d=0 sets up initial constants.
449  */
450 static
451 double A_los(double d, prop_type &prop)
452 {
453         static double wls;
454
455         if (d == 0.0) {
456                 // :18: prepare initial line-of-sight constants, page 8
457                 const double D1 = 47.7; // 47.7 m from [Alg 4.43]
458                 const double D2 = 10e3; // 10 km  from [Alg 4.43]
459                 const double D1R = 1 / D1;
460                 // weighting factor w
461                 wls = D1R / (D1R + prop.k * prop.delta_h / mymax(D2, prop.d_Ls)); // [Alg 4.43]
462                 return 0;
463         }
464
465         complex<double> prop_zgnd(prop.Z_g_real, prop.Z_g_imag);
466         complex<double> r;
467         double s, sps, q;
468
469         // :19: compute line of sight attentuation, page 8
470         const double D = 50e3; // 50 km from [Alg 3.9]
471         const double H = 16.0; // 16 m  from [Alg 3.10]
472         q = (1.0 - 0.8 * exp(-d / D)) * prop.delta_h;     // \Delta h(d), [Alg 3.9]
473         s = 0.78 * q * exp(-pow(q / H, 0.25));       // \sigma_h(d), [Alg 3.10]
474         q = prop.h_e[0] + prop.h_e[1];
475         sps = q / sqrt(d * d + q * q);               // sin(\psi), [Alg 4.46]
476         r = (sps - prop_zgnd) / (sps + prop_zgnd) * exp(-mymin(10.0, prop.k * s * sps)); // [Alg 4.47]
477         q = abq_alos(r);
478
479         if (q < 0.25 || q < sps)                     // [Alg 4.48]
480                 r = r * sqrt(sps / q);
481
482         double alosv = prop.emd * d + prop.aed;    // [Alg 4.45]
483         q = prop.k * prop.h_e[0] * prop.h_e[1] * 2.0 / d; // [Alg 4.49]
484
485         // M_PI is pi, M_PI_2 is pi/2
486         if (q > M_PI_2)                              // [Alg 4.50]
487                 q = M_PI - (M_PI_2 * M_PI_2) / q;
488
489         // [Alg 4.51 and 4.44]
490         return (-4.343 * log(abq_alos(complex<double>(cos(q), -sin(q)) + r)) - alosv) * wls + alosv;
491 }
492
493
494 // :5: LRprop, page 2
495 /*
496  * The value mdp controls some of the program flow. When it equals -1 we are
497  * in point-to-point mode, when 1 we are beginning the area mode, and when 0
498  * we are continuing the area mode. The assumption is that when one uses the
499  * area mode, one will want a sequence of results for varying distances.
500  */
501 static
502 void lrprop(double d, prop_type &prop)
503 {
504         static bool wlos, wscat;
505         static double dmin, xae;
506         complex<double> prop_zgnd(prop.Z_g_real, prop.Z_g_imag);
507         double a0, a1, a2, a3, a4, a5, a6;
508         double d0, d1, d2, d3, d4, d5, d6;
509         bool wq;
510         double q;
511         int j;
512
513
514         if (prop.mdp != 0) {
515                 // :6: Do secondary parameters, page 3
516                 // [Alg 3.5]
517                 for (j = 0; j < 2; j++)
518                         prop.d_Lsj[j] = sqrt(2.0 * prop.h_e[j] / prop.gamma_e);
519
520                 // [Alg 3.6]
521                 prop.d_Ls = prop.d_Lsj[0] + prop.d_Lsj[1];
522                 // [Alg 3.7]
523                 prop.d_L = prop.d_Lj[0] + prop.d_Lj[1];
524                 // [Alg 3.8]
525                 prop.theta_e = mymax(prop.theta_ej[0] + prop.theta_ej[1], -prop.d_L * prop.gamma_e);
526                 wlos = false;
527                 wscat = false;
528
529                 // :7: Check parameters range, page 3
530
531                 // kwx is some kind of error indicator. Setting kwx to 1
532                 // means "results are slightly bad", while setting it to 4
533                 // means "my calculations will be bogus"
534
535                 // Frequency must be between 40 MHz and 10 GHz
536                 // Wave number (wn) = 2 * M_PI / lambda, or f/f0, where f0 = 47.7 MHz*m;
537                 // 0.838 => 40 MHz, 210 => 10 GHz
538                 if (prop.k < 0.838 || prop.k > 210.0) {
539                         set_warn("Frequency not optimal", 1);
540                         prop.kwx = mymax(prop.kwx, 1);
541                 }
542
543                 // Surface refractivity, see [Alg 1.2]
544                 if (prop.N_s < 250.0 || prop.N_s > 400.0) {
545                         set_warn("Surface refractivity out-of-bounds", 4);
546                         prop.kwx = 4;
547                 } else
548                 // Earth's effective curvature, see [Alg 1.3]
549                 if (prop.gamma_e < 75e-9 || prop.gamma_e > 250e-9) {
550                         set_warn("Earth's curvature out-of-bounds", 4);
551                         prop.kwx = 4;
552                 } else
553                 // Surface transfer impedance, see [Alg 1.4]
554                 if (prop_zgnd.real() <= abs(prop_zgnd.imag())) {
555                         set_warn("Surface transfer impedance out-of-bounds", 4);
556                         prop.kwx = 4;
557                 } else
558                 // Calulating outside of 20 MHz to 40 GHz is really bad
559                 if (prop.k < 0.419 || prop.k > 420.0) {
560                         set_warn("Frequency out-of-bounds", 4);
561                         prop.kwx = 4;
562                 } else
563                 for (j = 0; j < 2; j++) {
564                         // Antenna structural height should be between 1 and 1000 m
565                         if (prop.h_g[j] < 1.0 || prop.h_g[j] > 1000.0) {
566                                 set_warn("Antenna height not optimal", 1);
567                                 prop.kwx = mymax(prop.kwx, 1);
568                         }
569
570                         // Horizon elevation angle
571                         if (abs(prop.theta_ej[j]) > 200e-3) {
572                                 set_warn("Horizon elevation weird", 3);
573                                 prop.kwx = mymax(prop.kwx, 3);
574                         }
575
576                         // Horizon distance dl,
577                         // Smooth earth horizon distance dls
578                         if (prop.d_Lj[j] < 0.1 * prop.d_Lsj[j] ||
579                             prop.d_Lj[j] > 3.0 * prop.d_Lsj[j]) {
580                                 set_warn("Horizon distance weird", 3);
581                                 prop.kwx = mymax(prop.kwx, 3);
582                         }
583                         // Antenna structural height must between  0.5 and 3000 m
584                         if (prop.h_g[j] < 0.5 || prop.h_g[j] > 3000.0) {
585                                 set_warn("Antenna heights out-of-bounds", 4);
586                                 prop.kwx = 4;
587                         }
588                 }
589
590                 dmin = abs(prop.h_e[0] - prop.h_e[1]) / 200e-3;
591
592                 // :9: Diffraction coefficients, page 4
593                 /*
594                  * This is the region beyond the smooth-earth horizon at
595                  * D_Lsa and short of where isotropic scatter takes over. It
596                  * is a key to central region and associated coefficients
597                  * must always be computed.
598                  */
599                 q = adiff(0.0, prop);
600                 xae = pow(prop.k * prop.gamma_e * prop.gamma_e, -THIRD);  // [Alg 4.2]
601                 d3 = mymax(prop.d_Ls, 1.3787 * xae + prop.d_L);    // [Alg 4.3]
602                 d4 = d3 + 2.7574 * xae;                            // [Alg 4.4]
603                 a3 = adiff(d3, prop);                              // [Alg 4.5]
604                 a4 = adiff(d4, prop);                              // [Alg 4.7]
605
606                 prop.emd = (a4 - a3) / (d4 - d3);                  // [Alg 4.8]
607                 prop.aed = a3 - prop.emd * d3;
608         }
609
610         if (prop.mdp >= 0) {
611                 prop.mdp = 0;
612                 prop.d = d;
613         }
614
615         if (prop.d > 0.0) {
616                 // :8: Check distance, page 3
617
618                 // Distance above 1000 km is guessing
619                 if (prop.d > 1000e3) {
620                         set_warn("Distance not optimal", 1);
621                         prop.kwx = mymax(prop.kwx, 1);
622                 }
623
624                 // Distance too small, use some indoor algorithm :-)
625                 if (prop.d < dmin) {
626                         set_warn("Distance too small", 3);
627                         prop.kwx = mymax(prop.kwx, 3);
628                 }
629
630                 // Distance above 2000 km is really bad, don't do that
631                 if (prop.d < 1e3 || prop.d > 2000e3) {
632                         set_warn("Distance out-of-bounds", 4);
633                         prop.kwx = 4;
634                 }
635         }
636
637         if (prop.d < prop.d_Ls) {
638                 // :15: Line-of-sight calculations, page 7
639                 if (!wlos) {
640                         // :16: Line-of-sight coefficients, page 7
641                         q = A_los(0.0, prop);
642                         d2 = prop.d_Ls;
643                         a2 = prop.aed + d2 * prop.emd;
644                         d0 = 1.908 * prop.k * prop.h_e[0] * prop.h_e[1];                // [Alg 4.38]
645
646                         if (prop.aed >= 0.0) {
647                                 d0 = mymin(d0, 0.5 * prop.d_L);                       // [Alg 4.28]
648                                 d1 = d0 + 0.25 * (prop.d_L - d0);                     // [Alg 4.29]
649                         } else {
650                                 d1 = mymax(-prop.aed / prop.emd, 0.25 * prop.d_L);  // [Alg 4.39]
651                         }
652
653                         a1 = A_los(d1, prop);  // [Alg 4.31]
654                         wq = false;
655
656                         if (d0 < d1) {
657                                 a0 = A_los(d0, prop); // [Alg 4.30]
658                                 q = log(d2 / d0);
659                                 prop.ak2 = mymax(0.0, ((d2 - d0) * (a1 - a0) - (d1 - d0) * (a2 - a0)) / ((d2 - d0) * log(d1 / d0) - (d1 - d0) * q)); // [Alg 4.32]
660                                 wq = prop.aed >= 0.0 || prop.ak2 > 0.0;
661
662                                 if (wq) {
663                                         prop.ak1 = (a2 - a0 - prop.ak2 * q) / (d2 - d0); // [Alg 4.33]
664
665                                         if (prop.ak1 < 0.0) {
666                                                 prop.ak1 = 0.0;                      // [Alg 4.36]
667                                                 prop.ak2 = FORTRAN_DIM(a2, a0) / q;  // [Alg 4.35]
668
669                                                 if (prop.ak2 == 0.0)                 // [Alg 4.37]
670                                                         prop.ak1 = prop.emd;
671                                         }
672                                 }
673                         }
674
675                         if (!wq) {
676                                 prop.ak1 = FORTRAN_DIM(a2, a1) / (d2 - d1);   // [Alg 4.40]
677                                 prop.ak2 = 0.0;                               // [Alg 4.41]
678
679                                 if (prop.ak1 == 0.0)                          // [Alg 4.37]
680                                         prop.ak1 = prop.emd;
681                         }
682
683                         prop.ael = a2 - prop.ak1 * d2 - prop.ak2 * log(d2); // [Alg 4.42]
684                         wlos = true;
685                 }
686
687                 if (prop.d > 0.0)
688                         // [Alg 4.1]
689                         /*
690                          * The reference attenuation is determined as a function of the distance
691                          * d from 3 piecewise formulatios. This is part one.
692                          */
693                         prop.A_ref = prop.ael + prop.ak1 * prop.d + prop.ak2 * log(prop.d);
694         }
695
696         if (prop.d <= 0.0 || prop.d >= prop.d_Ls) {
697                 // :20: Troposcatter calculations, page 9
698                 if (!wscat) {
699                         // :21: Troposcatter coefficients
700                         const double DS = 200e3;             // 200 km from [Alg 4.53]
701                         const double HS = 47.7;              // 47.7 m from [Alg 4.59]
702                         q = A_scat(0.0, prop);
703                         d5 = prop.d_L + DS;                  // [Alg 4.52]
704                         d6 = d5 + DS;                        // [Alg 4.53]
705                         a6 = A_scat(d6, prop);                // [Alg 4.54]
706                         a5 = A_scat(d5, prop);                // [Alg 4.55]
707
708                         if (a5 < 1000.0) {
709                                 prop.ems = (a6 - a5) / DS;   // [Alg 4.57]
710                                 prop.dx = mymax(prop.d_Ls,   // [Alg 4.58]
711                                                  mymax(prop.d_L + 0.3 * xae * log(HS * prop.k),
712                                                  (a5 - prop.aed - prop.ems * d5) / (prop.emd - prop.ems)));
713                                 prop.aes = (prop.emd - prop.ems) * prop.dx + prop.aed; // [Alg 4.59]
714                         } else {
715                                 prop.ems = prop.emd;
716                                 prop.aes = prop.aed;
717                                 prop.dx = 10.e6;             // [Alg 4.56]
718                         }
719                         wscat = true;
720                 }
721
722                 // [Alg 4.1], part two and three.
723                 if (prop.d > prop.dx)
724                         prop.A_ref = prop.aes + prop.ems * prop.d;  // scatter region
725                 else
726                         prop.A_ref = prop.aed + prop.emd * prop.d;  // diffraction region
727         }
728
729         prop.A_ref = mymax(prop.A_ref, 0.0);
730 }
731
732
733 /*****************************************************************************/
734
735 // :27: Variablility parameters, page 13
736
737 struct propv_type {
738         // Input:
739         int lvar;    // control switch for initialisation and/or recalculation
740                      // 0: no recalculations needed
741                      // 1: distance changed
742                      // 2: antenna heights changed
743                      // 3: frequency changed
744                      // 4: mdvar changed
745                      // 5: climate changed, or initialize everything
746         int mdvar;   // desired mode of variability
747         int klim;    // climate indicator
748         // Output
749         double sgc;  // standard deviation of situation variability (confidence)
750 };
751
752
753 #ifdef WITH_QRLA
754 // :40: Prepare model for area prediction mode, page 21
755 /*
756  * This is used to prepare the model in the area prediction mode. Normally,
757  * one first calls qlrps() and then qlra(). Before calling the latter, one
758  * should have defined in prop_type the antenna heights @hg, the terrain
759  * irregularity parameter @dh , and (probably through qlrps() the variables
760  * @wn, @ens, @gme, and @zgnd. The input @kst will define siting criteria
761  * for the terminals, @klimx the climate, and @mdvarx the mode of
762  * variability. If @klimx <= 0 or @mdvarx < 0 the associated parameters
763  * remain unchanged.
764  */
765 static
766 void qlra(int kst[], int klimx, int mdvarx, prop_type &prop, propv_type &propv)
767 {
768         double q;
769
770         for (int j = 0; j < 2; ++j) {
771                 if (kst[j] <= 0)
772                         // [Alg 3.1]
773                         prop.h_e[j] = prop.h_g[j];
774                 else {
775                         q = 4.0;
776
777                         if (kst[j] != 1)
778                                 q = 9.0;
779
780                         if (prop.h_g[j] < 5.0)
781                                 q *= sin(0.3141593 * prop.h_g[j]);
782
783                         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)));
784                 }
785
786                 // [Alg 3.3], upper function. q is d_Ls_j
787                 const double H_3 = 5; // 5m from [Alg 3.3]
788                 q = sqrt(2.0 * prop.h_e[j] / prop.gamma_e);
789                 prop.d_Lj[j] = q * exp(-0.07 * sqrt(prop.delta_h / mymax(prop.h_e[j], H_3)));
790                 // [Alg 3.4]
791                 prop.theta_ej[j] = (0.65 * prop.delta_h * (q / prop.d_Lj[j] - 1.0) - 2.0 * prop.h_e[j]) / q;
792         }
793
794         prop.mdp = 1;
795         propv.lvar = mymax(propv.lvar, 3);
796
797         if (mdvarx >= 0) {
798                 propv.mdvar = mdvarx;
799                 propv.lvar = mymax(propv.lvar, 4);
800         }
801
802         if (klimx > 0) {
803                 propv.klim = klimx;
804                 propv.lvar = 5;
805         }
806 }
807 #endif
808
809
810 // :51: Inverse of standard normal complementary probability
811 /*
812  * The approximation is due to C. Hastings, Jr. ("Approximations for digital
813  * computers," Princeton Univ. Press, 1955) and the maximum error should be
814  * 4.5e-4.
815  */
816 #if 1
817 static
818 double qerfi(double q)
819 {
820         double x, t, v;
821         const double c0 = 2.515516698;
822         const double c1 = 0.802853;
823         const double c2 = 0.010328;
824         const double d1 = 1.432788;
825         const double d2 = 0.189269;
826         const double d3 = 0.001308;
827
828         x = 0.5 - q;
829         t = mymax(0.5 - fabs(x), 0.000001);
830         t = sqrt(-2.0 * log(t));
831         v = t - ((c2 * t + c1) * t + c0) / (((d3 * t + d2) * t + d1) * t + 1.0);
832
833         if (x < 0.0)
834                 v = -v;
835
836         return v;
837 }
838 #endif
839
840
841 // :41: preparatory routine, page 20
842 /*
843  * This subroutine converts
844  *   frequency @fmhz
845  *   surface refractivity reduced to sea level @en0
846  *   general system elevation @zsys
847  *   polarization and ground constants @eps, @sgm
848  * to
849  *   wave number @wn,
850  *   surface refractivity @ens
851  *   effective earth curvature @gme
852  *   surface impedance @zgnd
853  * It may be used with either the area prediction or the point-to-point mode.
854  */
855 #if 1
856 static
857 void qlrps(double fmhz, double zsys, double en0, int ipol, double eps, double sgm, prop_type &prop)
858
859 {
860         const double Z_1 = 9.46e3; // 9.46 km       from [Alg 1.2]
861         const double gma = 157e-9; // 157e-9 1/m    from [Alg 1.3]
862         const double N_1 = 179.3;  // 179.3 N-units from [Alg 1.3]
863         const double Z_0 = 376.62; // 376.62 Ohm    from [Alg 1.5]
864
865         // Frequecy -> Wave number
866         prop.k = fmhz / f_0;                                      // [Alg 1.1]
867
868         // Surface refractivity
869         prop.N_s = en0;
870         if (zsys != 0.0)
871                 prop.N_s *= exp(-zsys / Z_1);                      // [Alg 1.2]
872
873         // Earths effective curvature
874         prop.gamma_e = gma * (1.0 - 0.04665 * exp(prop.N_s / N_1));    // [Alg 1.3]
875
876         // Surface transfer impedance
877         complex<double> zq, prop_zgnd(prop.Z_g_real, prop.Z_g_imag);
878         zq = complex<double> (eps, Z_0 * sgm / prop.k);        // [Alg 1.5]
879         prop_zgnd = sqrt(zq - 1.0);
880
881         // adjust surface transfer impedance for Polarization
882         if (ipol != 0.0)
883                 prop_zgnd = prop_zgnd / zq;                        // [Alg 1.4]
884
885         prop.Z_g_real = prop_zgnd.real();
886         prop.Z_g_imag = prop_zgnd.imag();
887 }
888 #endif
889
890
891 // :30: Function curv, page 15
892 static
893 double curve(double const &c1, double const &c2, double const &x1, double const &x2, double const &x3, double const &de)
894 {
895         double temp1, temp2;
896
897         temp1 = (de - x2) / x3;
898         temp2 = de / x1;
899
900         temp1 *= temp1;
901         temp2 *= temp2;
902
903         return (c1 + c2 / (1.0 + temp1)) * temp2 / (1.0 + temp2);
904 }
905
906
907 // :28: Area variablity, page 14
908 #if 1
909 static
910 double avar(double zzt, double zzl, double zzc, prop_type &prop, propv_type &propv)
911 {
912         static int kdv;
913         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;
914
915         // :29: Climatic constants, page 15
916         // Indexes are:
917         //   0: equatorial
918         //   1: continental suptropical
919         //   2: maritime subtropical
920         //   3: desert
921         //   4: continental temperature
922         //   5: maritime over land
923         //   6: maritime over sea
924         //                      equator  contsup maritsup   desert conttemp mariland  marisea
925         const double bv1[7]  = {  -9.67,   -0.62,    1.26,   -9.21,   -0.62,   -0.39,    3.15};
926         const double bv2[7]  = {   12.7,    9.19,    15.5,    9.05,    9.19,    2.86,   857.9};
927         const double xv1[7]  = {144.9e3, 228.9e3, 262.6e3,  84.1e3, 228.9e3, 141.7e3, 2222.e3};
928         const double xv2[7]  = {190.3e3, 205.2e3, 185.2e3, 101.1e3, 205.2e3, 315.9e3, 164.8e3};
929         const double xv3[7]  = {133.8e3, 143.6e3,  99.8e3,  98.6e3, 143.6e3, 167.4e3, 116.3e3};
930         const double bsm1[7] = {   2.13,    2.66,    6.11,    1.98,    2.68,    6.86,    8.51};
931         const double bsm2[7] = {  159.5,    7.67,    6.65,   13.11,    7.16,   10.38,   169.8};
932         const double xsm1[7] = {762.2e3, 100.4e3, 138.2e3, 139.1e3,  93.7e3, 187.8e3, 609.8e3};
933         const double xsm2[7] = {123.6e3, 172.5e3, 242.2e3, 132.7e3, 186.8e3, 169.6e3, 119.9e3};
934         const double xsm3[7] = { 94.5e3, 136.4e3, 178.6e3, 193.5e3, 133.5e3, 108.9e3, 106.6e3};
935         const double bsp1[7] = {   2.11,    6.87,   10.08,    3.68,    4.75,    8.58,    8.43};
936         const double bsp2[7] = {  102.3,   15.53,    9.60,   159.3,    8.12,   13.97,    8.19};
937         const double xsp1[7] = {636.9e3, 138.7e3, 165.3e3, 464.4e3,  93.2e3, 216.0e3, 136.2e3};
938         const double xsp2[7] = {134.8e3, 143.7e3, 225.7e3,  93.1e3, 135.9e3, 152.0e3, 188.5e3};
939         const double xsp3[7] = { 95.6e3,  98.6e3, 129.7e3,  94.2e3, 113.4e3, 122.7e3, 122.9e3};
940         // bds1 -> is similar to C_D from table 5.1 at [Alg 5.8]
941         // bzd1 -> is similar to z_D from table 5.1 at [Alg 5.8]
942         const double bsd1[7] = {  1.224,   0.801,   1.380,   1.000,   1.224,   1.518,   1.518};
943         const double bzd1[7] = {  1.282,   2.161,   1.282,    20.0,   1.282,   1.282,   1.282};
944         const double bfm1[7] = {    1.0,     1.0,     1.0,     1.0,    0.92,     1.0,    1.0};
945         const double bfm2[7] = {    0.0,     0.0,     0.0,     0.0,    0.25,     0.0,    0.0};
946         const double bfm3[7] = {    0.0,     0.0,     0.0,     0.0,    1.77,     0.0,    0.0};
947         const double bfp1[7] = {    1.0,    0.93,     1.0,    0.93,    0.93,     1.0,    1.0};
948         const double bfp2[7] = {    0.0,    0.31,     0.0,    0.19,    0.31,     0.0,    0.0};
949         const double bfp3[7] = {    0.0,    2.00,     0.0,    1.79,    2.00,     0.0,    0.0};
950         const double rt = 7.8, rl = 24.0;
951         static bool no_location_variability, no_situation_variability;
952         double avarv, q, vs, zt, zl, zc;
953         double sgt, yr;
954         int temp_klim;
955
956         if (propv.lvar > 0) {
957                 // :31: Setup variablity constants, page 16
958                 switch (propv.lvar) {
959                 default:
960                         // Initialization or climate change
961
962                         // if climate is wrong, use some "continental temperature" as default
963                         // and set error indicator
964                         if (propv.klim <= 0 || propv.klim > 7) {
965                                 propv.klim = 5;
966                                 prop.kwx = mymax(prop.kwx, 2);
967                                 set_warn("Climate index set to continental", 2);
968                         }
969
970                         // convert climate number into index into the climate tables
971                         temp_klim = propv.klim - 1;
972
973                         // :32: Climatic coefficients, page 17
974                         cv1 = bv1[temp_klim];
975                         cv2 = bv2[temp_klim];
976                         yv1 = xv1[temp_klim];
977                         yv2 = xv2[temp_klim];
978                         yv3 = xv3[temp_klim];
979                         csm1 = bsm1[temp_klim];
980                         csm2 = bsm2[temp_klim];
981                         ysm1 = xsm1[temp_klim];
982                         ysm2 = xsm2[temp_klim];
983                         ysm3 = xsm3[temp_klim];
984                         csp1 = bsp1[temp_klim];
985                         csp2 = bsp2[temp_klim];
986                         ysp1 = xsp1[temp_klim];
987                         ysp2 = xsp2[temp_klim];
988                         ysp3 = xsp3[temp_klim];
989                         csd1 = bsd1[temp_klim];
990                         zd = bzd1[temp_klim];
991                         cfm1 = bfm1[temp_klim];
992                         cfm2 = bfm2[temp_klim];
993                         cfm3 = bfm3[temp_klim];
994                         cfp1 = bfp1[temp_klim];
995                         cfp2 = bfp2[temp_klim];
996                         cfp3 = bfp3[temp_klim];
997                         // fall throught
998
999                 case 4:
1000                         // :33: Mode of variablity coefficients, page 17
1001
1002                         // This code means that propv.mdvar can be
1003                         //  0 ..  3
1004                         // 10 .. 13, then no_location_variability is set              (no location variability)
1005                         // 20 .. 23, then no_situation_variability is set             (no situatian variability)
1006                         // 30 .. 33, then no_location_variability and no_situation_variability is set
1007                         kdv = propv.mdvar;
1008                         no_situation_variability = kdv >= 20;
1009                         if (no_situation_variability)
1010                                 kdv -= 20;
1011
1012                         no_location_variability = kdv >= 10;
1013                         if (no_location_variability)
1014                                 kdv -= 10;
1015
1016                         if (kdv < 0 || kdv > 3) {
1017                                 kdv = 0;
1018                                 set_warn("kdv set to 0", 2);
1019                                 prop.kwx = mymax(prop.kwx, 2);
1020                         }
1021
1022                         // fall throught
1023
1024                 case 3:
1025                         // :34: Frequency coefficients, page 18
1026                         q = log(0.133 * prop.k);
1027                         gm = cfm1 + cfm2 / ((cfm3 * q * cfm3 * q) + 1.0);
1028                         gp = cfp1 + cfp2 / ((cfp3 * q * cfp3 * q) + 1.0);
1029                         // fall throught
1030
1031                 case 2:
1032                         // :35: System coefficients, page 18
1033                         // [Alg 5.3], effective distance
1034                         {
1035                         const double a_1 = 9000e3; // 9000 km from [[Alg 5.3]
1036                         //XXX I don't have any idea how they made up the third summand,
1037                         //XXX text says    a_1 * pow(k * D_1, -THIRD)
1038                         //XXX with const double D_1 = 1266; // 1266 km
1039                         dexa = sqrt(2*a_1 * prop.h_e[0]) +
1040                                sqrt(2*a_1 * prop.h_e[1]) +
1041                                pow((575.7e12 / prop.k), THIRD);
1042                         }
1043                         // fall throught
1044
1045                 case 1:
1046                         // :36: Distance coefficients, page 18
1047                         {
1048                         // [Alg 5.4]
1049                         const double D_0 = 130e3; // 130 km from [Alg 5.4]
1050                         if (prop.d < dexa)
1051                                 de = D_0 * prop.d / dexa;
1052                         else
1053                                 de = D_0 + prop.d - dexa;
1054                         }
1055                 }
1056                 /*
1057                  * Quantiles of time variability are computed using a
1058                  * variation of the methods described in Section 10 and
1059                  * Annex III.7 of NBS~TN101, and also in CCIR
1060                  * Report {238-3}. Those methods speak of eight or nine
1061                  * discrete radio climates, of which seven have been
1062                  * documented with corresponding empirical curves. It is
1063                  * these empirical curves to which we refer below. They are
1064                  * all curves of quantiles of deviations versus the
1065                  * effective distance @de.
1066                  */
1067                 vmd = curve(cv1, cv2, yv1, yv2, yv3, de);             // [Alg 5.5]
1068                 // [Alg 5.7], the slopes or "pseudo-standard deviations":
1069                 // sgtm -> \sigma T-
1070                 // sgtp -> \sigma T+
1071                 sgtm = curve(csm1, csm2, ysm1, ysm2, ysm3, de) * gm;
1072                 sgtp = curve(csp1, csp2, ysp1, ysp2, ysp3, de) * gp;
1073                 // [Alg 5.8], ducting, "sgtd" -> \sigma TD
1074                 sgtd = sgtp * csd1;
1075                 tgtd = (sgtp - sgtd) * zd;
1076
1077                 // Location variability
1078                 if (no_location_variability) {
1079                         sgl = 0.0;
1080                 } else {
1081                         // Alg [3.9]
1082                         q = (1.0 - 0.8 * exp(-prop.d / 50e3)) * prop.delta_h * prop.k;
1083                         // [Alg 5.9]
1084                         sgl = 10.0 * q / (q + 13.0);
1085                 }
1086
1087                 // Situation variability
1088                 if (no_situation_variability) {
1089                         vs0 = 0.0;
1090                 } else {
1091                         const double D = 100e3; // 100 km
1092                         vs0 = (5.0 + 3.0 * exp(-de / D));         // [Alg 5.10]
1093                         vs0 *= vs0;
1094                 }
1095
1096                 // Mark all constants as initialized
1097                 propv.lvar = 0;
1098         }
1099
1100
1101         // :37: Correct normal deviates, page 19
1102         zt = zzt;
1103         zl = zzl;
1104         zc = zzc;
1105         // kdv is derived from prop.mdvar
1106         switch (kdv) {
1107         case 0:
1108                 // Single message mode. Time, location and situation variability
1109                 // are combined to form a confidence level.
1110                 zt = zc;
1111                 zl = zc;
1112                 break;
1113
1114         case 1:
1115                 // Individual mode. Reliability is given by time
1116                 // variability. Confidence. is a combination of location
1117                 // and situation variability.
1118                 zl = zc;
1119                 break;
1120
1121         case 2:
1122                 // Mobile modes. Reliability is a combination of time and
1123                 // location variability. Confidence. is given by the
1124                 // situation variability.
1125                 zl = zt;
1126                 // case 3: Broadcast mode. like avar(zt, zl, zc).
1127                 // Reliability is given by the two-fold statement of at
1128                 // least qT of the time in qL of the locations. Confidence.
1129                 // is given by the situation variability.
1130         }
1131         if (fabs(zt) > 3.1 || fabs(zl) > 3.1 || fabs(zc) > 3.1) {
1132                 set_warn("Situations variables not optimal", 1);
1133                 prop.kwx = mymax(prop.kwx, 1);
1134         }
1135
1136
1137         // :38: Resolve standard deviations, page 19
1138         if (zt < 0.0)
1139                 sgt = sgtm;
1140         else
1141         if (zt <= zd)
1142                 sgt = sgtp;
1143         else
1144                 sgt = sgtd + tgtd / zt;
1145         // [Alg 5.11], situation variability
1146         vs = vs0 + (sgt * zt * sgt * zt) / (rt + zc * zc) + (sgl * zl * sgl * zl) / (rl + zc * zc);
1147
1148
1149         // :39: Resolve deviations, page 19
1150         if (kdv == 0) {
1151                 yr = 0.0;
1152                 propv.sgc = sqrt(sgt * sgt + sgl * sgl + vs);
1153         } else
1154         if (kdv == 1) {
1155                 yr = sgt * zt;
1156                 propv.sgc = sqrt(sgl * sgl + vs);
1157         } else
1158         if (kdv == 2) {
1159                 yr = sqrt(sgt * sgt + sgl * sgl) * zt;
1160                 propv.sgc = sqrt(vs);
1161         } else {
1162                 yr = sgt * zt + sgl * zl;
1163                 propv.sgc = sqrt(vs);
1164         }
1165
1166         // [Alg 5.1], area variability
1167         avarv = prop.A_ref - vmd - yr - propv.sgc * zc;
1168
1169         // [Alg 5.2]
1170         if (avarv < 0.0)
1171                 avarv = avarv * (29.0 - avarv) / (29.0 - 10.0 * avarv);
1172
1173         return avarv;
1174 }
1175 #endif
1176
1177 // :45: Find to horizons, page 24
1178 /*
1179  * Here we use the terrain profile @pfl to find the two horizons. Output
1180  * consists of the horizon distances @dl and the horizon take-off angles
1181  * @the. If the path is line-of-sight, the routine sets both horizon
1182  * distances equal to @dist.
1183  */
1184 static
1185 void hzns(double pfl[], prop_type &prop)
1186 {
1187         bool wq;
1188         int np;
1189         double xi, za, zb, qc, q, sb, sa;
1190
1191         np = (int)pfl[0];
1192         xi = pfl[1];
1193         za = pfl[2] + prop.h_g[0];
1194         zb = pfl[np+2] + prop.h_g[1];
1195         qc = 0.5 * prop.gamma_e;
1196         q = qc * prop.d;
1197         prop.theta_ej[1] = (zb - za) / prop.d;
1198         prop.theta_ej[0] = prop.theta_ej[1] - q;
1199         prop.theta_ej[1] = -prop.theta_ej[1] - q;
1200         prop.d_Lj[0] = prop.d;
1201         prop.d_Lj[1] = prop.d;
1202
1203         if (np >= 2) {
1204                 sa = 0.0;
1205                 sb = prop.d;
1206                 wq = true;
1207
1208                 for (int i = 1; i < np; i++) {
1209                         sa += xi;
1210                         sb -= xi;
1211                         q = pfl[i+2] - (qc * sa + prop.theta_ej[0]) * sa - za;
1212
1213                         if (q > 0.0) {
1214                                 prop.theta_ej[0] += q / sa;
1215                                 prop.d_Lj[0] = sa;
1216                                 wq = false;
1217                         }
1218
1219                         if (!wq) {
1220                                 q = pfl[i+2] - (qc * sb + prop.theta_ej[1]) * sb - zb;
1221
1222                                 if (q > 0.0) {
1223                                         prop.theta_ej[1] += q / sb;
1224                                         prop.d_Lj[1] = sb;
1225                                 }
1226                         }
1227                 }
1228         }
1229 }
1230
1231
1232 // :53: Linear least square fit, page 28
1233 static
1234 void zlsq1(double z[], const double &x1, const double &x2, double& z0, double& zn)
1235
1236 {
1237         double xn, xa, xb, x, a, b;
1238         int n, ja, jb;
1239
1240         xn = z[0];
1241         xa = int(FORTRAN_DIM(x1 / z[1], 0.0));
1242         xb = xn - int(FORTRAN_DIM(xn, x2 / z[1]));
1243
1244         if (xb <= xa) {
1245                 xa = FORTRAN_DIM(xa, 1.0);
1246                 xb = xn - FORTRAN_DIM(xn, xb + 1.0);
1247         }
1248
1249         ja = (int)xa;
1250         jb = (int)xb;
1251         n = jb - ja;
1252         xa = xb - xa;
1253         x = -0.5 * xa;
1254         xb += x;
1255         a = 0.5 * (z[ja+2] + z[jb+2]);
1256         b = 0.5 * (z[ja+2] - z[jb+2]) * x;
1257
1258         for (int i = 2; i <= n; ++i) {
1259                 ++ja;
1260                 x += 1.0;
1261                 a += z[ja+2];
1262                 b += z[ja+2] * x;
1263         }
1264
1265         a /= xa;
1266         b = b * 12.0 / ((xa * xa + 2.0) * xa);
1267         z0 = a - b * xb;
1268         zn = a + b * (xn - xb);
1269 }
1270
1271
1272 // :52: Provide a quantile and reorders array @a, page 27
1273 static
1274 double qtile(const int &nn, double a[], const int &ir)
1275 {
1276         double q = 0.0, r;                   /* Initializations by KD2BD */
1277         int m, n, i, j, j1 = 0, i0 = 0, k;   /* Initializations by KD2BD */
1278         bool done = false;
1279         bool goto10 = true;
1280
1281         m = 0;
1282         n = nn;
1283         k = mymin(mymax(0, ir), n);
1284
1285         while (!done) {
1286                 if (goto10) {
1287                         q = a[k];
1288                         i0 = m;
1289                         j1 = n;
1290                 }
1291
1292                 i = i0;
1293
1294                 while (i <= n && a[i] >= q)
1295                         i++;
1296
1297                 if (i > n)
1298                         i = n;
1299                 j = j1;
1300
1301                 while (j >= m && a[j] <= q)
1302                         j--;
1303
1304                 if (j < m)
1305                         j = m;
1306
1307                 if (i < j) {
1308                         r = a[i];
1309                         a[i] = a[j];
1310                         a[j] = r;
1311                         i0 = i + 1;
1312                         j1 = j - 1;
1313                         goto10 = false;
1314                 } else
1315                 if (i < k) {
1316                         a[k] = a[i];
1317                         a[i] = q;
1318                         m = i + 1;
1319                         goto10 = true;
1320                 } else
1321                 if (j > k) {
1322                         a[k] = a[j];
1323                         a[j] = q;
1324                         n = j - 1;
1325                         goto10 = true;
1326                 } else {
1327                         done = true;
1328                 }
1329         }
1330
1331         return q;
1332 }
1333
1334
1335 #ifdef WITH_QERF
1336 // :50: Standard normal complementary probability, page 26
1337 static
1338 double qerf(const double &z)
1339 {
1340         const double b1 = 0.319381530, b2 = -0.356563782, b3 = 1.781477937;
1341         const double b4 = -1.821255987, b5 = 1.330274429;
1342         const double rp = 4.317008, rrt2pi = 0.398942280;
1343         double t, x, qerfv;
1344
1345         x = z;
1346         t = fabs(x);
1347
1348         if (t >= 10.0)
1349                 qerfv = 0.0;
1350         else {
1351                 t = rp / (t + rp);
1352                 qerfv = exp(-0.5 * x * x) * rrt2pi * ((((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t;
1353         }
1354
1355         if (x < 0.0)
1356                 qerfv = 1.0 - qerfv;
1357
1358         return qerfv;
1359 }
1360 #endif
1361
1362
1363 // :48: Find interdecile range of elevations, page 25
1364 /*
1365  * Using the terrain profile @pfl we find \Delta h, the interdecile range of
1366  * elevations between the two points @x1 and @x2.
1367  */
1368 static
1369 double dlthx(double pfl[], const double &x1, const double &x2)
1370 {
1371         int np, ka, kb, n, k, j;
1372         double dlthxv, sn, xa, xb;
1373         double *s;
1374
1375         np = (int)pfl[0];
1376         xa = x1 / pfl[1];
1377         xb = x2 / pfl[1];
1378         dlthxv = 0.0;
1379
1380         if (xb - xa < 2.0) // exit out
1381                 return dlthxv;
1382
1383         ka = (int)(0.1 * (xb - xa + 8.0));
1384         ka = mymin(mymax(4, ka), 25);
1385         n = 10 * ka - 5;
1386         kb = n - ka + 1;
1387         sn = n - 1;
1388         s = new double[n+2];
1389         s[0] = sn;
1390         s[1] = 1.0;
1391         xb = (xb - xa) / sn;
1392         k = (int)(xa + 1.0);
1393         xa -= (double)k;
1394
1395         for (j = 0; j < n; j++) {
1396                 // Reduce
1397                 while (xa > 0.0 && k < np) {
1398                         xa -= 1.0;
1399                         ++k;
1400                 }
1401
1402                 s[j+2] = pfl[k+2] + (pfl[k+2] - pfl[k+1]) * xa;
1403                 xa = xa + xb;
1404         }
1405
1406         zlsq1(s, 0.0, sn, xa, xb);
1407         xb = (xb - xa) / sn;
1408
1409         for (j = 0; j < n; j++) {
1410                 s[j+2] -= xa;
1411                 xa = xa + xb;
1412         }
1413
1414         dlthxv = qtile(n - 1, s + 2, ka - 1) - qtile(n - 1, s + 2, kb - 1);
1415         dlthxv /= 1.0 - 0.8 * exp(-(x2 - x1) / 50.0e3);
1416         delete[] s;
1417
1418         return dlthxv;
1419 }
1420
1421
1422 // :41: Prepare model for point-to-point operation, page 22
1423 /*
1424  * This mode requires the terrain profile lying between the terminals. This
1425  * should be a sequence of surface elevations at points along the great
1426  * circle path joining the two points. It should start at the ground beneath
1427  * the first terminal and end at the ground beneath the second. In the
1428  * present routine it is assumed that the elevations are equispaced along
1429  * the path. They are stored in the array @pfl along with two defining
1430  * parameters.
1431  *
1432  * We will have
1433  *   pfl[0] = np, the number of points in the path
1434  *   pfl[1] = xi, the length of each increment
1435 */
1436 #if 1
1437 static
1438 void qlrpfl(double pfl[], int klimx, int mdvarx, prop_type &prop, propv_type &propv)
1439 {
1440         int np, j;
1441         double xl[2], q, za, zb;
1442
1443         prop.d = pfl[0] * pfl[1];
1444         np = (int)pfl[0];
1445
1446         // :44: determine horizons and dh from pfl, page 23
1447         hzns(pfl, prop);
1448         for (j = 0; j < 2; j++)
1449                 xl[j] = mymin(15.0 * prop.h_g[j], 0.1 * prop.d_Lj[j]);
1450
1451         xl[1] = prop.d - xl[1];
1452         prop.delta_h = dlthx(pfl, xl[0], xl[1]);
1453
1454         if (prop.d_Lj[0] + prop.d_Lj[1] > 1.5 * prop.d) {
1455                 // :45: Redo line-of-sight horizons, page 23
1456
1457                 /*
1458                  * If the path is line-of-sight, we still need to know where
1459                  * the horizons might have been, and so we turn to
1460                  * techniques used in area prediction mode.
1461                  */
1462                 zlsq1(pfl, xl[0], xl[1], za, zb);
1463                 prop.h_e[0] = prop.h_g[0] + FORTRAN_DIM(pfl[2], za);
1464                 prop.h_e[1] = prop.h_g[1] + FORTRAN_DIM(pfl[np+2], zb);
1465
1466                 for (j = 0; j < 2; j++)
1467                         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)));
1468
1469                 q = prop.d_Lj[0] + prop.d_Lj[1];
1470
1471                 if (q <= prop.d) {
1472                         q = ((prop.d / q) * (prop.d / q));
1473
1474                         for (j = 0; j < 2; j++) {
1475                                 prop.h_e[j] *= q;
1476                                 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)));
1477                         }
1478                 }
1479
1480                 for (j = 0; j < 2; j++) {
1481                         q = sqrt(2.0 * prop.h_e[j] / prop.gamma_e);
1482                         prop.theta_ej[j] = (0.65 * prop.delta_h * (q / prop.d_Lj[j] - 1.0) - 2.0 * prop.h_e[j]) / q;
1483                 }
1484         } else {
1485                 // :46: Transhorizon effective heights, page 23
1486                 zlsq1(pfl, xl[0], 0.9 * prop.d_Lj[0], za, q);
1487                 zlsq1(pfl, prop.d - 0.9 * prop.d_Lj[1], xl[1], q, zb);
1488                 prop.h_e[0] = prop.h_g[0] + FORTRAN_DIM(pfl[2], za);
1489                 prop.h_e[1] = prop.h_g[1] + FORTRAN_DIM(pfl[np+2], zb);
1490         }
1491
1492         prop.mdp = -1;
1493         propv.lvar = mymax(propv.lvar, 3);
1494
1495         if (mdvarx >= 0) {
1496                 propv.mdvar = mdvarx;
1497                 propv.lvar = mymax(propv.lvar, 4);
1498         }
1499
1500         if (klimx > 0) {
1501                 propv.klim = klimx;
1502                 propv.lvar = 5;
1503         }
1504
1505         lrprop(0.0, prop);
1506 }
1507 #endif
1508
1509
1510 //********************************************************
1511 //* Point-To-Point Mode Calculations                     *
1512 //********************************************************
1513
1514
1515 #ifdef WITH_POINT_TO_POINT
1516 #include <string.h>
1517 static
1518 void point_to_point(double elev[],
1519                     double tht_m,              // Transceiver above ground level
1520                     double rht_m,              // Receiver above groud level
1521                     double eps_dielect,        // Earth dielectric constant (rel. permittivity)
1522                     double sgm_conductivity,   // Earth conductivity (Siemens/m)
1523                     double eno,                // Atmospheric bending const, n-Units
1524                     double frq_mhz,
1525                     int radio_climate,         // 1..7
1526                     int pol,                   // 0 horizontal, 1 vertical
1527                     double conf,               // 0.01 .. .99, Fractions of situations
1528                     double rel,                // 0.01 .. .99, Fractions of time
1529                     double &dbloss,
1530                     char *strmode,
1531                     int &p_mode,                                // propagation mode selector
1532                     double (&horizons)[2],                      // horizon distances
1533                     int &errnum)
1534 {
1535         // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1536         //                4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1537         //                7-Maritime Temperate, Over Sea
1538         // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
1539         // errnum: 0- No Error.
1540         //         1- Warning: Some parameters are nearly out of range.
1541         //                     Results should be used with caution.
1542         //         2- Note: Default parameters have been substituted for impossible ones.
1543         //         3- Warning: A combination of parameters is out of range.
1544         //                     Results are probably invalid.
1545         //         Other-  Warning: Some parameters are out of range.
1546         //                          Results are probably invalid.
1547
1548         prop_type   prop;
1549         propv_type  propv;
1550
1551         double zsys = 0;
1552         double zc, zr;
1553         double q = eno;
1554         long ja, jb, i, np;
1555         double dkm, xkm;
1556         double fs;
1557
1558         prop.h_g[0] = tht_m;          // Tx height above ground level
1559         prop.h_g[1] = rht_m;          // Rx height above ground level
1560         propv.klim = radio_climate;
1561         prop.kwx = 0;                // no error yet
1562         propv.lvar = 5;              // initialize all constants
1563         prop.mdp = -1;               // point-to-point
1564         zc = qerfi(conf);
1565         zr = qerfi(rel);
1566         np = (long)elev[0];
1567         dkm = (elev[1] * elev[0]) / 1000.0;
1568         xkm = elev[1] / 1000.0;
1569
1570         ja = (long)(3.0 + 0.1 * elev[0]);
1571         jb = np - ja + 6;
1572         for (i = ja - 1; i < jb; ++i)
1573                 zsys += elev[i];
1574         zsys /= (jb - ja + 1);
1575
1576         propv.mdvar = 12;
1577         qlrps(frq_mhz, zsys, q, pol, eps_dielect, sgm_conductivity, prop);
1578         qlrpfl(elev, propv.klim, propv.mdvar, prop, propv);
1579         fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
1580         q = prop.d - prop.d_L;
1581
1582         horizons[0] = 0.0;
1583         horizons[1] = 0.0;
1584         if (int(q) < 0.0) {
1585                 strcpy(strmode, "Line-Of-Sight Mode");
1586                 p_mode = 0;
1587         } else {
1588                 if (int(q) == 0.0) {
1589                         strcpy(strmode, "Single Horizon");
1590                         horizons[0] = prop.d_Lj[0];
1591                         p_mode = 1;
1592                 }
1593
1594                 else {
1595                         if (int(q) > 0.0) {
1596                                 strcpy(strmode, "Double Horizon");
1597                                 horizons[0] = prop.d_Lj[0];
1598                                 horizons[1] = prop.d_Lj[1];
1599                                 p_mode = 1;
1600                         }
1601                 }
1602
1603                 if (prop.d <= prop.d_Ls || prop.d <= prop.dx) {
1604                         strcat(strmode, ", Diffraction Dominant");
1605                         p_mode = 1;
1606                 }
1607
1608                 else {
1609                         if (prop.d > prop.dx) {
1610                                 strcat(strmode, ", Troposcatter Dominant");
1611                                 p_mode = 2;
1612                         }
1613                 }
1614         }
1615
1616         dbloss = avar(zr, 0.0, zc, prop, propv) + fs;
1617         errnum = prop.kwx;
1618 }
1619 #endif
1620
1621
1622 #ifdef WITH_POINT_TO_POINT_MDH
1623 static
1624 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)
1625 {
1626         // pol: 0-Horizontal, 1-Vertical
1627         // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1628         //                4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1629         //                7-Maritime Temperate, Over Sea
1630         // timepct, locpct, confpct: .01 to .99
1631         // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
1632         // propmode:  Value   Mode
1633         //             -1     mode is undefined
1634         //              0     Line of Sight
1635         //              5     Single Horizon, Diffraction
1636         //              6     Single Horizon, Troposcatter
1637         //              9     Double Horizon, Diffraction
1638         //             10     Double Horizon, Troposcatter
1639         // errnum: 0- No Error.
1640         //         1- Warning: Some parameters are nearly out of range.
1641         //                     Results should be used with caution.
1642         //         2- Note: Default parameters have been substituted for impossible ones.
1643         //         3- Warning: A combination of parameters is out of range.
1644         //                     Results are probably invalid.
1645         //         Other-  Warning: Some parameters are out of range.
1646         //                          Results are probably invalid.
1647
1648         prop_type   prop;
1649         propv_type  propv;
1650         propa_type  propa;
1651         double zsys = 0;
1652         double ztime, zloc, zconf;
1653         double q = eno;
1654         long ja, jb, i, np;
1655         double dkm, xkm;
1656         double fs;
1657
1658         propmode = -1; // mode is undefined
1659         prop.h_g[0] = tht_m;
1660         prop.h_g[1] = rht_m;
1661         propv.klim = radio_climate;
1662         prop.kwx = 0;
1663         propv.lvar = 5;
1664         prop.mdp = -1;
1665         ztime = qerfi(timepct);
1666         zloc = qerfi(locpct);
1667         zconf = qerfi(confpct);
1668
1669         np = (long)elev[0];
1670         dkm = (elev[1] * elev[0]) / 1000.0;
1671         xkm = elev[1] / 1000.0;
1672
1673         ja = (long)(3.0 + 0.1 * elev[0]);
1674         jb = np - ja + 6;
1675         for (i = ja - 1; i < jb; ++i)
1676                 zsys += elev[i];
1677         zsys /= (jb - ja + 1);
1678
1679         propv.mdvar = 12;
1680         qlrps(frq_mhz, zsys, q, pol, eps_dielect, sgm_conductivity, prop);
1681         qlrpfl(elev, propv.klim, propv.mdvar, prop, propa, propv);
1682         fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
1683         deltaH = prop.delta_h;
1684         q = prop.d - prop.d_L;
1685
1686         if (int(q) < 0.0) {
1687                 propmode = 0; // Line-Of-Sight Mode
1688         } else {
1689                 if (int(q) == 0.0)
1690                         propmode = 4; // Single Horizon
1691                 else
1692                         if (int(q) > 0.0)
1693                                 propmode = 8; // Double Horizon
1694
1695                 if (prop.d <= prop.d_Ls || prop.d <= prop.dx)
1696                         propmode += 1; // Diffraction Dominant
1697
1698                 else
1699                 if (prop.d > prop.dx)
1700                         propmode += 2; // Troposcatter Dominant
1701         }
1702
1703         dbloss = avar(ztime, zloc, zconf, prop, propv) + fs;  //avar(time,location,confidence)
1704         errnum = prop.kwx;
1705 }
1706 #endif
1707
1708
1709 #ifdef WITH_POINT_TO_POINT_DH
1710 static
1711 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)
1712 {
1713         // pol: 0-Horizontal, 1-Vertical
1714         // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1715         //                4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1716         //                7-Maritime Temperate, Over Sea
1717         // conf, rel: .01 to .99
1718         // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
1719         // errnum: 0- No Error.
1720         //         1- Warning: Some parameters are nearly out of range.
1721         //                     Results should be used with caution.
1722         //         2- Note: Default parameters have been substituted for impossible ones.
1723         //         3- Warning: A combination of parameters is out of range.
1724         //                     Results are probably invalid.
1725         //         Other-  Warning: Some parameters are out of range.
1726         //                          Results are probably invalid.
1727
1728         char strmode[100];
1729         prop_type   prop;
1730         propv_type  propv;
1731         propa_type  propa;
1732         double zsys = 0;
1733         double zc, zr;
1734         double q = eno;
1735         long ja, jb, i, np;
1736         double dkm, xkm;
1737         double fs;
1738
1739         prop.h_g[0] = tht_m;
1740         prop.h_g[1] = rht_m;
1741         propv.klim = radio_climate;
1742         prop.kwx = 0;
1743         propv.lvar = 5;
1744         prop.mdp = -1;
1745         zc = qerfi(conf);
1746         zr = qerfi(rel);
1747         np = (long)elev[0];
1748         dkm = (elev[1] * elev[0]) / 1000.0;
1749         xkm = elev[1] / 1000.0;
1750
1751         ja = (long)(3.0 + 0.1 * elev[0]);
1752         jb = np - ja + 6;
1753         for (i = ja - 1; i < jb; ++i)
1754                 zsys += elev[i];
1755         zsys /= (jb - ja + 1);
1756
1757         propv.mdvar = 12;
1758         qlrps(frq_mhz, zsys, q, pol, eps_dielect, sgm_conductivity, prop);
1759         qlrpfl(elev, propv.klim, propv.mdvar, prop, propa, propv);
1760         fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
1761         deltaH = prop.delta_h;
1762         q = prop.d - prop.d_L;
1763
1764         if (int(q) < 0.0)
1765                 strcpy(strmode, "Line-Of-Sight Mode");
1766         else {
1767                 if (int(q) == 0.0)
1768                         strcpy(strmode, "Single Horizon");
1769
1770                 else
1771                 if (int(q) > 0.0)
1772                         strcpy(strmode, "Double Horizon");
1773
1774                 if (prop.d <= prop.d_Ls || prop.d <= prop.dx)
1775                         strcat(strmode, ", Diffraction Dominant");
1776
1777                 else
1778                 if (prop.d > prop.dx)
1779                         strcat(strmode, ", Troposcatter Dominant");
1780         }
1781
1782         dbloss = avar(zr, 0.0, zc, prop, propv) + fs; //avar(time,location,confidence)
1783         errnum = prop.kwx;
1784 }
1785 #endif
1786
1787
1788
1789 //********************************************************
1790 //* Area Mode Calculations                               *
1791 //********************************************************
1792
1793
1794 #ifdef WITH_AREA
1795 static
1796 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)
1797 {
1798         // pol: 0-Horizontal, 1-Vertical
1799         // TSiteCriteria, RSiteCriteria:
1800         //                 0 - random, 1 - careful, 2 - very careful
1801         // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1802         //                4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1803         //                7-Maritime Temperate, Over Sea
1804         // ModVar: 0 - Single: pctConf is "Time/Situation/Location", pctTime, pctLoc not used
1805         //         1 - Individual: pctTime is "Situation/Location", pctConf is "Confidence", pctLoc not used
1806         //         2 - Mobile: pctTime is "Time/Locations (Reliability)", pctConf is "Confidence", pctLoc not used
1807         //         3 - Broadcast: pctTime is "Time", pctLoc is "Location", pctConf is "Confidence"
1808         // pctTime, pctLoc, pctConf: .01 to .99
1809         // errnum: 0- No Error.
1810         //         1- Warning: Some parameters are nearly out of range.
1811         //                     Results should be used with caution.
1812         //         2- Note: Default parameters have been substituted for impossible ones.
1813         //         3- Warning: A combination of parameters is out of range.
1814         //                     Results are probably invalid.
1815         //         Other-  Warning: Some parameters are out of range.
1816         //                          Results are probably invalid.
1817         // NOTE: strmode is not used at this time.
1818
1819         prop_type prop;
1820         propv_type propv;
1821         propa_type propa;
1822         double zt, zl, zc, xlb;
1823         double fs;
1824         long ivar;
1825         double eps, sgm;
1826         long ipol;
1827         int kst[2];
1828
1829         kst[0] = (int)TSiteCriteria;
1830         kst[1] = (int)RSiteCriteria;
1831         zt = qerfi(pctTime);
1832         zl = qerfi(pctLoc);
1833         zc = qerfi(pctConf);
1834         eps = eps_dielect;
1835         sgm = sgm_conductivity;
1836         prop.delta_h = deltaH;
1837         prop.h_g[0] = tht_m;
1838         prop.h_g[1] = rht_m;
1839         propv.klim = (long)radio_climate;
1840         prop.N_s = eno;
1841         prop.kwx = 0;
1842         ivar = (long)ModVar;
1843         ipol = (long)pol;
1844         qlrps(frq_mhz, 0.0, eno, ipol, eps, sgm, prop);
1845         qlra(kst, propv.klim, ivar, prop, propv);
1846
1847         if (propv.lvar < 1)
1848                 propv.lvar = 1;
1849
1850         lrprop(dist_km * 1000.0, prop, propa);
1851         fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
1852         xlb = fs + avar(zt, zl, zc, prop, propv);
1853         dbloss = xlb;
1854
1855         if (prop.kwx == 0)
1856                 errnum = 0;
1857         else
1858                 errnum = prop.kwx;
1859 }
1860 #endif
1861
1862 }