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