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 // *************************************
42 const float THIRD = (1.0/3.0);
43 const float f_0 = 47.7; // 47.7 MHz from [Alg 1.1], to convert frequency into wavenumber and vica versa
51 double d; // distance between the two terminals
52 double h_g[2]; // antenna structural heights (above ground)
53 double k; // wave number (= radio frequency)
54 double delta_h; // terrain irregularity parameter
55 double N_s; // minimum monthly surface refractivity, n-Units
56 double gamma_e; // earth's effective curvature
57 double Z_g_real; // surface transfer impedance of the ground
59 // Additional input for point-to-point mode
60 double h_e[2]; // antenna effective heights
61 double d_Lj[2]; // horizon distances
62 double theta_ej[2];// horizontal elevation angles
63 int mdp; // controlling mode: -1: point to point, 1 start of area, 0 area continuation
65 int kwx; // error indicator
66 double A_ref; // reference attenuation
67 // used to be propa_type, computed in lrprop()
68 double dx; // scatter distance
69 double ael; // line-of-sight coefficients
72 double aed; // diffraction coefficients
74 double aes; // scatter coefficients
76 double d_Lsj[2]; // smooth earth horizon distances
77 double d_Ls; // d_Lsj[] accumulated
78 double d_L; // d_Lj[] accumulated
79 double theta_e; // theta_ej[] accumulated, total bending angle
84 int mymin(const int &i, const int &j)
94 int mymax(const int &i, const int &j)
104 double mymin(const double &a, const double &b)
114 double mymax(const double &a, const double &b)
124 double FORTRAN_DIM(const double &x, const double &y)
126 // This performs the FORTRAN DIM function.
127 // result is x-y if x is greater than y; otherwise result is 0.0
135 #define set_warn(txt, err) printf("%s:%d %s\n", __func__, __LINE__, txt);
138 // :13: single-knife attenuation, page 6
140 * The attenuation due to a sigle knife edge -- this is an approximation of
141 * a Fresnel integral as a function of v^2. The non-approximated integral
142 * is documented as [Alg 4.21]
144 * Now, what is "single knife attenuation" ? Googling found some paper
145 * http://www.its.bldrdoc.gov/pub/ntia-rpt/81-86/81-86.pdf, which actually
146 * talks about multi-knife attenuation calculation. However, it says there
147 * that single-knife attenuation models attenuation over the edge of one
150 * Note that the arguments to this function aren't v, but v^2
153 double Fn(const double &v_square)
158 if (v_square <= 5.76) // this is the 2.40 from the text, but squared
159 a = 6.02 + 9.11 * sqrt(v_square) - 1.27 * v_square;
161 a = 12.953 + 4.343 * log(v_square);
169 * The heigh-gain over a smooth spherical earth -- to be used in the "three
170 * radii" mode. The approximation is that given in in [Alg 6.4ff].
173 double F(const double& x, const double& K)
177 // F = F_2(x, L), which is defined in [Alg 6.6]
181 // XXX the text says "or x * w^3 > 450"
182 if (K < 1e-5 || (x * w * w * w) > 5495.0) {
183 // F_2(x, k) = F_1(x), which is defined in [Alg 6.5]
184 // XXX but this isn't the same as in itm_alg.pdf
187 fhtv = 17.372 * log(x) + fhtv;
189 // [Alg 6.6], lower part
190 fhtv = 2.5e-5 * x * x / K - 8.686 * w - 15.0;
193 // [Alg 6.3] and [Alg 6.4], lower part, which is G(x)
194 fhtv = 0.05751 * x - 4.343 * log(x);
196 // [Alg 6.4], middle part, but again XXX this doesn't match totally
198 w = 0.0134 * x * exp(-0.005 * x);
199 fhtv = (1.0 - w) * fhtv + w * (17.372 * log(x) - 117.0);
207 // :25: Tropospheric scatter frequency gain, [Alg 6.10ff], page 12
209 double H_0(double r, double et)
211 // constants from [Alg 6.13]
212 const double a[5] = {25.0, 80.0, 177.0, 395.0, 705.0};
213 const double b[5] = {24.0, 45.0, 68.0, 80.0, 105.0};
233 // [Alg 6.13], calculates something like H_01(r,j), but not really XXX
234 h0fv = 4.343 * log((1.0 + a[it-1] * x + b[it-1]) * x);
236 // XXX not sure what this means
238 h0fv = (1.0 - q) * h0fv + q * 4.343 * log((a[it] * x + b[it]) * x + 1.0);
244 // :25: This is the F(\Theta d) function for scatter fields, page 12
246 double F_0(double td)
249 if (td <= 10e3) // below 10 km
250 return 133.4 + 104.6 * td + 71.8 * log(td);
252 if (td <= 70e3) // between 10 km and 70 km
253 return 0.332e-3 + 0.212e-3 * td + 0.157e-3 * log(td);
255 return-4.343 + -1.086 * td + 2.171 * log(td);
259 // :10: Diffraction attenuation, page 4
261 double adiff(double s, prop_type &prop)
264 * The function adiff finds the "Diffraction attenuation" at the
265 * distance s. It uses a convex combination of smooth earth
266 * diffraction and knife-edge diffraction.
268 static double wd1, xd1, A_fo, qk, aht, xht;
269 const double A = 151.03; // dimensionles constant from [Alg 4.20]
270 const double D = 50e3; // 50 km from [Alg 3.9], scale distance for \delta_h(s)
271 const double H = 16; // 16 m from [Alg 3.10]
272 const double ALPHA = 4.77e-4; // from [Alg 4.10]
275 complex<double> prop_zgnd(prop.Z_g_real, prop.Z_g_imag);
277 // :11: Prepare initial diffraction constants, page 5
278 double q = prop.h_g[0] * prop.h_g[1];
279 double qk = prop.h_e[0] * prop.h_e[1] - q;
282 q += 10.0; // "C" from [Alg 4.9]
284 // wd1 and xd1 are parts of Q in [Alg 4.10], but I cannot find this there
285 wd1 = sqrt(1.0 + qk / q);
286 xd1 = prop.d_L + prop.theta_e / prop.gamma_e; // [Alg 4.9] upper right
287 // \delta_h(s), [Alg 3.9]
288 q = (1.0 - 0.8 * exp(-prop.d_Ls / D)) * prop.delta_h;
289 // \sigma_h(s), [Alg 3.10]
290 q *= 0.78 * exp(-pow(q / H, 0.25));
292 // A_fo is the "clutter factor"
293 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]
294 // qk is part of the K_j calculation from [Alg 4.17]
295 qk = 1.0 / abs(prop_zgnd);
296 aht = 20.0; // 20 dB approximation for C_1(K) from [Alg 6.7], see also [Alg 4.25]
299 for (int j = 0; j < 2; ++j) {
300 double gamma_j_recip, alpha, K;
301 // [Alg 4.15], a is reciproke of gamma_j
302 gamma_j_recip = 0.5 * (prop.d_Lj[j] * prop.d_Lj[j]) / prop.h_e[j];
304 alpha = pow(gamma_j_recip * prop.k, THIRD);
307 // [Alg 4.18 and 6.2]
308 q = A * (1.607 - K) * alpha * prop.d_Lj[j] / gamma_j_recip;
309 // [Alg 4.19, high gain part]
311 // [Alg 4.20] ?, F(x, k) is in [Alg 6.4]
317 double gamma_o_recip, q, K, ds, theta, alpha, A_r, w;
319 // :12: Compute diffraction attenuation, page 5
321 theta = prop.theta_e + s * prop.gamma_e;
322 // XXX this is not [Alg 4.13]
324 q = 0.0795775 * prop.k * ds * theta * theta;
325 // [Alg 4.14], double knife edge attenuation
326 // Note that the arguments to Fn() are not v, but v^2
327 double A_k = Fn(q * prop.d_Lj[0] / (ds + prop.d_Lj[0])) +
328 Fn(q * prop.d_Lj[1] / (ds + prop.d_Lj[1]));
329 // kinda [Alg 4.15], just so that gamma_o is 1/a
330 gamma_o_recip = ds / theta;
332 alpha = pow(gamma_o_recip * prop.k, THIRD);
333 // [Alg 4.17], note that qk is "1.0 / abs(prop_zgnd)" from above
335 // [Alg 4.19], q is now X_o
336 q = A * (1.607 - K) * alpha * theta + xht;
337 // looks a bit like [Alg 4.20], rounded earth attenuation, or??
338 // note that G(x) should be "0.05751 * x - 10 * log(q)"
339 A_r = 0.05751 * q - 4.343 * log(q) - aht;
340 // I'm very unsure if this has anything to do with [Alg 4.9] or not
341 q = (wd1 + xd1 / s) * mymin(((1.0 - 0.8 * exp(-s / 50e3)) * prop.delta_h * prop.k), 6283.2);
342 // XXX this is NOT the same as the weighting factor from [Alg 4.9]
343 w = 25.1 / (25.1 + sqrt(q));
345 return (1.0 - w) * A_k + w * A_r + A_fo;
350 * :22: Scatter attenuation, page 9
352 * The function ascat finds the "scatter attenuation" at the distance d. It
353 * uses an approximation to the methods of NBS TN101 with check for
354 * inadmissable situations. For proper operation, the larger distance (d =
355 * d6) must be the first called. A call with d = 0 sets up initial
358 * One needs to get TN101, especially chaper 9, to understand this function.
361 double A_scat(double s, prop_type &prop)
363 static double ad, rr, etq, h0s;
366 // :23: Prepare initial scatter constants, page 10
367 ad = prop.d_Lj[0] - prop.d_Lj[1];
368 rr = prop.h_e[1] / prop.h_e[0];
375 etq = (5.67e-6 * prop.N_s - 2.32e-3) * prop.N_s + 0.031; // part of [Alg 4.67]
380 double h0, r1, r2, z0, ss, et, ett, theta_tick, q, temp;
382 // :24: Compute scatter attenuation, page 11
387 theta_tick = prop.theta_ej[0] + prop.theta_ej[1] + prop.gamma_e * s;
389 r2 = 2.0 * prop.k * theta_tick;
390 r1 = r2 * prop.h_e[0];
393 if (r1 < 0.2 && r2 < 0.2)
394 // The function is undefined
397 // XXX not like [Alg 4.65]
398 ss = (s - ad) / (s + ad);
401 q = mymin(mymax(0.1, q), 10.0);
402 // XXX not like [Alg 4.66]
403 z0 = (s - ad) * (s + ad) * theta_tick * 0.25 / s;
405 temp = mymin(1.7, z0 / 8.0e3);
406 temp = temp * temp * temp * temp * temp * temp;
407 et = (etq * exp(-temp) + 1.0) * z0 / 1.7556e3;
409 ett = mymax(et, 1.0);
410 h0 = (H_0(r1, ett) + H_0(r2, ett)) * 0.5; // [Alg 6.12]
411 h0 += mymin(h0, (1.38 - log(ett)) * log(ss) * log(q) * 0.49); // [Alg 6.10 and 6.11]
412 h0 = FORTRAN_DIM(h0, 0.0);
416 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));
418 if (h0 > 15.0 && h0s >= 0.0)
423 double theta = prop.theta_e + s * prop.gamma_e; // [Alg 4.60]
424 // [Alg 4.63 and 6.8]
425 const double D_0 = 40e3; // 40 km from [Alg 6.8]
426 const double H = 47.7; // 47.7 m from [Alg 4.63]
427 return 4.343 * log(prop.k * H * theta * theta * theta * theta)
429 - 0.1 * (prop.N_s - 301.0) * exp(-theta * s / D_0)
435 double abq_alos(complex<double> r)
437 return r.real() * r.real() + r.imag() * r.imag();
442 * :17: line-of-sight attenuation, page 8
444 * The function alos finds the "line-of-sight attenuation" at the distance
445 * d. It uses a convex combination of plane earth fields and diffracted
446 * fields. A call with d=0 sets up initial constants.
449 double A_los(double d, prop_type &prop)
454 // :18: prepare initial line-of-sight constants, page 8
455 const double D1 = 47.7; // 47.7 m from [Alg 4.43]
456 const double D2 = 10e3; // 10 km from [Alg 4.43]
457 const double D1R = 1 / D1;
458 // weighting factor w
459 wls = D1R / (D1R + prop.k * prop.delta_h / mymax(D2, prop.d_Ls)); // [Alg 4.43]
463 complex<double> prop_zgnd(prop.Z_g_real, prop.Z_g_imag);
467 // :19: compute line of sight attentuation, page 8
468 const double D = 50e3; // 50 km from [Alg 3.9]
469 const double H = 16.0; // 16 m from [Alg 3.10]
470 q = (1.0 - 0.8 * exp(-d / D)) * prop.delta_h; // \Delta h(d), [Alg 3.9]
471 s = 0.78 * q * exp(-pow(q / H, 0.25)); // \sigma_h(d), [Alg 3.10]
472 q = prop.h_e[0] + prop.h_e[1];
473 sps = q / sqrt(d * d + q * q); // sin(\psi), [Alg 4.46]
474 r = (sps - prop_zgnd) / (sps + prop_zgnd) * exp(-mymin(10.0, prop.k * s * sps)); // [Alg 4.47]
477 if (q < 0.25 || q < sps) // [Alg 4.48]
478 r = r * sqrt(sps / q);
480 double alosv = prop.emd * d + prop.aed; // [Alg 4.45]
481 q = prop.k * prop.h_e[0] * prop.h_e[1] * 2.0 / d; // [Alg 4.49]
483 // M_PI is pi, M_PI_2 is pi/2
484 if (q > M_PI_2) // [Alg 4.50]
485 q = M_PI - (M_PI_2 * M_PI_2) / q;
487 // [Alg 4.51 and 4.44]
488 return (-4.343 * log(abq_alos(complex<double>(cos(q), -sin(q)) + r)) - alosv) * wls + alosv;
492 // :5: LRprop, page 2
494 * The value mdp controls some of the program flow. When it equals -1 we are
495 * in point-to-point mode, when 1 we are beginning the area mode, and when 0
496 * we are continuing the area mode. The assumption is that when one uses the
497 * area mode, one will want a sequence of results for varying distances.
500 void lrprop(double d, prop_type &prop)
502 static bool wlos, wscat;
503 static double dmin, xae;
504 complex<double> prop_zgnd(prop.Z_g_real, prop.Z_g_imag);
505 double a0, a1, a2, a3, a4, a5, a6;
506 double d0, d1, d2, d3, d4, d5, d6;
513 // :6: Do secondary parameters, page 3
515 for (j = 0; j < 2; j++)
516 prop.d_Lsj[j] = sqrt(2.0 * prop.h_e[j] / prop.gamma_e);
519 prop.d_Ls = prop.d_Lsj[0] + prop.d_Lsj[1];
521 prop.d_L = prop.d_Lj[0] + prop.d_Lj[1];
523 prop.theta_e = mymax(prop.theta_ej[0] + prop.theta_ej[1], -prop.d_L * prop.gamma_e);
527 // :7: Check parameters range, page 3
529 // kwx is some kind of error indicator. Setting kwx to 1
530 // means "results are slightly bad", while setting it to 4
531 // means "my calculations will be bogus"
533 // Frequency must be between 40 MHz and 10 GHz
534 // Wave number (wn) = 2 * M_PI / lambda, or f/f0, where f0 = 47.7 MHz*m;
535 // 0.838 => 40 MHz, 210 => 10 GHz
536 if (prop.k < 0.838 || prop.k > 210.0) {
537 set_warn("Frequency not optimal", 1);
538 prop.kwx = mymax(prop.kwx, 1);
541 // Surface refractivity, see [Alg 1.2]
542 if (prop.N_s < 250.0 || prop.N_s > 400.0) {
543 set_warn("Surface refractivity out-of-bounds", 4);
546 // Earth's effective curvature, see [Alg 1.3]
547 if (prop.gamma_e < 75e-9 || prop.gamma_e > 250e-9) {
548 set_warn("Earth's curvature out-of-bounds", 4);
551 // Surface transfer impedance, see [Alg 1.4]
552 if (prop_zgnd.real() <= abs(prop_zgnd.imag())) {
553 set_warn("Surface transfer impedance out-of-bounds", 4);
556 // Calulating outside of 20 MHz to 40 GHz is really bad
557 if (prop.k < 0.419 || prop.k > 420.0) {
558 set_warn("Frequency out-of-bounds", 4);
561 for (j = 0; j < 2; j++) {
562 // Antenna structural height should be between 1 and 1000 m
563 if (prop.h_g[j] < 1.0 || prop.h_g[j] > 1000.0) {
564 set_warn("Antenna height not optimal", 1);
565 prop.kwx = mymax(prop.kwx, 1);
568 // Horizon elevation angle
569 if (abs(prop.theta_ej[j]) > 200e-3) {
570 set_warn("Horizon elevation weird", 3);
571 prop.kwx = mymax(prop.kwx, 3);
574 // Horizon distance dl,
575 // Smooth earth horizon distance dls
576 if (prop.d_Lj[j] < 0.1 * prop.d_Lsj[j] ||
577 prop.d_Lj[j] > 3.0 * prop.d_Lsj[j]) {
578 set_warn("Horizon distance weird", 3);
579 prop.kwx = mymax(prop.kwx, 3);
581 // Antenna structural height must between 0.5 and 3000 m
582 if (prop.h_g[j] < 0.5 || prop.h_g[j] > 3000.0) {
583 set_warn("Antenna heights out-of-bounds", 4);
588 dmin = abs(prop.h_e[0] - prop.h_e[1]) / 200e-3;
590 // :9: Diffraction coefficients, page 4
592 * This is the region beyond the smooth-earth horizon at
593 * D_Lsa and short of where isotropic scatter takes over. It
594 * is a key to central region and associated coefficients
595 * must always be computed.
597 q = adiff(0.0, prop);
598 xae = pow(prop.k * prop.gamma_e * prop.gamma_e, -THIRD); // [Alg 4.2]
599 d3 = mymax(prop.d_Ls, 1.3787 * xae + prop.d_L); // [Alg 4.3]
600 d4 = d3 + 2.7574 * xae; // [Alg 4.4]
601 a3 = adiff(d3, prop); // [Alg 4.5]
602 a4 = adiff(d4, prop); // [Alg 4.7]
604 prop.emd = (a4 - a3) / (d4 - d3); // [Alg 4.8]
605 prop.aed = a3 - prop.emd * d3;
614 // :8: Check distance, page 3
616 // Distance above 1000 km is guessing
617 if (prop.d > 1000e3) {
618 set_warn("Distance not optimal", 1);
619 prop.kwx = mymax(prop.kwx, 1);
622 // Distance too small, use some indoor algorithm :-)
624 set_warn("Distance too small", 3);
625 prop.kwx = mymax(prop.kwx, 3);
628 // Distance above 2000 km is really bad, don't do that
629 if (prop.d < 1e3 || prop.d > 2000e3) {
630 set_warn("Distance out-of-bounds", 4);
635 if (prop.d < prop.d_Ls) {
636 // :15: Line-of-sight calculations, page 7
638 // :16: Line-of-sight coefficients, page 7
639 q = A_los(0.0, prop);
641 a2 = prop.aed + d2 * prop.emd;
642 d0 = 1.908 * prop.k * prop.h_e[0] * prop.h_e[1]; // [Alg 4.38]
644 if (prop.aed >= 0.0) {
645 d0 = mymin(d0, 0.5 * prop.d_L); // [Alg 4.28]
646 d1 = d0 + 0.25 * (prop.d_L - d0); // [Alg 4.29]
648 d1 = mymax(-prop.aed / prop.emd, 0.25 * prop.d_L); // [Alg 4.39]
651 a1 = A_los(d1, prop); // [Alg 4.31]
655 a0 = A_los(d0, prop); // [Alg 4.30]
657 prop.ak2 = mymax(0.0, ((d2 - d0) * (a1 - a0) - (d1 - d0) * (a2 - a0)) / ((d2 - d0) * log(d1 / d0) - (d1 - d0) * q)); // [Alg 4.32]
658 wq = prop.aed >= 0.0 || prop.ak2 > 0.0;
661 prop.ak1 = (a2 - a0 - prop.ak2 * q) / (d2 - d0); // [Alg 4.33]
663 if (prop.ak1 < 0.0) {
664 prop.ak1 = 0.0; // [Alg 4.36]
665 prop.ak2 = FORTRAN_DIM(a2, a0) / q; // [Alg 4.35]
667 if (prop.ak2 == 0.0) // [Alg 4.37]
674 prop.ak1 = FORTRAN_DIM(a2, a1) / (d2 - d1); // [Alg 4.40]
675 prop.ak2 = 0.0; // [Alg 4.41]
677 if (prop.ak1 == 0.0) // [Alg 4.37]
681 prop.ael = a2 - prop.ak1 * d2 - prop.ak2 * log(d2); // [Alg 4.42]
688 * The reference attenuation is determined as a function of the distance
689 * d from 3 piecewise formulatios. This is part one.
691 prop.A_ref = prop.ael + prop.ak1 * prop.d + prop.ak2 * log(prop.d);
694 if (prop.d <= 0.0 || prop.d >= prop.d_Ls) {
695 // :20: Troposcatter calculations, page 9
697 // :21: Troposcatter coefficients
698 const double DS = 200e3; // 200 km from [Alg 4.53]
699 const double HS = 47.7; // 47.7 m from [Alg 4.59]
700 q = A_scat(0.0, prop);
701 d5 = prop.d_L + DS; // [Alg 4.52]
702 d6 = d5 + DS; // [Alg 4.53]
703 a6 = A_scat(d6, prop); // [Alg 4.54]
704 a5 = A_scat(d5, prop); // [Alg 4.55]
707 prop.ems = (a6 - a5) / DS; // [Alg 4.57]
708 prop.dx = mymax(prop.d_Ls, // [Alg 4.58]
709 mymax(prop.d_L + 0.3 * xae * log(HS * prop.k),
710 (a5 - prop.aed - prop.ems * d5) / (prop.emd - prop.ems)));
711 prop.aes = (prop.emd - prop.ems) * prop.dx + prop.aed; // [Alg 4.59]
715 prop.dx = 10.e6; // [Alg 4.56]
720 // [Alg 4.1], part two and three.
721 if (prop.d > prop.dx)
722 prop.A_ref = prop.aes + prop.ems * prop.d; // scatter region
724 prop.A_ref = prop.aed + prop.emd * prop.d; // diffraction region
727 prop.A_ref = mymax(prop.A_ref, 0.0);
731 /*****************************************************************************/
733 // :27: Variablility parameters, page 13
737 int lvar; // control switch for initialisation and/or recalculation
738 // 0: no recalculations needed
739 // 1: distance changed
740 // 2: antenna heights changed
741 // 3: frequency changed
743 // 5: climate changed, or initialize everything
744 int mdvar; // desired mode of variability
745 int klim; // climate indicator
747 double sgc; // standard deviation of situation variability (confidence)
752 // :40: Prepare model for area prediction mode, page 21
754 * This is used to prepare the model in the area prediction mode. Normally,
755 * one first calls qlrps() and then qlra(). Before calling the latter, one
756 * should have defined in prop_type the antenna heights @hg, the terrain
757 * irregularity parameter @dh , and (probably through qlrps() the variables
758 * @wn, @ens, @gme, and @zgnd. The input @kst will define siting criteria
759 * for the terminals, @klimx the climate, and @mdvarx the mode of
760 * variability. If @klimx <= 0 or @mdvarx < 0 the associated parameters
764 void qlra(int kst[], int klimx, int mdvarx, prop_type &prop, propv_type &propv)
768 for (int j = 0; j < 2; ++j) {
771 prop.h_e[j] = prop.h_g[j];
778 if (prop.h_g[j] < 5.0)
779 q *= sin(0.3141593 * prop.h_g[j]);
781 prop.h_e[j] = prop.h_g[j] + (1.0 + q) * exp(-mymin(20.0, 2.0 * prop.h_g[j] / mymax(1e-3, prop.delta_h)));
784 // [Alg 3.3], upper function. q is d_Ls_j
785 const double H_3 = 5; // 5m from [Alg 3.3]
786 q = sqrt(2.0 * prop.h_e[j] / prop.gamma_e);
787 prop.d_Lj[j] = q * exp(-0.07 * sqrt(prop.delta_h / mymax(prop.h_e[j], H_3)));
789 prop.theta_ej[j] = (0.65 * prop.delta_h * (q / prop.d_Lj[j] - 1.0) - 2.0 * prop.h_e[j]) / q;
793 propv.lvar = mymax(propv.lvar, 3);
796 propv.mdvar = mdvarx;
797 propv.lvar = mymax(propv.lvar, 4);
808 // :51: Inverse of standard normal complementary probability
810 * The approximation is due to C. Hastings, Jr. ("Approximations for digital
811 * computers," Princeton Univ. Press, 1955) and the maximum error should be
815 double qerfi(double q)
818 const double c0 = 2.515516698;
819 const double c1 = 0.802853;
820 const double c2 = 0.010328;
821 const double d1 = 1.432788;
822 const double d2 = 0.189269;
823 const double d3 = 0.001308;
826 t = mymax(0.5 - fabs(x), 0.000001);
827 t = sqrt(-2.0 * log(t));
828 v = t - ((c2 * t + c1) * t + c0) / (((d3 * t + d2) * t + d1) * t + 1.0);
837 // :41: preparatory routine, page 20
839 * This subroutine converts
841 * surface refractivity reduced to sea level @en0
842 * general system elevation @zsys
843 * polarization and ground constants @eps, @sgm
846 * surface refractivity @ens
847 * effective earth curvature @gme
848 * surface impedance @zgnd
849 * It may be used with either the area prediction or the point-to-point mode.
852 void qlrps(double fmhz, double zsys, double en0, int ipol, double eps, double sgm, prop_type &prop)
855 const double Z_1 = 9.46e3; // 9.46 km from [Alg 1.2]
856 const double gma = 157e-9; // 157e-9 1/m from [Alg 1.3]
857 const double N_1 = 179.3; // 179.3 N-units from [Alg 1.3]
858 const double Z_0 = 376.62; // 376.62 Ohm from [Alg 1.5]
860 // Frequecy -> Wave number
861 prop.k = fmhz / f_0; // [Alg 1.1]
863 // Surface refractivity
866 prop.N_s *= exp(-zsys / Z_1); // [Alg 1.2]
868 // Earths effective curvature
869 prop.gamma_e = gma * (1.0 - 0.04665 * exp(prop.N_s / N_1)); // [Alg 1.3]
871 // Surface transfer impedance
872 complex<double> zq, prop_zgnd(prop.Z_g_real, prop.Z_g_imag);
873 zq = complex<double> (eps, Z_0 * sgm / prop.k); // [Alg 1.5]
874 prop_zgnd = sqrt(zq - 1.0);
876 // adjust surface transfer impedance for Polarization
878 prop_zgnd = prop_zgnd / zq; // [Alg 1.4]
880 prop.Z_g_real = prop_zgnd.real();
881 prop.Z_g_imag = prop_zgnd.imag();
885 // :30: Function curv, page 15
887 double curve(double const &c1, double const &c2, double const &x1, double const &x2, double const &x3, double const &de)
891 temp1 = (de - x2) / x3;
897 return (c1 + c2 / (1.0 + temp1)) * temp2 / (1.0 + temp2);
901 // :28: Area variablity, page 14
903 double avar(double zzt, double zzl, double zzc, prop_type &prop, propv_type &propv)
906 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;
908 // :29: Climatic constants, page 15
911 // 1: continental suptropical
912 // 2: maritime subtropical
914 // 4: continental temperature
915 // 5: maritime over land
916 // 6: maritime over sea
917 // equator contsup maritsup desert conttemp mariland marisea
918 const double bv1[7] = { -9.67, -0.62, 1.26, -9.21, -0.62, -0.39, 3.15};
919 const double bv2[7] = { 12.7, 9.19, 15.5, 9.05, 9.19, 2.86, 857.9};
920 const double xv1[7] = {144.9e3, 228.9e3, 262.6e3, 84.1e3, 228.9e3, 141.7e3, 2222.e3};
921 const double xv2[7] = {190.3e3, 205.2e3, 185.2e3, 101.1e3, 205.2e3, 315.9e3, 164.8e3};
922 const double xv3[7] = {133.8e3, 143.6e3, 99.8e3, 98.6e3, 143.6e3, 167.4e3, 116.3e3};
923 const double bsm1[7] = { 2.13, 2.66, 6.11, 1.98, 2.68, 6.86, 8.51};
924 const double bsm2[7] = { 159.5, 7.67, 6.65, 13.11, 7.16, 10.38, 169.8};
925 const double xsm1[7] = {762.2e3, 100.4e3, 138.2e3, 139.1e3, 93.7e3, 187.8e3, 609.8e3};
926 const double xsm2[7] = {123.6e3, 172.5e3, 242.2e3, 132.7e3, 186.8e3, 169.6e3, 119.9e3};
927 const double xsm3[7] = { 94.5e3, 136.4e3, 178.6e3, 193.5e3, 133.5e3, 108.9e3, 106.6e3};
928 const double bsp1[7] = { 2.11, 6.87, 10.08, 3.68, 4.75, 8.58, 8.43};
929 const double bsp2[7] = { 102.3, 15.53, 9.60, 159.3, 8.12, 13.97, 8.19};
930 const double xsp1[7] = {636.9e3, 138.7e3, 165.3e3, 464.4e3, 93.2e3, 216.0e3, 136.2e3};
931 const double xsp2[7] = {134.8e3, 143.7e3, 225.7e3, 93.1e3, 135.9e3, 152.0e3, 188.5e3};
932 const double xsp3[7] = { 95.6e3, 98.6e3, 129.7e3, 94.2e3, 113.4e3, 122.7e3, 122.9e3};
933 // bds1 -> is similar to C_D from table 5.1 at [Alg 5.8]
934 // bzd1 -> is similar to z_D from table 5.1 at [Alg 5.8]
935 const double bsd1[7] = { 1.224, 0.801, 1.380, 1.000, 1.224, 1.518, 1.518};
936 const double bzd1[7] = { 1.282, 2.161, 1.282, 20.0, 1.282, 1.282, 1.282};
937 const double bfm1[7] = { 1.0, 1.0, 1.0, 1.0, 0.92, 1.0, 1.0};
938 const double bfm2[7] = { 0.0, 0.0, 0.0, 0.0, 0.25, 0.0, 0.0};
939 const double bfm3[7] = { 0.0, 0.0, 0.0, 0.0, 1.77, 0.0, 0.0};
940 const double bfp1[7] = { 1.0, 0.93, 1.0, 0.93, 0.93, 1.0, 1.0};
941 const double bfp2[7] = { 0.0, 0.31, 0.0, 0.19, 0.31, 0.0, 0.0};
942 const double bfp3[7] = { 0.0, 2.00, 0.0, 1.79, 2.00, 0.0, 0.0};
943 const double rt = 7.8, rl = 24.0;
944 static bool no_location_variability, no_situation_variability;
945 double avarv, q, vs, zt, zl, zc;
949 if (propv.lvar > 0) {
950 // :31: Setup variablity constants, page 16
951 switch (propv.lvar) {
953 // Initialization or climate change
955 // if climate is wrong, use some "continental temperature" as default
956 // and set error indicator
957 if (propv.klim <= 0 || propv.klim > 7) {
959 prop.kwx = mymax(prop.kwx, 2);
960 set_warn("Climate index set to continental", 2);
963 // convert climate number into index into the climate tables
964 temp_klim = propv.klim - 1;
966 // :32: Climatic coefficients, page 17
967 cv1 = bv1[temp_klim];
968 cv2 = bv2[temp_klim];
969 yv1 = xv1[temp_klim];
970 yv2 = xv2[temp_klim];
971 yv3 = xv3[temp_klim];
972 csm1 = bsm1[temp_klim];
973 csm2 = bsm2[temp_klim];
974 ysm1 = xsm1[temp_klim];
975 ysm2 = xsm2[temp_klim];
976 ysm3 = xsm3[temp_klim];
977 csp1 = bsp1[temp_klim];
978 csp2 = bsp2[temp_klim];
979 ysp1 = xsp1[temp_klim];
980 ysp2 = xsp2[temp_klim];
981 ysp3 = xsp3[temp_klim];
982 csd1 = bsd1[temp_klim];
983 zd = bzd1[temp_klim];
984 cfm1 = bfm1[temp_klim];
985 cfm2 = bfm2[temp_klim];
986 cfm3 = bfm3[temp_klim];
987 cfp1 = bfp1[temp_klim];
988 cfp2 = bfp2[temp_klim];
989 cfp3 = bfp3[temp_klim];
993 // :33: Mode of variablity coefficients, page 17
995 // This code means that propv.mdvar can be
997 // 10 .. 13, then no_location_variability is set (no location variability)
998 // 20 .. 23, then no_situation_variability is set (no situatian variability)
999 // 30 .. 33, then no_location_variability and no_situation_variability is set
1001 no_situation_variability = kdv >= 20;
1002 if (no_situation_variability)
1003 no_situation_variability -= 20;
1005 no_location_variability = kdv >= 10;
1006 if (no_location_variability)
1009 if (kdv < 0 || kdv > 3) {
1011 set_warn("kdv set to 0", 2);
1012 prop.kwx = mymax(prop.kwx, 2);
1018 // :34: Frequency coefficients, page 18
1019 q = log(0.133 * prop.k);
1020 gm = cfm1 + cfm2 / ((cfm3 * q * cfm3 * q) + 1.0);
1021 gp = cfp1 + cfp2 / ((cfp3 * q * cfp3 * q) + 1.0);
1025 // :35: System coefficients, page 18
1026 // [Alg 5.3], effective distance
1028 const double a_1 = 9000e3; // 9000 km from [[Alg 5.3]
1029 //XXX I don't have any idea how they made up the third summand,
1030 //XXX text says a_1 * pow(k * D_1, -THIRD)
1031 //XXX with const double D_1 = 1266; // 1266 km
1032 dexa = sqrt(2*a_1 * prop.h_e[0]) +
1033 sqrt(2*a_1 * prop.h_e[1]) +
1034 pow((575.7e12 / prop.k), THIRD);
1039 // :36: Distance coefficients, page 18
1042 const double D_0 = 130e3; // 130 km from [Alg 5.4]
1044 de = D_0 * prop.d / dexa;
1046 de = D_0 + prop.d - dexa;
1050 * Quantiles of time variability are computed using a
1051 * variation of the methods described in Section 10 and
1052 * Annex III.7 of NBS~TN101, and also in CCIR
1053 * Report {238-3}. Those methods speak of eight or nine
1054 * discrete radio climates, of which seven have been
1055 * documented with corresponding empirical curves. It is
1056 * these empirical curves to which we refer below. They are
1057 * all curves of quantiles of deviations versus the
1058 * effective distance @de.
1060 vmd = curve(cv1, cv2, yv1, yv2, yv3, de); // [Alg 5.5]
1061 // [Alg 5.7], the slopes or "pseudo-standard deviations":
1062 // sgtm -> \sigma T-
1063 // sgtp -> \sigma T+
1064 sgtm = curve(csm1, csm2, ysm1, ysm2, ysm3, de) * gm;
1065 sgtp = curve(csp1, csp2, ysp1, ysp2, ysp3, de) * gp;
1066 // [Alg 5.8], ducting, "sgtd" -> \sigma TD
1068 tgtd = (sgtp - sgtd) * zd;
1070 // Location variability
1071 if (no_location_variability) {
1075 q = (1.0 - 0.8 * exp(-prop.d / 50e3)) * prop.delta_h * prop.k;
1077 sgl = 10.0 * q / (q + 13.0);
1080 // Situation variability
1081 if (no_situation_variability) {
1084 const double D = 100e3; // 100 km
1085 vs0 = (5.0 + 3.0 * exp(-de / D)); // [Alg 5.10]
1089 // Mark all constants as initialized
1094 // :37: Correct normal deviates, page 19
1098 // kdv is derived from prop.mdvar
1101 // Single message mode. Time, location and situation variability
1102 // are combined to form a confidence level.
1108 // Individual mode. Reliability is given by time
1109 // variability. Confidence. is a combination of location
1110 // and situation variability.
1115 // Mobile modes. Reliability is a combination of time and
1116 // location variability. Confidence. is given by the
1117 // situation variability.
1119 // case 3: Broadcast mode. like avar(zt, zl, zc).
1120 // Reliability is given by the two-fold statement of at
1121 // least qT of the time in qL of the locations. Confidence.
1122 // is given by the situation variability.
1124 if (fabs(zt) > 3.1 || fabs(zl) > 3.1 || fabs(zc) > 3.1) {
1125 set_warn("Situations variables not optimal", 1);
1126 prop.kwx = mymax(prop.kwx, 1);
1130 // :38: Resolve standard deviations, page 19
1137 sgt = sgtd + tgtd / zt;
1138 // [Alg 5.11], situation variability
1139 vs = vs0 + (sgt * zt * sgt * zt) / (rt + zc * zc) + (sgl * zl * sgl * zl) / (rl + zc * zc);
1142 // :39: Resolve deviations, page 19
1145 propv.sgc = sqrt(sgt * sgt + sgl * sgl + vs);
1149 propv.sgc = sqrt(sgl * sgl + vs);
1152 yr = sqrt(sgt * sgt + sgl * sgl) * zt;
1153 propv.sgc = sqrt(vs);
1155 yr = sgt * zt + sgl * zl;
1156 propv.sgc = sqrt(vs);
1159 // [Alg 5.1], area variability
1160 avarv = prop.A_ref - vmd - yr - propv.sgc * zc;
1164 avarv = avarv * (29.0 - avarv) / (29.0 - 10.0 * avarv);
1170 // :45: Find to horizons, page 24
1172 * Here we use the terrain profile @pfl to find the two horizons. Output
1173 * consists of the horizon distances @dl and the horizon take-off angles
1174 * @the. If the path is line-of-sight, the routine sets both horizon
1175 * distances equal to @dist.
1178 void hzns(double pfl[], prop_type &prop)
1182 double xi, za, zb, qc, q, sb, sa;
1186 za = pfl[2] + prop.h_g[0];
1187 zb = pfl[np+2] + prop.h_g[1];
1188 qc = 0.5 * prop.gamma_e;
1190 prop.theta_ej[1] = (zb - za) / prop.d;
1191 prop.theta_ej[0] = prop.theta_ej[1] - q;
1192 prop.theta_ej[1] = -prop.theta_ej[1] - q;
1193 prop.d_Lj[0] = prop.d;
1194 prop.d_Lj[1] = prop.d;
1201 for (int i = 1; i < np; i++) {
1204 q = pfl[i+2] - (qc * sa + prop.theta_ej[0]) * sa - za;
1207 prop.theta_ej[0] += q / sa;
1213 q = pfl[i+2] - (qc * sb + prop.theta_ej[1]) * sb - zb;
1216 prop.theta_ej[1] += q / sb;
1225 // :53: Linear least square fit, page 28
1227 void zlsq1(double z[], const double &x1, const double &x2, double& z0, double& zn)
1230 double xn, xa, xb, x, a, b;
1234 xa = int(FORTRAN_DIM(x1 / z[1], 0.0));
1235 xb = xn - int(FORTRAN_DIM(xn, x2 / z[1]));
1238 xa = FORTRAN_DIM(xa, 1.0);
1239 xb = xn - FORTRAN_DIM(xn, xb + 1.0);
1248 a = 0.5 * (z[ja+2] + z[jb+2]);
1249 b = 0.5 * (z[ja+2] - z[jb+2]) * x;
1251 for (int i = 2; i <= n; ++i) {
1259 b = b * 12.0 / ((xa * xa + 2.0) * xa);
1261 zn = a + b * (xn - xb);
1265 // :52: Provide a quantile and reorders array @a, page 27
1267 double qtile(const int &nn, double a[], const int &ir)
1269 double q = 0.0, r; /* Initializations by KD2BD */
1270 int m, n, i, j, j1 = 0, i0 = 0, k; /* Initializations by KD2BD */
1276 k = mymin(mymax(0, ir), n);
1287 while (i <= n && a[i] >= q)
1294 while (j >= m && a[j] <= q)
1329 // :50: Standard normal complementary probability, page 26
1331 double qerf(const double &z)
1333 const double b1 = 0.319381530, b2 = -0.356563782, b3 = 1.781477937;
1334 const double b4 = -1.821255987, b5 = 1.330274429;
1335 const double rp = 4.317008, rrt2pi = 0.398942280;
1345 qerfv = exp(-0.5 * x * x) * rrt2pi * ((((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t;
1349 qerfv = 1.0 - qerfv;
1356 // :48: Find interdecile range of elevations, page 25
1358 * Using the terrain profile @pfl we find \Delta h, the interdecile range of
1359 * elevations between the two points @x1 and @x2.
1362 double dlthx(double pfl[], const double &x1, const double &x2)
1364 int np, ka, kb, n, k, j;
1365 double dlthxv, sn, xa, xb;
1373 if (xb - xa < 2.0) // exit out
1376 ka = (int)(0.1 * (xb - xa + 8.0));
1377 ka = mymin(mymax(4, ka), 25);
1381 assert((s = new double[n+2]) != 0);
1384 xb = (xb - xa) / sn;
1385 k = (int)(xa + 1.0);
1388 for (j = 0; j < n; j++) {
1390 while (xa > 0.0 && k < np) {
1395 s[j+2] = pfl[k+2] + (pfl[k+2] - pfl[k+1]) * xa;
1399 zlsq1(s, 0.0, sn, xa, xb);
1400 xb = (xb - xa) / sn;
1402 for (j = 0; j < n; j++) {
1407 dlthxv = qtile(n - 1, s + 2, ka - 1) - qtile(n - 1, s + 2, kb - 1);
1408 dlthxv /= 1.0 - 0.8 * exp(-(x2 - x1) / 50.0e3);
1415 // :41: Prepare model for point-to-point operation, page 22
1417 * This mode requires the terrain profile lying between the terminals. This
1418 * should be a sequence of surface elevations at points along the great
1419 * circle path joining the two points. It should start at the ground beneath
1420 * the first terminal and end at the ground beneath the second. In the
1421 * present routine it is assumed that the elevations are equispaced along
1422 * the path. They are stored in the array @pfl along with two defining
1426 * pfl[0] = np, the number of points in the path
1427 * pfl[1] = xi, the length of each increment
1430 void qlrpfl(double pfl[], int klimx, int mdvarx, prop_type &prop, propv_type &propv)
1433 double xl[2], q, za, zb;
1435 prop.d = pfl[0] * pfl[1];
1438 // :44: determine horizons and dh from pfl, page 23
1440 for (j = 0; j < 2; j++)
1441 xl[j] = mymin(15.0 * prop.h_g[j], 0.1 * prop.d_Lj[j]);
1443 xl[1] = prop.d - xl[1];
1444 prop.delta_h = dlthx(pfl, xl[0], xl[1]);
1446 if (prop.d_Lj[0] + prop.d_Lj[1] > 1.5 * prop.d) {
1447 // :45: Redo line-of-sight horizons, page 23
1450 * If the path is line-of-sight, we still need to know where
1451 * the horizons might have been, and so we turn to
1452 * techniques used in area prediction mode.
1454 zlsq1(pfl, xl[0], xl[1], za, zb);
1455 prop.h_e[0] = prop.h_g[0] + FORTRAN_DIM(pfl[2], za);
1456 prop.h_e[1] = prop.h_g[1] + FORTRAN_DIM(pfl[np+2], zb);
1458 for (j = 0; j < 2; j++)
1459 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)));
1461 q = prop.d_Lj[0] + prop.d_Lj[1];
1464 q = ((prop.d / q) * (prop.d / q));
1466 for (j = 0; j < 2; j++) {
1468 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)));
1472 for (j = 0; j < 2; j++) {
1473 q = sqrt(2.0 * prop.h_e[j] / prop.gamma_e);
1474 prop.theta_ej[j] = (0.65 * prop.delta_h * (q / prop.d_Lj[j] - 1.0) - 2.0 * prop.h_e[j]) / q;
1477 // :46: Transhorizon effective heights, page 23
1478 zlsq1(pfl, xl[0], 0.9 * prop.d_Lj[0], za, q);
1479 zlsq1(pfl, prop.d - 0.9 * prop.d_Lj[1], xl[1], q, zb);
1480 prop.h_e[0] = prop.h_g[0] + FORTRAN_DIM(pfl[2], za);
1481 prop.h_e[1] = prop.h_g[1] + FORTRAN_DIM(pfl[np+2], zb);
1485 propv.lvar = mymax(propv.lvar, 3);
1488 propv.mdvar = mdvarx;
1489 propv.lvar = mymax(propv.lvar, 4);
1501 //********************************************************
1502 //* Point-To-Point Mode Calculations *
1503 //********************************************************
1506 #ifdef WITH_POINT_TO_POINT
1509 void point_to_point(double elev[],
1510 double tht_m, // Transceiver above ground level
1511 double rht_m, // Receiver above groud level
1512 double eps_dielect, // Earth dielectric constant (rel. permittivity)
1513 double sgm_conductivity, // Earth conductivity (Siemens/m)
1514 double eno, // Atmospheric bending const, n-Units
1516 int radio_climate, // 1..7
1517 int pol, // 0 horizontal, 1 vertical
1518 double conf, // 0.01 .. .99, Fractions of situations
1519 double rel, // 0.01 .. .99, Fractions of time
1524 // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1525 // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1526 // 7-Maritime Temperate, Over Sea
1527 // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
1528 // errnum: 0- No Error.
1529 // 1- Warning: Some parameters are nearly out of range.
1530 // Results should be used with caution.
1531 // 2- Note: Default parameters have been substituted for impossible ones.
1532 // 3- Warning: A combination of parameters is out of range.
1533 // Results are probably invalid.
1534 // Other- Warning: Some parameters are out of range.
1535 // Results are probably invalid.
1547 prop.h_g[0] = tht_m; // Tx height above ground level
1548 prop.h_g[1] = rht_m; // Rx height above ground level
1549 propv.klim = radio_climate;
1550 prop.kwx = 0; // no error yet
1551 propv.lvar = 5; // initialize all constants
1552 prop.mdp = -1; // point-to-point
1556 dkm = (elev[1] * elev[0]) / 1000.0;
1557 xkm = elev[1] / 1000.0;
1559 ja = (long)(3.0 + 0.1 * elev[0]);
1561 for (i = ja - 1; i < jb; ++i)
1563 zsys /= (jb - ja + 1);
1566 qlrps(frq_mhz, zsys, q, pol, eps_dielect, sgm_conductivity, prop);
1567 qlrpfl(elev, propv.klim, propv.mdvar, prop, propv);
1568 fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
1569 q = prop.d - prop.d_L;
1572 strcpy(strmode, "Line-Of-Sight Mode");
1575 strcpy(strmode, "Single Horizon");
1579 strcpy(strmode, "Double Horizon");
1581 if (prop.d <= prop.d_Ls || prop.d <= prop.dx)
1582 strcat(strmode, ", Diffraction Dominant");
1585 if (prop.d > prop.dx)
1586 strcat(strmode, ", Troposcatter Dominant");
1589 dbloss = avar(zr, 0.0, zc, prop, propv) + fs;
1595 #ifdef WITH_POINT_TO_POINT_MDH
1597 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)
1599 // pol: 0-Horizontal, 1-Vertical
1600 // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1601 // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1602 // 7-Maritime Temperate, Over Sea
1603 // timepct, locpct, confpct: .01 to .99
1604 // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
1605 // propmode: Value Mode
1606 // -1 mode is undefined
1608 // 5 Single Horizon, Diffraction
1609 // 6 Single Horizon, Troposcatter
1610 // 9 Double Horizon, Diffraction
1611 // 10 Double Horizon, Troposcatter
1612 // errnum: 0- No Error.
1613 // 1- Warning: Some parameters are nearly out of range.
1614 // Results should be used with caution.
1615 // 2- Note: Default parameters have been substituted for impossible ones.
1616 // 3- Warning: A combination of parameters is out of range.
1617 // Results are probably invalid.
1618 // Other- Warning: Some parameters are out of range.
1619 // Results are probably invalid.
1625 double ztime, zloc, zconf;
1631 propmode = -1; // mode is undefined
1632 prop.h_g[0] = tht_m;
1633 prop.h_g[1] = rht_m;
1634 propv.klim = radio_climate;
1638 ztime = qerfi(timepct);
1639 zloc = qerfi(locpct);
1640 zconf = qerfi(confpct);
1643 dkm = (elev[1] * elev[0]) / 1000.0;
1644 xkm = elev[1] / 1000.0;
1646 ja = (long)(3.0 + 0.1 * elev[0]);
1648 for (i = ja - 1; i < jb; ++i)
1650 zsys /= (jb - ja + 1);
1653 qlrps(frq_mhz, zsys, q, pol, eps_dielect, sgm_conductivity, prop);
1654 qlrpfl(elev, propv.klim, propv.mdvar, prop, propa, propv);
1655 fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
1656 deltaH = prop.delta_h;
1657 q = prop.d - prop.d_L;
1660 propmode = 0; // Line-Of-Sight Mode
1663 propmode = 4; // Single Horizon
1666 propmode = 8; // Double Horizon
1668 if (prop.d <= prop.d_Ls || prop.d <= prop.dx)
1669 propmode += 1; // Diffraction Dominant
1672 if (prop.d > prop.dx)
1673 propmode += 2; // Troposcatter Dominant
1676 dbloss = avar(ztime, zloc, zconf, prop, propv) + fs; //avar(time,location,confidence)
1682 #ifdef WITH_POINT_TO_POINT_DH
1684 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)
1686 // pol: 0-Horizontal, 1-Vertical
1687 // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1688 // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1689 // 7-Maritime Temperate, Over Sea
1690 // conf, rel: .01 to .99
1691 // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
1692 // errnum: 0- No Error.
1693 // 1- Warning: Some parameters are nearly out of range.
1694 // Results should be used with caution.
1695 // 2- Note: Default parameters have been substituted for impossible ones.
1696 // 3- Warning: A combination of parameters is out of range.
1697 // Results are probably invalid.
1698 // Other- Warning: Some parameters are out of range.
1699 // Results are probably invalid.
1712 prop.h_g[0] = tht_m;
1713 prop.h_g[1] = rht_m;
1714 propv.klim = radio_climate;
1721 dkm = (elev[1] * elev[0]) / 1000.0;
1722 xkm = elev[1] / 1000.0;
1724 ja = (long)(3.0 + 0.1 * elev[0]);
1726 for (i = ja - 1; i < jb; ++i)
1728 zsys /= (jb - ja + 1);
1731 qlrps(frq_mhz, zsys, q, pol, eps_dielect, sgm_conductivity, prop);
1732 qlrpfl(elev, propv.klim, propv.mdvar, prop, propa, propv);
1733 fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
1734 deltaH = prop.delta_h;
1735 q = prop.d - prop.d_L;
1738 strcpy(strmode, "Line-Of-Sight Mode");
1741 strcpy(strmode, "Single Horizon");
1745 strcpy(strmode, "Double Horizon");
1747 if (prop.d <= prop.d_Ls || prop.d <= prop.dx)
1748 strcat(strmode, ", Diffraction Dominant");
1751 if (prop.d > prop.dx)
1752 strcat(strmode, ", Troposcatter Dominant");
1755 dbloss = avar(zr, 0.0, zc, prop, propv) + fs; //avar(time,location,confidence)
1762 //********************************************************
1763 //* Area Mode Calculations *
1764 //********************************************************
1769 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)
1771 // pol: 0-Horizontal, 1-Vertical
1772 // TSiteCriteria, RSiteCriteria:
1773 // 0 - random, 1 - careful, 2 - very careful
1774 // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1775 // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1776 // 7-Maritime Temperate, Over Sea
1777 // ModVar: 0 - Single: pctConf is "Time/Situation/Location", pctTime, pctLoc not used
1778 // 1 - Individual: pctTime is "Situation/Location", pctConf is "Confidence", pctLoc not used
1779 // 2 - Mobile: pctTime is "Time/Locations (Reliability)", pctConf is "Confidence", pctLoc not used
1780 // 3 - Broadcast: pctTime is "Time", pctLoc is "Location", pctConf is "Confidence"
1781 // pctTime, pctLoc, pctConf: .01 to .99
1782 // errnum: 0- No Error.
1783 // 1- Warning: Some parameters are nearly out of range.
1784 // Results should be used with caution.
1785 // 2- Note: Default parameters have been substituted for impossible ones.
1786 // 3- Warning: A combination of parameters is out of range.
1787 // Results are probably invalid.
1788 // Other- Warning: Some parameters are out of range.
1789 // Results are probably invalid.
1790 // NOTE: strmode is not used at this time.
1795 double zt, zl, zc, xlb;
1802 kst[0] = (int)TSiteCriteria;
1803 kst[1] = (int)RSiteCriteria;
1804 zt = qerfi(pctTime);
1806 zc = qerfi(pctConf);
1808 sgm = sgm_conductivity;
1809 prop.delta_h = deltaH;
1810 prop.h_g[0] = tht_m;
1811 prop.h_g[1] = rht_m;
1812 propv.klim = (long)radio_climate;
1815 ivar = (long)ModVar;
1817 qlrps(frq_mhz, 0.0, eno, ipol, eps, sgm, prop);
1818 qlra(kst, propv.klim, ivar, prop, propv);
1823 lrprop(dist_km * 1000.0, prop, propa);
1824 fs = 32.45 + 20.0 * log10(frq_mhz) + 20.0 * log10(prop.d / 1000.0);
1825 xlb = fs + avar(zt, zl, zc, prop, propv);