1 /*****************************************************************************\
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. *
12 *****************************************************************************
14 * Added comments that refer to itm.pdf and itmalg.pdf in a way to easy *
17 * // :0: Blah, page 0 This is the implementation of the code from *
18 * itm.pdf, section "0". The description is *
20 * [Alg 0.0] please refer to algorithm 0.0 from itm_alg.pdf *
22 * Holger Schurig, DH3HS *
23 \*****************************************************************************/
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
32 // Irregular Terrain Model (ITM) (Longley-Rice)
33 // *************************************
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
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
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
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
74 double aed; // diffraction coefficients
76 double aes; // scatter coefficients
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
86 int mymin(const int &i, const int &j)
96 int mymax(const int &i, const int &j)
106 double mymin(const double &a, const double &b)
116 double mymax(const double &a, const double &b)
126 double FORTRAN_DIM(const double &x, const double &y)
128 // This performs the FORTRAN DIM function.
129 // result is x-y if x is greater than y; otherwise result is 0.0
137 //#define set_warn(txt, err) printf("%s:%d %s\n", __func__, __LINE__, txt);
138 #define set_warn(txt, err) 1;
140 // :13: single-knife attenuation, page 6
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]
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
152 * Note that the arguments to this function aren't v, but v^2
155 double Fn(const double &v_square)
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;
163 a = 12.953 + 4.343 * log(v_square);
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].
175 double F(const double& x, const double& K)
179 // F = F_2(x, L), which is defined in [Alg 6.6]
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
189 fhtv = 17.372 * log(x) + fhtv;
191 // [Alg 6.6], lower part
192 fhtv = 2.5e-5 * x * x / K - 8.686 * w - 15.0;
195 // [Alg 6.3] and [Alg 6.4], lower part, which is G(x)
196 fhtv = 0.05751 * x - 4.343 * log(x);
198 // [Alg 6.4], middle part, but again XXX this doesn't match totally
200 w = 0.0134 * x * exp(-0.005 * x);
201 fhtv = (1.0 - w) * fhtv + w * (17.372 * log(x) - 117.0);
209 // :25: Tropospheric scatter frequency gain, [Alg 6.10ff], page 12
211 double H_0(double r, double et)
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};
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);
238 // XXX not sure what this means
240 h0fv = (1.0 - q) * h0fv + q * 4.343 * log((a[it] * x + b[it]) * x + 1.0);
246 // :25: This is the F(\Theta d) function for scatter fields, page 12
248 double F_0(double td)
251 if (td <= 10e3) // below 10 km
252 return 133.4 + 104.6 * td + 71.8 * log(td);
254 if (td <= 70e3) // between 10 km and 70 km
255 return 0.332e-3 + 0.212e-3 * td + 0.157e-3 * log(td);
257 return-4.343 + -1.086 * td + 2.171 * log(td);
261 // :10: Diffraction attenuation, page 4
263 double adiff(double s, prop_type &prop)
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.
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]
277 complex<double> prop_zgnd(prop.Z_g_real, prop.Z_g_imag);
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;
284 q += 10.0; // "C" from [Alg 4.9]
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));
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]
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];
306 alpha = pow(gamma_j_recip * prop.k, THIRD);
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]
313 // [Alg 4.20] ?, F(x, k) is in [Alg 6.4]
319 double gamma_o_recip, q, K, ds, theta, alpha, A_r, w;
321 // :12: Compute diffraction attenuation, page 5
323 theta = prop.theta_e + s * prop.gamma_e;
324 // XXX this is not [Alg 4.13]
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;
334 alpha = pow(gamma_o_recip * prop.k, THIRD);
335 // [Alg 4.17], note that qk is "1.0 / abs(prop_zgnd)" from above
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));
347 return (1.0 - w) * A_k + w * A_r + A_fo;
352 * :22: Scatter attenuation, page 9
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
360 * One needs to get TN101, especially chaper 9, to understand this function.
363 double A_scat(double s, prop_type &prop)
365 static double ad, rr, etq, h0s;
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];
377 etq = (5.67e-6 * prop.N_s - 2.32e-3) * prop.N_s + 0.031; // part of [Alg 4.67]
382 double h0, r1, r2, z0, ss, et, ett, theta_tick, q, temp;
384 // :24: Compute scatter attenuation, page 11
389 theta_tick = prop.theta_ej[0] + prop.theta_ej[1] + prop.gamma_e * s;
391 r2 = 2.0 * prop.k * theta_tick;
392 r1 = r2 * prop.h_e[0];
395 if (r1 < 0.2 && r2 < 0.2)
396 // The function is undefined
399 // XXX not like [Alg 4.65]
400 ss = (s - ad) / (s + ad);
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;
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;
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);
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));
420 if (h0 > 15.0 && h0s >= 0.0)
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)
431 - 0.1 * (prop.N_s - 301.0) * exp(-theta * s / D_0)
437 double abq_alos(complex<double> r)
439 return r.real() * r.real() + r.imag() * r.imag();
444 * :17: line-of-sight attenuation, page 8
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.
451 double A_los(double d, prop_type &prop)
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]
465 complex<double> prop_zgnd(prop.Z_g_real, prop.Z_g_imag);
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]
479 if (q < 0.25 || q < sps) // [Alg 4.48]
480 r = r * sqrt(sps / q);
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]
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;
489 // [Alg 4.51 and 4.44]
490 return (-4.343 * log(abq_alos(complex<double>(cos(q), -sin(q)) + r)) - alosv) * wls + alosv;
494 // :5: LRprop, page 2
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.
502 void lrprop(double d, prop_type &prop)
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;
515 // :6: Do secondary parameters, page 3
517 for (j = 0; j < 2; j++)
518 prop.d_Lsj[j] = sqrt(2.0 * prop.h_e[j] / prop.gamma_e);
521 prop.d_Ls = prop.d_Lsj[0] + prop.d_Lsj[1];
523 prop.d_L = prop.d_Lj[0] + prop.d_Lj[1];
525 prop.theta_e = mymax(prop.theta_ej[0] + prop.theta_ej[1], -prop.d_L * prop.gamma_e);
529 // :7: Check parameters range, page 3
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"
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);
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);
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);
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);
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);
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);
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);
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);
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);
590 dmin = abs(prop.h_e[0] - prop.h_e[1]) / 200e-3;
592 // :9: Diffraction coefficients, page 4
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.
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]
606 prop.emd = (a4 - a3) / (d4 - d3); // [Alg 4.8]
607 prop.aed = a3 - prop.emd * d3;
616 // :8: Check distance, page 3
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);
624 // Distance too small, use some indoor algorithm :-)
626 set_warn("Distance too small", 3);
627 prop.kwx = mymax(prop.kwx, 3);
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);
637 if (prop.d < prop.d_Ls) {
638 // :15: Line-of-sight calculations, page 7
640 // :16: Line-of-sight coefficients, page 7
641 q = A_los(0.0, prop);
643 a2 = prop.aed + d2 * prop.emd;
644 d0 = 1.908 * prop.k * prop.h_e[0] * prop.h_e[1]; // [Alg 4.38]
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]
650 d1 = mymax(-prop.aed / prop.emd, 0.25 * prop.d_L); // [Alg 4.39]
653 a1 = A_los(d1, prop); // [Alg 4.31]
657 a0 = A_los(d0, prop); // [Alg 4.30]
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;
663 prop.ak1 = (a2 - a0 - prop.ak2 * q) / (d2 - d0); // [Alg 4.33]
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]
669 if (prop.ak2 == 0.0) // [Alg 4.37]
676 prop.ak1 = FORTRAN_DIM(a2, a1) / (d2 - d1); // [Alg 4.40]
677 prop.ak2 = 0.0; // [Alg 4.41]
679 if (prop.ak1 == 0.0) // [Alg 4.37]
683 prop.ael = a2 - prop.ak1 * d2 - prop.ak2 * log(d2); // [Alg 4.42]
690 * The reference attenuation is determined as a function of the distance
691 * d from 3 piecewise formulatios. This is part one.
693 prop.A_ref = prop.ael + prop.ak1 * prop.d + prop.ak2 * log(prop.d);
696 if (prop.d <= 0.0 || prop.d >= prop.d_Ls) {
697 // :20: Troposcatter calculations, page 9
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]
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]
717 prop.dx = 10.e6; // [Alg 4.56]
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
726 prop.A_ref = prop.aed + prop.emd * prop.d; // diffraction region
729 prop.A_ref = mymax(prop.A_ref, 0.0);
733 /*****************************************************************************/
735 // :27: Variablility parameters, page 13
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
745 // 5: climate changed, or initialize everything
746 int mdvar; // desired mode of variability
747 int klim; // climate indicator
749 double sgc; // standard deviation of situation variability (confidence)
754 // :40: Prepare model for area prediction mode, page 21
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
766 void qlra(int kst[], int klimx, int mdvarx, prop_type &prop, propv_type &propv)
770 for (int j = 0; j < 2; ++j) {
773 prop.h_e[j] = prop.h_g[j];
780 if (prop.h_g[j] < 5.0)
781 q *= sin(0.3141593 * prop.h_g[j]);
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)));
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)));
791 prop.theta_ej[j] = (0.65 * prop.delta_h * (q / prop.d_Lj[j] - 1.0) - 2.0 * prop.h_e[j]) / q;
795 propv.lvar = mymax(propv.lvar, 3);
798 propv.mdvar = mdvarx;
799 propv.lvar = mymax(propv.lvar, 4);
810 // :51: Inverse of standard normal complementary probability
812 * The approximation is due to C. Hastings, Jr. ("Approximations for digital
813 * computers," Princeton Univ. Press, 1955) and the maximum error should be
817 double qerfi(double q)
820 const double c0 = 2.515516698;
821 const double c1 = 0.802853;
822 const double c2 = 0.010328;
823 const double d1 = 1.432788;
824 const double d2 = 0.189269;
825 const double d3 = 0.001308;
828 t = mymax(0.5 - fabs(x), 0.000001);
829 t = sqrt(-2.0 * log(t));
830 v = t - ((c2 * t + c1) * t + c0) / (((d3 * t + d2) * t + d1) * t + 1.0);
839 // :41: preparatory routine, page 20
841 * This subroutine converts
843 * surface refractivity reduced to sea level @en0
844 * general system elevation @zsys
845 * polarization and ground constants @eps, @sgm
848 * surface refractivity @ens
849 * effective earth curvature @gme
850 * surface impedance @zgnd
851 * It may be used with either the area prediction or the point-to-point mode.
854 void qlrps(double fmhz, double zsys, double en0, int ipol, double eps, double sgm, prop_type &prop)
857 const double Z_1 = 9.46e3; // 9.46 km from [Alg 1.2]
858 const double gma = 157e-9; // 157e-9 1/m from [Alg 1.3]
859 const double N_1 = 179.3; // 179.3 N-units from [Alg 1.3]
860 const double Z_0 = 376.62; // 376.62 Ohm from [Alg 1.5]
862 // Frequecy -> Wave number
863 prop.k = fmhz / f_0; // [Alg 1.1]
865 // Surface refractivity
868 prop.N_s *= exp(-zsys / Z_1); // [Alg 1.2]
870 // Earths effective curvature
871 prop.gamma_e = gma * (1.0 - 0.04665 * exp(prop.N_s / N_1)); // [Alg 1.3]
873 // Surface transfer impedance
874 complex<double> zq, prop_zgnd(prop.Z_g_real, prop.Z_g_imag);
875 zq = complex<double> (eps, Z_0 * sgm / prop.k); // [Alg 1.5]
876 prop_zgnd = sqrt(zq - 1.0);
878 // adjust surface transfer impedance for Polarization
880 prop_zgnd = prop_zgnd / zq; // [Alg 1.4]
882 prop.Z_g_real = prop_zgnd.real();
883 prop.Z_g_imag = prop_zgnd.imag();
887 // :30: Function curv, page 15
889 double curve(double const &c1, double const &c2, double const &x1, double const &x2, double const &x3, double const &de)
893 temp1 = (de - x2) / x3;
899 return (c1 + c2 / (1.0 + temp1)) * temp2 / (1.0 + temp2);
903 // :28: Area variablity, page 14
905 double avar(double zzt, double zzl, double zzc, prop_type &prop, propv_type &propv)
908 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;
910 // :29: Climatic constants, page 15
913 // 1: continental suptropical
914 // 2: maritime subtropical
916 // 4: continental temperature
917 // 5: maritime over land
918 // 6: maritime over sea
919 // equator contsup maritsup desert conttemp mariland marisea
920 const double bv1[7] = { -9.67, -0.62, 1.26, -9.21, -0.62, -0.39, 3.15};
921 const double bv2[7] = { 12.7, 9.19, 15.5, 9.05, 9.19, 2.86, 857.9};
922 const double xv1[7] = {144.9e3, 228.9e3, 262.6e3, 84.1e3, 228.9e3, 141.7e3, 2222.e3};
923 const double xv2[7] = {190.3e3, 205.2e3, 185.2e3, 101.1e3, 205.2e3, 315.9e3, 164.8e3};
924 const double xv3[7] = {133.8e3, 143.6e3, 99.8e3, 98.6e3, 143.6e3, 167.4e3, 116.3e3};
925 const double bsm1[7] = { 2.13, 2.66, 6.11, 1.98, 2.68, 6.86, 8.51};
926 const double bsm2[7] = { 159.5, 7.67, 6.65, 13.11, 7.16, 10.38, 169.8};
927 const double xsm1[7] = {762.2e3, 100.4e3, 138.2e3, 139.1e3, 93.7e3, 187.8e3, 609.8e3};
928 const double xsm2[7] = {123.6e3, 172.5e3, 242.2e3, 132.7e3, 186.8e3, 169.6e3, 119.9e3};
929 const double xsm3[7] = { 94.5e3, 136.4e3, 178.6e3, 193.5e3, 133.5e3, 108.9e3, 106.6e3};
930 const double bsp1[7] = { 2.11, 6.87, 10.08, 3.68, 4.75, 8.58, 8.43};
931 const double bsp2[7] = { 102.3, 15.53, 9.60, 159.3, 8.12, 13.97, 8.19};
932 const double xsp1[7] = {636.9e3, 138.7e3, 165.3e3, 464.4e3, 93.2e3, 216.0e3, 136.2e3};
933 const double xsp2[7] = {134.8e3, 143.7e3, 225.7e3, 93.1e3, 135.9e3, 152.0e3, 188.5e3};
934 const double xsp3[7] = { 95.6e3, 98.6e3, 129.7e3, 94.2e3, 113.4e3, 122.7e3, 122.9e3};
935 // bds1 -> is similar to C_D from table 5.1 at [Alg 5.8]
936 // bzd1 -> is similar to z_D from table 5.1 at [Alg 5.8]
937 const double bsd1[7] = { 1.224, 0.801, 1.380, 1.000, 1.224, 1.518, 1.518};
938 const double bzd1[7] = { 1.282, 2.161, 1.282, 20.0, 1.282, 1.282, 1.282};
939 const double bfm1[7] = { 1.0, 1.0, 1.0, 1.0, 0.92, 1.0, 1.0};
940 const double bfm2[7] = { 0.0, 0.0, 0.0, 0.0, 0.25, 0.0, 0.0};
941 const double bfm3[7] = { 0.0, 0.0, 0.0, 0.0, 1.77, 0.0, 0.0};
942 const double bfp1[7] = { 1.0, 0.93, 1.0, 0.93, 0.93, 1.0, 1.0};
943 const double bfp2[7] = { 0.0, 0.31, 0.0, 0.19, 0.31, 0.0, 0.0};
944 const double bfp3[7] = { 0.0, 2.00, 0.0, 1.79, 2.00, 0.0, 0.0};
945 const double rt = 7.8, rl = 24.0;
946 static bool no_location_variability, no_situation_variability;
947 double avarv, q, vs, zt, zl, zc;
951 if (propv.lvar > 0) {
952 // :31: Setup variablity constants, page 16
953 switch (propv.lvar) {
955 // Initialization or climate change
957 // if climate is wrong, use some "continental temperature" as default
958 // and set error indicator
959 if (propv.klim <= 0 || propv.klim > 7) {
961 prop.kwx = mymax(prop.kwx, 2);
962 set_warn("Climate index set to continental", 2);
965 // convert climate number into index into the climate tables
966 temp_klim = propv.klim - 1;
968 // :32: Climatic coefficients, page 17
969 cv1 = bv1[temp_klim];
970 cv2 = bv2[temp_klim];
971 yv1 = xv1[temp_klim];
972 yv2 = xv2[temp_klim];
973 yv3 = xv3[temp_klim];
974 csm1 = bsm1[temp_klim];
975 csm2 = bsm2[temp_klim];
976 ysm1 = xsm1[temp_klim];
977 ysm2 = xsm2[temp_klim];
978 ysm3 = xsm3[temp_klim];
979 csp1 = bsp1[temp_klim];
980 csp2 = bsp2[temp_klim];
981 ysp1 = xsp1[temp_klim];
982 ysp2 = xsp2[temp_klim];
983 ysp3 = xsp3[temp_klim];
984 csd1 = bsd1[temp_klim];
985 zd = bzd1[temp_klim];
986 cfm1 = bfm1[temp_klim];
987 cfm2 = bfm2[temp_klim];
988 cfm3 = bfm3[temp_klim];
989 cfp1 = bfp1[temp_klim];
990 cfp2 = bfp2[temp_klim];
991 cfp3 = bfp3[temp_klim];
995 // :33: Mode of variablity coefficients, page 17
997 // This code means that propv.mdvar can be
999 // 10 .. 13, then no_location_variability is set (no location variability)
1000 // 20 .. 23, then no_situation_variability is set (no situatian variability)
1001 // 30 .. 33, then no_location_variability and no_situation_variability is set
1003 no_situation_variability = kdv >= 20;
1004 if (no_situation_variability)
1007 no_location_variability = kdv >= 10;
1008 if (no_location_variability)
1011 if (kdv < 0 || kdv > 3) {
1013 set_warn("kdv set to 0", 2);
1014 prop.kwx = mymax(prop.kwx, 2);
1020 // :34: Frequency coefficients, page 18
1021 q = log(0.133 * prop.k);
1022 gm = cfm1 + cfm2 / ((cfm3 * q * cfm3 * q) + 1.0);
1023 gp = cfp1 + cfp2 / ((cfp3 * q * cfp3 * q) + 1.0);
1027 // :35: System coefficients, page 18
1028 // [Alg 5.3], effective distance
1030 const double a_1 = 9000e3; // 9000 km from [[Alg 5.3]
1031 //XXX I don't have any idea how they made up the third summand,
1032 //XXX text says a_1 * pow(k * D_1, -THIRD)
1033 //XXX with const double D_1 = 1266; // 1266 km
1034 dexa = sqrt(2*a_1 * prop.h_e[0]) +
1035 sqrt(2*a_1 * prop.h_e[1]) +
1036 pow((575.7e12 / prop.k), THIRD);
1041 // :36: Distance coefficients, page 18
1044 const double D_0 = 130e3; // 130 km from [Alg 5.4]
1046 de = D_0 * prop.d / dexa;
1048 de = D_0 + prop.d - dexa;
1052 * Quantiles of time variability are computed using a
1053 * variation of the methods described in Section 10 and
1054 * Annex III.7 of NBS~TN101, and also in CCIR
1055 * Report {238-3}. Those methods speak of eight or nine
1056 * discrete radio climates, of which seven have been
1057 * documented with corresponding empirical curves. It is
1058 * these empirical curves to which we refer below. They are
1059 * all curves of quantiles of deviations versus the
1060 * effective distance @de.
1062 vmd = curve(cv1, cv2, yv1, yv2, yv3, de); // [Alg 5.5]
1063 // [Alg 5.7], the slopes or "pseudo-standard deviations":
1064 // sgtm -> \sigma T-
1065 // sgtp -> \sigma T+
1066 sgtm = curve(csm1, csm2, ysm1, ysm2, ysm3, de) * gm;
1067 sgtp = curve(csp1, csp2, ysp1, ysp2, ysp3, de) * gp;
1068 // [Alg 5.8], ducting, "sgtd" -> \sigma TD
1070 tgtd = (sgtp - sgtd) * zd;
1072 // Location variability
1073 if (no_location_variability) {
1077 q = (1.0 - 0.8 * exp(-prop.d / 50e3)) * prop.delta_h * prop.k;
1079 sgl = 10.0 * q / (q + 13.0);
1082 // Situation variability
1083 if (no_situation_variability) {
1086 const double D = 100e3; // 100 km
1087 vs0 = (5.0 + 3.0 * exp(-de / D)); // [Alg 5.10]
1091 // Mark all constants as initialized
1096 // :37: Correct normal deviates, page 19
1100 // kdv is derived from prop.mdvar
1103 // Single message mode. Time, location and situation variability
1104 // are combined to form a confidence level.
1110 // Individual mode. Reliability is given by time
1111 // variability. Confidence. is a combination of location
1112 // and situation variability.
1117 // Mobile modes. Reliability is a combination of time and
1118 // location variability. Confidence. is given by the
1119 // situation variability.
1121 // case 3: Broadcast mode. like avar(zt, zl, zc).
1122 // Reliability is given by the two-fold statement of at
1123 // least qT of the time in qL of the locations. Confidence.
1124 // is given by the situation variability.
1126 if (fabs(zt) > 3.1 || fabs(zl) > 3.1 || fabs(zc) > 3.1) {
1127 set_warn("Situations variables not optimal", 1);
1128 prop.kwx = mymax(prop.kwx, 1);
1132 // :38: Resolve standard deviations, page 19
1139 sgt = sgtd + tgtd / zt;
1140 // [Alg 5.11], situation variability
1141 vs = vs0 + (sgt * zt * sgt * zt) / (rt + zc * zc) + (sgl * zl * sgl * zl) / (rl + zc * zc);
1144 // :39: Resolve deviations, page 19
1147 propv.sgc = sqrt(sgt * sgt + sgl * sgl + vs);
1151 propv.sgc = sqrt(sgl * sgl + vs);
1154 yr = sqrt(sgt * sgt + sgl * sgl) * zt;
1155 propv.sgc = sqrt(vs);
1157 yr = sgt * zt + sgl * zl;
1158 propv.sgc = sqrt(vs);
1161 // [Alg 5.1], area variability
1162 avarv = prop.A_ref - vmd - yr - propv.sgc * zc;
1166 avarv = avarv * (29.0 - avarv) / (29.0 - 10.0 * avarv);
1172 // :45: Find to horizons, page 24
1174 * Here we use the terrain profile @pfl to find the two horizons. Output
1175 * consists of the horizon distances @dl and the horizon take-off angles
1176 * @the. If the path is line-of-sight, the routine sets both horizon
1177 * distances equal to @dist.
1180 void hzns(double pfl[], prop_type &prop)
1184 double xi, za, zb, qc, q, sb, sa;
1188 za = pfl[2] + prop.h_g[0];
1189 zb = pfl[np+2] + prop.h_g[1];
1190 qc = 0.5 * prop.gamma_e;
1192 prop.theta_ej[1] = (zb - za) / prop.d;
1193 prop.theta_ej[0] = prop.theta_ej[1] - q;
1194 prop.theta_ej[1] = -prop.theta_ej[1] - q;
1195 prop.d_Lj[0] = prop.d;
1196 prop.d_Lj[1] = prop.d;
1203 for (int i = 1; i < np; i++) {
1206 q = pfl[i+2] - (qc * sa + prop.theta_ej[0]) * sa - za;
1209 prop.theta_ej[0] += q / sa;
1215 q = pfl[i+2] - (qc * sb + prop.theta_ej[1]) * sb - zb;
1218 prop.theta_ej[1] += q / sb;
1227 // :53: Linear least square fit, page 28
1229 void zlsq1(double z[], const double &x1, const double &x2, double& z0, double& zn)
1232 double xn, xa, xb, x, a, b;
1236 xa = int(FORTRAN_DIM(x1 / z[1], 0.0));
1237 xb = xn - int(FORTRAN_DIM(xn, x2 / z[1]));
1240 xa = FORTRAN_DIM(xa, 1.0);
1241 xb = xn - FORTRAN_DIM(xn, xb + 1.0);
1250 a = 0.5 * (z[ja+2] + z[jb+2]);
1251 b = 0.5 * (z[ja+2] - z[jb+2]) * x;
1253 for (int i = 2; i <= n; ++i) {
1261 b = b * 12.0 / ((xa * xa + 2.0) * xa);
1263 zn = a + b * (xn - xb);
1267 // :52: Provide a quantile and reorders array @a, page 27
1269 double qtile(const int &nn, double a[], const int &ir)
1271 double q = 0.0, r; /* Initializations by KD2BD */
1272 int m, n, i, j, j1 = 0, i0 = 0, k; /* Initializations by KD2BD */
1278 k = mymin(mymax(0, ir), n);
1289 while (i <= n && a[i] >= q)
1296 while (j >= m && a[j] <= q)
1331 // :50: Standard normal complementary probability, page 26
1333 double qerf(const double &z)
1335 const double b1 = 0.319381530, b2 = -0.356563782, b3 = 1.781477937;
1336 const double b4 = -1.821255987, b5 = 1.330274429;
1337 const double rp = 4.317008, rrt2pi = 0.398942280;
1347 qerfv = exp(-0.5 * x * x) * rrt2pi * ((((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t;
1351 qerfv = 1.0 - qerfv;
1358 // :48: Find interdecile range of elevations, page 25
1360 * Using the terrain profile @pfl we find \Delta h, the interdecile range of
1361 * elevations between the two points @x1 and @x2.
1364 double dlthx(double pfl[], const double &x1, const double &x2)
1366 int np, ka, kb, n, k, j;
1367 double dlthxv, sn, xa, xb;
1375 if (xb - xa < 2.0) // exit out
1378 ka = (int)(0.1 * (xb - xa + 8.0));
1379 ka = mymin(mymax(4, ka), 25);
1383 s = new double[n+2];
1386 xb = (xb - xa) / sn;
1387 k = (int)(xa + 1.0);
1390 for (j = 0; j < n; j++) {
1392 while (xa > 0.0 && k < np) {
1397 s[j+2] = pfl[k+2] + (pfl[k+2] - pfl[k+1]) * xa;
1401 zlsq1(s, 0.0, sn, xa, xb);
1402 xb = (xb - xa) / sn;
1404 for (j = 0; j < n; j++) {
1409 dlthxv = qtile(n - 1, s + 2, ka - 1) - qtile(n - 1, s + 2, kb - 1);
1410 dlthxv /= 1.0 - 0.8 * exp(-(x2 - x1) / 50.0e3);
1417 // :41: Prepare model for point-to-point operation, page 22
1419 * This mode requires the terrain profile lying between the terminals. This
1420 * should be a sequence of surface elevations at points along the great
1421 * circle path joining the two points. It should start at the ground beneath
1422 * the first terminal and end at the ground beneath the second. In the
1423 * present routine it is assumed that the elevations are equispaced along
1424 * the path. They are stored in the array @pfl along with two defining
1428 * pfl[0] = np, the number of points in the path
1429 * pfl[1] = xi, the length of each increment
1432 void qlrpfl(double pfl[], int klimx, int mdvarx, prop_type &prop, propv_type &propv)
1435 double xl[2], q, za, zb;
1437 prop.d = pfl[0] * pfl[1];
1440 // :44: determine horizons and dh from pfl, page 23
1442 for (j = 0; j < 2; j++)
1443 xl[j] = mymin(15.0 * prop.h_g[j], 0.1 * prop.d_Lj[j]);
1445 xl[1] = prop.d - xl[1];
1446 prop.delta_h = dlthx(pfl, xl[0], xl[1]);
1448 if (prop.d_Lj[0] + prop.d_Lj[1] > 1.5 * prop.d) {
1449 // :45: Redo line-of-sight horizons, page 23
1452 * If the path is line-of-sight, we still need to know where
1453 * the horizons might have been, and so we turn to
1454 * techniques used in area prediction mode.
1456 zlsq1(pfl, xl[0], xl[1], za, zb);
1457 prop.h_e[0] = prop.h_g[0] + FORTRAN_DIM(pfl[2], za);
1458 prop.h_e[1] = prop.h_g[1] + FORTRAN_DIM(pfl[np+2], zb);
1460 for (j = 0; j < 2; j++)
1461 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)));
1463 q = prop.d_Lj[0] + prop.d_Lj[1];
1466 q = ((prop.d / q) * (prop.d / q));
1468 for (j = 0; j < 2; j++) {
1470 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 for (j = 0; j < 2; j++) {
1475 q = sqrt(2.0 * prop.h_e[j] / prop.gamma_e);
1476 prop.theta_ej[j] = (0.65 * prop.delta_h * (q / prop.d_Lj[j] - 1.0) - 2.0 * prop.h_e[j]) / q;
1479 // :46: Transhorizon effective heights, page 23
1480 zlsq1(pfl, xl[0], 0.9 * prop.d_Lj[0], za, q);
1481 zlsq1(pfl, prop.d - 0.9 * prop.d_Lj[1], xl[1], q, zb);
1482 prop.h_e[0] = prop.h_g[0] + FORTRAN_DIM(pfl[2], za);
1483 prop.h_e[1] = prop.h_g[1] + FORTRAN_DIM(pfl[np+2], zb);
1487 propv.lvar = mymax(propv.lvar, 3);
1490 propv.mdvar = mdvarx;
1491 propv.lvar = mymax(propv.lvar, 4);
1503 //********************************************************
1504 //* Point-To-Point Mode Calculations *
1505 //********************************************************
1508 #ifdef WITH_POINT_TO_POINT
1511 void point_to_point(double elev[],
1512 double tht_m, // Transceiver above ground level
1513 double rht_m, // Receiver above groud level
1514 double eps_dielect, // Earth dielectric constant (rel. permittivity)
1515 double sgm_conductivity, // Earth conductivity (Siemens/m)
1516 double eno, // Atmospheric bending const, n-Units
1518 int radio_climate, // 1..7
1519 int pol, // 0 horizontal, 1 vertical
1520 double conf, // 0.01 .. .99, Fractions of situations
1521 double rel, // 0.01 .. .99, Fractions of time
1524 int &p_mode, // propagation mode selector
1525 double (&horizons)[2], // horizon distances
1528 // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1529 // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1530 // 7-Maritime Temperate, Over Sea
1531 // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
1532 // errnum: 0- No Error.
1533 // 1- Warning: Some parameters are nearly out of range.
1534 // Results should be used with caution.
1535 // 2- Note: Default parameters have been substituted for impossible ones.
1536 // 3- Warning: A combination of parameters is out of range.
1537 // Results are probably invalid.
1538 // Other- Warning: Some parameters are out of range.
1539 // Results are probably invalid.
1551 prop.h_g[0] = tht_m; // Tx height above ground level
1552 prop.h_g[1] = rht_m; // Rx height above ground level
1553 propv.klim = radio_climate;
1554 prop.kwx = 0; // no error yet
1555 propv.lvar = 5; // initialize all constants
1556 prop.mdp = -1; // point-to-point
1560 dkm = (elev[1] * elev[0]) / 1000.0;
1561 xkm = elev[1] / 1000.0;
1563 ja = (long)(3.0 + 0.1 * elev[0]);
1565 for (i = ja - 1; i < jb; ++i)
1567 zsys /= (jb - ja + 1);
1570 qlrps(frq_mhz, zsys, q, pol, eps_dielect, sgm_conductivity, prop);
1571 qlrpfl(elev, propv.klim, propv.mdvar, prop, propv);
1572 fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
1573 q = prop.d - prop.d_L;
1578 strcpy(strmode, "Line-Of-Sight Mode");
1581 if (int(q) == 0.0) {
1582 strcpy(strmode, "Single Horizon");
1583 horizons[0] = prop.d_Lj[0];
1589 strcpy(strmode, "Double Horizon");
1590 horizons[0] = prop.d_Lj[0];
1591 horizons[1] = prop.d_Lj[1];
1596 if (prop.d <= prop.d_Ls || prop.d <= prop.dx) {
1597 strcat(strmode, ", Diffraction Dominant");
1602 if (prop.d > prop.dx) {
1603 strcat(strmode, ", Troposcatter Dominant");
1609 dbloss = avar(zr, 0.0, zc, prop, propv) + fs;
1615 #ifdef WITH_POINT_TO_POINT_MDH
1617 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)
1619 // pol: 0-Horizontal, 1-Vertical
1620 // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1621 // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1622 // 7-Maritime Temperate, Over Sea
1623 // timepct, locpct, confpct: .01 to .99
1624 // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
1625 // propmode: Value Mode
1626 // -1 mode is undefined
1628 // 5 Single Horizon, Diffraction
1629 // 6 Single Horizon, Troposcatter
1630 // 9 Double Horizon, Diffraction
1631 // 10 Double Horizon, Troposcatter
1632 // errnum: 0- No Error.
1633 // 1- Warning: Some parameters are nearly out of range.
1634 // Results should be used with caution.
1635 // 2- Note: Default parameters have been substituted for impossible ones.
1636 // 3- Warning: A combination of parameters is out of range.
1637 // Results are probably invalid.
1638 // Other- Warning: Some parameters are out of range.
1639 // Results are probably invalid.
1645 double ztime, zloc, zconf;
1651 propmode = -1; // mode is undefined
1652 prop.h_g[0] = tht_m;
1653 prop.h_g[1] = rht_m;
1654 propv.klim = radio_climate;
1658 ztime = qerfi(timepct);
1659 zloc = qerfi(locpct);
1660 zconf = qerfi(confpct);
1663 dkm = (elev[1] * elev[0]) / 1000.0;
1664 xkm = elev[1] / 1000.0;
1666 ja = (long)(3.0 + 0.1 * elev[0]);
1668 for (i = ja - 1; i < jb; ++i)
1670 zsys /= (jb - ja + 1);
1673 qlrps(frq_mhz, zsys, q, pol, eps_dielect, sgm_conductivity, prop);
1674 qlrpfl(elev, propv.klim, propv.mdvar, prop, propa, propv);
1675 fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
1676 deltaH = prop.delta_h;
1677 q = prop.d - prop.d_L;
1680 propmode = 0; // Line-Of-Sight Mode
1683 propmode = 4; // Single Horizon
1686 propmode = 8; // Double Horizon
1688 if (prop.d <= prop.d_Ls || prop.d <= prop.dx)
1689 propmode += 1; // Diffraction Dominant
1692 if (prop.d > prop.dx)
1693 propmode += 2; // Troposcatter Dominant
1696 dbloss = avar(ztime, zloc, zconf, prop, propv) + fs; //avar(time,location,confidence)
1702 #ifdef WITH_POINT_TO_POINT_DH
1704 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)
1706 // pol: 0-Horizontal, 1-Vertical
1707 // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1708 // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1709 // 7-Maritime Temperate, Over Sea
1710 // conf, rel: .01 to .99
1711 // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
1712 // errnum: 0- No Error.
1713 // 1- Warning: Some parameters are nearly out of range.
1714 // Results should be used with caution.
1715 // 2- Note: Default parameters have been substituted for impossible ones.
1716 // 3- Warning: A combination of parameters is out of range.
1717 // Results are probably invalid.
1718 // Other- Warning: Some parameters are out of range.
1719 // Results are probably invalid.
1732 prop.h_g[0] = tht_m;
1733 prop.h_g[1] = rht_m;
1734 propv.klim = radio_climate;
1741 dkm = (elev[1] * elev[0]) / 1000.0;
1742 xkm = elev[1] / 1000.0;
1744 ja = (long)(3.0 + 0.1 * elev[0]);
1746 for (i = ja - 1; i < jb; ++i)
1748 zsys /= (jb - ja + 1);
1751 qlrps(frq_mhz, zsys, q, pol, eps_dielect, sgm_conductivity, prop);
1752 qlrpfl(elev, propv.klim, propv.mdvar, prop, propa, propv);
1753 fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
1754 deltaH = prop.delta_h;
1755 q = prop.d - prop.d_L;
1758 strcpy(strmode, "Line-Of-Sight Mode");
1761 strcpy(strmode, "Single Horizon");
1765 strcpy(strmode, "Double Horizon");
1767 if (prop.d <= prop.d_Ls || prop.d <= prop.dx)
1768 strcat(strmode, ", Diffraction Dominant");
1771 if (prop.d > prop.dx)
1772 strcat(strmode, ", Troposcatter Dominant");
1775 dbloss = avar(zr, 0.0, zc, prop, propv) + fs; //avar(time,location,confidence)
1782 //********************************************************
1783 //* Area Mode Calculations *
1784 //********************************************************
1789 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)
1791 // pol: 0-Horizontal, 1-Vertical
1792 // TSiteCriteria, RSiteCriteria:
1793 // 0 - random, 1 - careful, 2 - very careful
1794 // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1795 // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1796 // 7-Maritime Temperate, Over Sea
1797 // ModVar: 0 - Single: pctConf is "Time/Situation/Location", pctTime, pctLoc not used
1798 // 1 - Individual: pctTime is "Situation/Location", pctConf is "Confidence", pctLoc not used
1799 // 2 - Mobile: pctTime is "Time/Locations (Reliability)", pctConf is "Confidence", pctLoc not used
1800 // 3 - Broadcast: pctTime is "Time", pctLoc is "Location", pctConf is "Confidence"
1801 // pctTime, pctLoc, pctConf: .01 to .99
1802 // errnum: 0- No Error.
1803 // 1- Warning: Some parameters are nearly out of range.
1804 // Results should be used with caution.
1805 // 2- Note: Default parameters have been substituted for impossible ones.
1806 // 3- Warning: A combination of parameters is out of range.
1807 // Results are probably invalid.
1808 // Other- Warning: Some parameters are out of range.
1809 // Results are probably invalid.
1810 // NOTE: strmode is not used at this time.
1815 double zt, zl, zc, xlb;
1822 kst[0] = (int)TSiteCriteria;
1823 kst[1] = (int)RSiteCriteria;
1824 zt = qerfi(pctTime);
1826 zc = qerfi(pctConf);
1828 sgm = sgm_conductivity;
1829 prop.delta_h = deltaH;
1830 prop.h_g[0] = tht_m;
1831 prop.h_g[1] = rht_m;
1832 propv.klim = (long)radio_climate;
1835 ivar = (long)ModVar;
1837 qlrps(frq_mhz, 0.0, eno, ipol, eps, sgm, prop);
1838 qlra(kst, propv.klim, ivar, prop, propv);
1843 lrprop(dist_km * 1000.0, prop, propa);
1844 fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
1845 xlb = fs + avar(zt, zl, zc, prop, propv);