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 \*****************************************************************************/
15 // *************************************
16 // C++ routines for this program are taken from
17 // a translation of the FORTRAN code written by
18 // U.S. Department of Commerce NTIA/ITS
19 // Institute for Telecommunication Sciences
21 // Irregular Terrain Model (ITM) (Longley-Rice)
22 // *************************************
29 #define THIRD (1.0/3.0)
81 int mymin(const int &i, const int &j)
89 int mymax(const int &i, const int &j)
97 double mymin(const double &a, const double &b)
105 double mymax(const double &a, const double &b)
113 double FORTRAN_DIM(const double &x, const double &y)
115 // This performs the FORTRAN DIM function.
116 // result is x-y if x is greater than y; otherwise result is 0.0
124 double aknfe(const double &v2)
129 a=6.02+9.11*sqrt(v2)-1.27*v2;
131 a=12.953+4.343*log(v2);
136 double fht(const double& x, const double& pk)
143 /* if (pk < 1e-5 || x*pow(w,3.0) > 5495.0 ) */
145 if (pk < 1e-5 || (x*w*w*w) > 5495.0 )
150 fhtv=17.372*log(x)+fhtv;
154 fhtv=2.5e-5*x*x/pk-8.686*w-15.0;
159 fhtv=0.05751*x-4.343*log(x);
163 w=0.0134*x*exp(-0.005*x);
164 fhtv=(1.0-w)*fhtv+w*(17.372*log(x)-117.0);
171 double h0f(double r, double et)
173 double a[5]={25.0, 80.0, 177.0, 395.0, 705.0};
174 double b[5]={24.0, 45.0, 68.0, 80.0, 105.0};
196 /* x=pow(1.0/r,2.0); */
199 h0fv=4.343*log((a[it-1]*x+b[it-1])*x+1.0);
202 h0fv=(1.0-q)*h0fv+q*4.343*log((a[it]*x+b[it])*x+1.0);
207 double ahd(double td)
210 double a[3] = { 133.4, 104.6, 71.8};
211 double b[3] = {0.332e-3, 0.212e-3, 0.157e-3};
212 double c[3] = { -4.343, -1.086, 2.171};
223 return a[i]+b[i]*td+c[i]*log(td);
226 double adiff( double d, prop_type &prop, propa_type &propa)
228 complex<double> prop_zgnd(prop.zgndreal,prop.zgndimag);
229 static double wd1, xd1, afo, qk, aht, xht;
230 double a, q, pk, ds, th, wa, ar, wd, adiffv;
234 q=prop.hg[0]*prop.hg[1];
235 qk=prop.he[0]*prop.he[1]-q;
241 xd1=propa.dla+propa.tha/prop.gme;
242 q=(1.0-0.8*exp(-propa.dlsa/50e3))*prop.dh;
243 q*=0.78*exp(-pow(q/16.0,0.25));
244 afo=mymin(15.0,2.171*log(1.0+4.77e-4*prop.hg[0]*prop.hg[1]*prop.wn*q));
245 qk=1.0/abs(prop_zgnd);
249 for (int j=0; j<2; ++j)
251 /* a=0.5*pow(prop.dl[j],2.0)/prop.he[j]; */
252 a=0.5*(prop.dl[j]*prop.dl[j])/prop.he[j];
253 wa=pow(a*prop.wn,THIRD);
255 q=(1.607-pk)*151.0*wa*prop.dl[j]/a;
265 th=propa.tha+d*prop.gme;
267 /* q=0.0795775*prop.wn*ds*pow(th,2.0); */
268 q=0.0795775*prop.wn*ds*th*th;
269 adiffv=aknfe(q*prop.dl[0]/(ds+prop.dl[0]))+aknfe(q*prop.dl[1]/(ds+prop.dl[1]));
271 wa=pow(a*prop.wn,THIRD);
273 q=(1.607-pk)*151.0*wa*th+xht;
274 ar=0.05751*q-4.343*log(q)-aht;
275 q=(wd1+xd1/d)*mymin(((1.0-0.8*exp(-d/50e3))*prop.dh*prop.wn),6283.2);
276 wd=25.1/(25.1+sqrt(q));
277 adiffv=ar*wd+(1.0-wd)*adiffv+afo;
283 double ascat( double d, prop_type &prop, propa_type &propa)
285 static double ad, rr, etq, h0s;
286 double h0, r1, r2, z0, ss, et, ett, th, q;
291 ad=prop.dl[0]-prop.dl[1];
292 rr=prop.he[1]/prop.he[0];
300 etq=(5.67e-6*prop.ens-2.32e-3)*prop.ens+0.031;
311 th=prop.the[0]+prop.the[1]+d*prop.gme;
316 if (r1<0.2 && r2<0.2)
317 return 1001.0; // <==== early return
322 q=mymin(mymax(0.1,q),10.0);
323 z0=(d-ad)*(d+ad)*th*0.25/d;
324 /* et=(etq*exp(-pow(mymin(1.7,z0/8.0e3),6.0))+1.0)*z0/1.7556e3; */
325 temp=mymin(1.7,z0/8.0e3);
326 temp=temp*temp*temp*temp*temp*temp;
327 et=(etq*exp(-temp)+1.0)*z0/1.7556e3;
330 h0=(h0f(r1,ett)+h0f(r2,ett))*0.5;
331 h0+=mymin(h0,(1.38-log(ett))*log(ss)*log(q)*0.49);
332 h0=FORTRAN_DIM(h0,0.0);
335 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));
337 if (h0>15.0 && h0s>=0.0)
342 th=propa.tha+d*prop.gme;
343 /* ascatv=ahd(th*d)+4.343*log(47.7*prop.wn*pow(th,4.0))-0.1*(prop.ens-301.0)*exp(-th*d/40e3)+h0; */
345 ascatv=ahd(th*d)+4.343*log(47.7*prop.wn*th*th*th*th)-0.1*(prop.ens-301.0)*exp(-th*d/40e3)+h0;
352 double qerfi (double q)
355 double c0=2.515516698;
363 t=mymax(0.5-fabs(x),0.000001);
365 v=t-((c2*t+c1)*t+c0)/(((d3*t+d2)*t+d1)*t+1.0);
373 void qlrps( double fmhz, double zsys, double en0, int ipol, double eps, double sgm, prop_type &prop)
382 prop.ens*=exp(-zsys/9460.0);
384 prop.gme=gma*(1.0-0.04665*exp(prop.ens/179.3));
385 complex<double> zq, prop_zgnd(prop.zgndreal,prop.zgndimag);
386 zq=complex<double> (eps,376.62*sgm/prop.wn);
387 prop_zgnd=sqrt(zq-1.0);
390 prop_zgnd=prop_zgnd/zq;
392 prop.zgndreal=prop_zgnd.real();
393 prop.zgndimag=prop_zgnd.imag();
396 double abq_alos(complex<double> r)
398 return r.real()*r.real()+r.imag()*r.imag();
401 double alos(double d, prop_type &prop, propa_type &propa)
403 complex<double> prop_zgnd(prop.zgndreal,prop.zgndimag);
411 wls=0.021/(0.021+prop.wn*prop.dh/mymax(10e3,propa.dlsa));
417 q=(1.0-0.8*exp(-d/50e3))*prop.dh;
418 s=0.78*q*exp(-pow(q/16.0,0.25));
419 q=prop.he[0]+prop.he[1];
421 r=(sps-prop_zgnd)/(sps+prop_zgnd)*exp(-mymin(10.0,prop.wn*s*sps));
427 alosv=propa.emd*d+propa.aed;
428 q=prop.wn*prop.he[0]*prop.he[1]*2.0/d;
433 alosv=(-4.343*log(abq_alos(complex<double>(cos(q),-sin(q))+r))-alosv)*wls+alosv;
440 void qlra(int kst[], int klimx, int mdvarx, prop_type &prop, propv_type &propv)
444 for (int j=0; j<2; ++j)
447 prop.he[j]=prop.hg[j];
456 q*=sin(0.3141593*prop.hg[j]);
458 prop.he[j]=prop.hg[j]+(1.0+q)*exp(-mymin(20.0,2.0*prop.hg[j]/mymax(1e-3,prop.dh)));
461 q=sqrt(2.0*prop.he[j]/prop.gme);
462 prop.dl[j]=q*exp(-0.07*sqrt(prop.dh/mymax(prop.he[j],5.0)));
463 prop.the[j]=(0.65*prop.dh*(q/prop.dl[j]-1.0)-2.0*prop.he[j])/q;
467 propv.lvar=mymax(propv.lvar,3);
472 propv.lvar=mymax(propv.lvar,4);
482 void lrprop (double d, prop_type &prop, propa_type &propa) // PaulM_lrprop
484 static bool wlos, wscat;
485 static double dmin, xae;
486 complex<double> prop_zgnd(prop.zgndreal,prop.zgndimag);
487 double a0, a1, a2, a3, a4, a5, a6;
488 double d0, d1, d2, d3, d4, d5, d6;
496 propa.dls[j]=sqrt(2.0*prop.he[j]/prop.gme);
498 propa.dlsa=propa.dls[0]+propa.dls[1];
499 propa.dla=prop.dl[0]+prop.dl[1];
500 propa.tha=mymax(prop.the[0]+prop.the[1],-propa.dla*prop.gme);
504 if (prop.wn<0.838 || prop.wn>210.0)
505 prop.kwx=mymax(prop.kwx,1);
508 if (prop.hg[j]<1.0 || prop.hg[j]>1000.0)
509 prop.kwx=mymax(prop.kwx,1);
512 if (abs(prop.the[j]) >200e-3 || prop.dl[j]<0.1*propa.dls[j] || prop.dl[j]>3.0*propa.dls[j])
513 prop.kwx=mymax(prop.kwx,3);
515 if (prop.ens < 250.0 || prop.ens > 400.0 || prop.gme < 75e-9 || prop.gme > 250e-9 || prop_zgnd.real() <= abs(prop_zgnd.imag()) || prop.wn < 0.419 || prop.wn > 420.0)
521 if (prop.hg[j]<0.5 || prop.hg[j]>3000.0)
526 dmin=abs(prop.he[0]-prop.he[1])/200e-3;
527 q=adiff(0.0,prop,propa);
528 /* xae=pow(prop.wn*pow(prop.gme,2),-THIRD); */
529 xae=pow(prop.wn*prop.gme*prop.gme,-THIRD);
530 d3=mymax(propa.dlsa,1.3787*xae+propa.dla);
532 a3=adiff(d3,prop,propa);
533 a4=adiff(d4,prop,propa);
534 propa.emd=(a4-a3)/(d4-d3);
535 propa.aed=a3-propa.emd*d3;
546 if (prop.dist>1000e3)
547 prop.kwx=mymax(prop.kwx,1);
550 prop.kwx=mymax(prop.kwx,3);
552 if (prop.dist<1e3 || prop.dist>2000e3)
556 if (prop.dist<propa.dlsa)
560 q=alos(0.0,prop,propa);
562 a2=propa.aed+d2*propa.emd;
563 d0=1.908*prop.wn*prop.he[0]*prop.he[1];
567 d0=mymin(d0,0.5*propa.dla);
568 d1=d0+0.25*(propa.dla-d0);
572 d1=mymax(-propa.aed/propa.emd,0.25*propa.dla);
574 a1=alos(d1,prop,propa);
579 a0=alos(d0,prop,propa);
581 propa.ak2=mymax(0.0,((d2-d0)*(a1-a0)-(d1-d0)*(a2-a0))/((d2-d0)*log(d1/d0)-(d1-d0)*q));
582 wq=propa.aed>=0.0 || propa.ak2>0.0;
586 propa.ak1=(a2-a0-propa.ak2*q)/(d2-d0);
591 propa.ak2=FORTRAN_DIM(a2,a0)/q;
601 propa.ak1=FORTRAN_DIM(a2,a1)/(d2-d1);
608 propa.ael=a2-propa.ak1*d2-propa.ak2*log(d2);
613 prop.aref=propa.ael+propa.ak1*prop.dist+propa.ak2*log(prop.dist);
616 if (prop.dist<=0.0 || prop.dist>=propa.dlsa)
620 q=ascat(0.0,prop,propa);
623 a6=ascat(d6,prop,propa);
624 a5=ascat(d5,prop,propa);
628 propa.ems=(a6-a5)/200e3;
629 propa.dx=mymax(propa.dlsa,mymax(propa.dla+0.3*xae*log(47.7*prop.wn),(a5-propa.aed-propa.ems*d5)/(propa.emd-propa.ems)));
631 propa.aes=(propa.emd-propa.ems)*propa.dx+propa.aed;
644 if (prop.dist>propa.dx)
645 prop.aref=propa.aes+propa.ems*prop.dist;
647 prop.aref=propa.aed+propa.emd*prop.dist;
650 prop.aref=mymax(prop.aref,0.0);
653 double curve (double const &c1, double const &c2, double const &x1, double const &x2, double const &x3, double const &de)
655 /* return (c1+c2/(1.0+pow((de-x2)/x3,2.0)))*pow(de/x1,2.0)/(1.0+pow(de/x1,2.0)); */
665 return (c1+c2/(1.0+temp1))*temp2/(1.0+temp2);
668 double avar(double zzt, double zzl, double zzc,prop_type &prop, propv_type &propv)
671 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;
673 double bv1[7]={-9.67,-0.62,1.26,-9.21,-0.62,-0.39,3.15};
674 double bv2[7]={12.7,9.19,15.5,9.05,9.19,2.86,857.9};
675 double xv1[7]={144.9e3,228.9e3,262.6e3,84.1e3,228.9e3,141.7e3,2222.e3};
676 double xv2[7]={190.3e3,205.2e3,185.2e3,101.1e3,205.2e3,315.9e3,164.8e3};
677 double xv3[7]={133.8e3,143.6e3,99.8e3,98.6e3,143.6e3,167.4e3,116.3e3};
678 double bsm1[7]={2.13,2.66,6.11,1.98,2.68,6.86,8.51};
679 double bsm2[7]={159.5,7.67,6.65,13.11,7.16,10.38,169.8};
680 double xsm1[7]={762.2e3,100.4e3,138.2e3,139.1e3,93.7e3,187.8e3,609.8e3};
681 double xsm2[7]={123.6e3,172.5e3,242.2e3,132.7e3,186.8e3,169.6e3,119.9e3};
682 double xsm3[7]={94.5e3,136.4e3,178.6e3,193.5e3,133.5e3,108.9e3,106.6e3};
683 double bsp1[7]={2.11,6.87,10.08,3.68,4.75,8.58,8.43};
684 double bsp2[7]={102.3,15.53,9.60,159.3,8.12,13.97,8.19};
685 double xsp1[7]={636.9e3,138.7e3,165.3e3,464.4e3,93.2e3,216.0e3,136.2e3};
686 double xsp2[7]={134.8e3,143.7e3,225.7e3,93.1e3,135.9e3,152.0e3,188.5e3};
687 double xsp3[7]={95.6e3,98.6e3,129.7e3,94.2e3,113.4e3,122.7e3,122.9e3};
688 double bsd1[7]={1.224,0.801,1.380,1.000,1.224,1.518,1.518};
689 double bzd1[7]={1.282,2.161,1.282,20.,1.282,1.282,1.282};
690 double bfm1[7]={1.0,1.0,1.0,1.0,0.92,1.0,1.0};
691 double bfm2[7]={0.0,0.0,0.0,0.0,0.25,0.0,0.0};
692 double bfm3[7]={0.0,0.0,0.0,0.0,1.77,0.0,0.0};
693 double bfp1[7]={1.0,0.93,1.0,0.93,0.93,1.0,1.0};
694 double bfp2[7]={0.0,0.31,0.0,0.19,0.31,0.0,0.0};
695 double bfp3[7]={0.0,2.00,0.0,1.79,2.00,0.0,0.0};
697 double rt=7.8, rl=24.0, avarv, q, vs, zt, zl, zc;
699 int temp_klim=propv.klim-1;
706 if (propv.klim<=0 || propv.klim>7)
710 prop.kwx=mymax(prop.kwx,2);
713 cv1 = bv1[temp_klim];
714 cv2 = bv2[temp_klim];
715 yv1 = xv1[temp_klim];
716 yv2 = xv2[temp_klim];
717 yv3 = xv3[temp_klim];
718 csm1=bsm1[temp_klim];
719 csm2=bsm2[temp_klim];
720 ysm1=xsm1[temp_klim];
721 ysm2=xsm2[temp_klim];
722 ysm3=xsm3[temp_klim];
723 csp1=bsp1[temp_klim];
724 csp2=bsp2[temp_klim];
725 ysp1=xsp1[temp_klim];
726 ysp2=xsp2[temp_klim];
727 ysp3=xsp3[temp_klim];
728 csd1=bsd1[temp_klim];
730 cfm1=bfm1[temp_klim];
731 cfm2=bfm2[temp_klim];
732 cfm3=bfm3[temp_klim];
733 cfp1=bfp1[temp_klim];
734 cfp2=bfp2[temp_klim];
735 cfp3=bfp3[temp_klim];
752 prop.kwx=mymax(prop.kwx,2);
756 q=log(0.133*prop.wn);
758 /* gm=cfm1+cfm2/(pow(cfm3*q,2.0)+1.0); */
759 /* gp=cfp1+cfp2/(pow(cfp3*q,2.0)+1.0); */
761 gm=cfm1+cfm2/((cfm3*q*cfm3*q)+1.0);
762 gp=cfp1+cfp2/((cfp3*q*cfp3*q)+1.0);
765 dexa=sqrt(18e6*prop.he[0])+sqrt(18e6*prop.he[1])+pow((575.7e12/prop.wn),THIRD);
769 de=130e3*prop.dist/dexa;
771 de=130e3+prop.dist-dexa;
774 vmd=curve(cv1,cv2,yv1,yv2,yv3,de);
775 sgtm=curve(csm1,csm2,ysm1,ysm2,ysm3,de)*gm;
776 sgtp=curve(csp1,csp2,ysp1,ysp2,ysp3,de)*gp;
784 q=(1.0-0.8*exp(-prop.dist/50e3))*prop.dh*prop.wn;
792 /* vs0=pow(5.0+3.0*exp(-de/100e3),2.0); */
793 vs0=(5.0+3.0*exp(-de/100e3));
819 if (fabs(zt)>3.1 || fabs(zl)>3.1 || fabs(zc)>3.1)
820 prop.kwx=mymax(prop.kwx,1);
831 /* vs=vs0+pow(sgt*zt,2.0)/(rt+zc*zc)+pow(sgl*zl,2.0)/(rl+zc*zc); */
832 vs=vs0+(sgt*zt*sgt*zt)/(rt+zc*zc)+(sgl*zl*sgl*zl)/(rl+zc*zc);
837 propv.sgc=sqrt(sgt*sgt+sgl*sgl+vs);
843 propv.sgc=sqrt(sgl*sgl+vs);
848 yr=sqrt(sgt*sgt+sgl*sgl)*zt;
858 avarv=prop.aref-vmd-yr-propv.sgc*zc;
861 avarv=avarv*(29.0-avarv)/(29.0-10.0*avarv);
866 void hzns(double pfl[], prop_type &prop)
870 double xi, za, zb, qc, q, sb, sa;
874 za=pfl[2]+prop.hg[0];
875 zb=pfl[np+2]+prop.hg[1];
878 prop.the[1]=(zb-za)/prop.dist;
879 prop.the[0]=prop.the[1]-q;
880 prop.the[1]=-prop.the[1]-q;
881 prop.dl[0]=prop.dist;
882 prop.dl[1]=prop.dist;
890 for (int i=1; i<np; i++)
894 q=pfl[i+2]-(qc*sa+prop.the[0])*sa-za;
905 q=pfl[i+2]-(qc*sb+prop.the[1])*sb-zb;
917 void z1sq1 (double z[], const double &x1, const double &x2, double& z0, double& zn)
920 double xn, xa, xb, x, a, b;
924 xa=int(FORTRAN_DIM(x1/z[1],0.0));
925 xb=xn-int(FORTRAN_DIM(xn,x2/z[1]));
929 xa=FORTRAN_DIM(xa,1.0);
930 xb=xn-FORTRAN_DIM(xn,xb+1.0);
939 a=0.5*(z[ja+2]+z[jb+2]);
940 b=0.5*(z[ja+2]-z[jb+2])*x;
942 for (int i=2; i<=n; ++i)
951 b=b*12.0/((xa*xa+2.0)*xa);
956 double qtile (const int &nn, double a[], const int &ir)
958 double q=0.0, r; /* Initializations by KD2BD */
959 int m, n, i, j, j1=0, i0=0, k; /* Initializations by KD2BD */
965 k=mymin(mymax(0,ir),n);
978 while (i<=n && a[i]>=q)
985 while (j>=m && a[j]<=q)
1024 double qerf(const double &z)
1026 double b1=0.319381530, b2=-0.356563782, b3=1.781477937;
1027 double b4=-1.821255987, b5=1.330274429;
1028 double rp=4.317008, rrt2pi=0.398942280;
1039 qerfv=exp(-0.5*x*x)*rrt2pi*((((b5*t+b4)*t+b3)*t+b2)*t+b1)*t;
1048 double d1thx(double pfl[], const double &x1, const double &x2)
1050 int np, ka, kb, n, k, j;
1051 double d1thxv, sn, xa, xb;
1059 if (xb-xa<2.0) // exit out
1062 ka=(int)(0.1*(xb-xa+8.0));
1063 ka=mymin(mymax(4,ka),25);
1067 assert((s=new double[n+2])!=0);
1076 while (xa>0.0 && k<np)
1082 s[j+2]=pfl[k+2]+(pfl[k+2]-pfl[k+1])*xa;
1086 z1sq1(s,0.0,sn,xa,xb);
1095 d1thxv=qtile(n-1,s+2,ka-1)-qtile(n-1,s+2,kb-1);
1096 d1thxv/=1.0-0.8*exp(-(x2-x1)/50.0e3);
1102 void qlrpfl (double pfl[], int klimx, int mdvarx, prop_type &prop, propa_type &propa, propv_type &propv)
1105 double xl[2], q, za, zb;
1107 prop.dist=pfl[0]*pfl[1];
1112 xl[j]=mymin(15.0*prop.hg[j],0.1*prop.dl[j]);
1114 xl[1]=prop.dist-xl[1];
1115 prop.dh=d1thx(pfl,xl[0],xl[1]);
1117 if (prop.dl[0]+prop.dl[1]>1.5*prop.dist)
1119 z1sq1(pfl,xl[0],xl[1],za,zb);
1120 prop.he[0]=prop.hg[0]+FORTRAN_DIM(pfl[2],za);
1121 prop.he[1]=prop.hg[1]+FORTRAN_DIM(pfl[np+2],zb);
1124 prop.dl[j]=sqrt(2.0*prop.he[j]/prop.gme)*exp(-0.07*sqrt(prop.dh/mymax(prop.he[j],5.0)));
1126 q=prop.dl[0]+prop.dl[1];
1130 /* q=pow(prop.dist/q,2.0); */
1131 q=((prop.dist/q)*(prop.dist/q));
1136 prop.dl[j]=sqrt(2.0*prop.he[j]/prop.gme)*exp(-0.07*sqrt(prop.dh/mymax(prop.he[j],5.0)));
1142 q=sqrt(2.0*prop.he[j]/prop.gme);
1143 prop.the[j]=(0.65*prop.dh*(q/prop.dl[j]-1.0)-2.0*prop.he[j])/q;
1149 z1sq1(pfl,xl[0],0.9*prop.dl[0],za,q);
1150 z1sq1(pfl,prop.dist-0.9*prop.dl[1],xl[1],q,zb);
1151 prop.he[0]=prop.hg[0]+FORTRAN_DIM(pfl[2],za);
1152 prop.he[1]=prop.hg[1]+FORTRAN_DIM(pfl[np+2],zb);
1156 propv.lvar=mymax(propv.lvar,3);
1161 propv.lvar=mymax(propv.lvar,4);
1170 lrprop(0.0,prop,propa);
1173 double deg2rad(double d)
1175 return d*3.1415926535897/180.0;
1178 //********************************************************
1179 //* Point-To-Point Mode Calculations *
1180 //********************************************************
1182 void point_to_point(double elev[], double tht_m, double rht_m, double eps_dielect, double sgm_conductivity, double eno_ns_surfref, double frq_mhz, int radio_climate, int pol, double conf, double rel, double &dbloss, char *strmode, int &errnum)
1184 // pol: 0-Horizontal, 1-Vertical
1185 // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1186 // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1187 // 7-Maritime Temperate, Over Sea
1188 // conf, rel: .01 to .99
1189 // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
1190 // errnum: 0- No Error.
1191 // 1- Warning: Some parameters are nearly out of range.
1192 // Results should be used with caution.
1193 // 2- Note: Default parameters have been substituted for impossible ones.
1194 // 3- Warning: A combination of parameters is out of range.
1195 // Results are probably invalid.
1196 // Other- Warning: Some parameters are out of range.
1197 // Results are probably invalid.
1205 double eno, enso, q;
1212 propv.klim=radio_climate;
1219 dkm=(elev[1]*elev[0])/1000.0;
1227 /* ja = 3.0 + 0.1 * elev[0]; */
1228 ja=(long)(3.0+0.1* elev[0]);
1232 for (i=ja-1; i<jb; ++i)
1240 qlrps(frq_mhz,zsys,q,pol,eps_dielect,sgm_conductivity,prop);
1241 qlrpfl(elev,propv.klim,propv.mdvar,prop,propa,propv);
1242 fs=32.45+20.0*log10(frq_mhz)+20.0*log10(prop.dist/1000.0);
1243 q=prop.dist-propa.dla;
1246 strcpy(strmode,"Line-Of-Sight Mode");
1251 strcpy(strmode,"Single Horizon");
1253 else if (int(q)>0.0)
1254 strcpy(strmode,"Double Horizon");
1256 if (prop.dist<=propa.dlsa || prop.dist <= propa.dx)
1257 strcat(strmode,", Diffraction Dominant");
1259 else if (prop.dist>propa.dx)
1260 strcat(strmode, ", Troposcatter Dominant");
1263 dbloss=avar(zr,0.0,zc,prop,propv)+fs;
1267 void point_to_pointMDH (double elev[], double tht_m, double rht_m, double eps_dielect, double sgm_conductivity, double eno_ns_surfref, double frq_mhz, int radio_climate, int pol, double timepct, double locpct, double confpct, double &dbloss, int &propmode, double &deltaH, int &errnum)
1269 // pol: 0-Horizontal, 1-Vertical
1270 // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1271 // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1272 // 7-Maritime Temperate, Over Sea
1273 // timepct, locpct, confpct: .01 to .99
1274 // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
1275 // propmode: Value Mode
1276 // -1 mode is undefined
1278 // 5 Single Horizon, Diffraction
1279 // 6 Single Horizon, Troposcatter
1280 // 9 Double Horizon, Diffraction
1281 // 10 Double Horizon, Troposcatter
1282 // errnum: 0- No Error.
1283 // 1- Warning: Some parameters are nearly out of range.
1284 // Results should be used with caution.
1285 // 2- Note: Default parameters have been substituted for impossible ones.
1286 // 3- Warning: A combination of parameters is out of range.
1287 // Results are probably invalid.
1288 // Other- Warning: Some parameters are out of range.
1289 // Results are probably invalid.
1295 double ztime, zloc, zconf;
1296 double eno, enso, q;
1301 propmode=-1; // mode is undefined
1304 propv.klim=radio_climate;
1308 ztime=qerfi(timepct);
1310 zconf=qerfi(confpct);
1313 dkm=(elev[1]*elev[0])/1000.0;
1321 /* ja = 3.0 + 0.1 * elev[0]; */
1322 ja=(long)(3.0+0.1*elev[0]);
1325 for (i=ja-1; i<jb; ++i)
1333 qlrps(frq_mhz,zsys,q,pol,eps_dielect,sgm_conductivity,prop);
1334 qlrpfl(elev,propv.klim,propv.mdvar,prop,propa,propv);
1335 fs=32.45+20.0*log10(frq_mhz)+20.0*log10(prop.dist/1000.0);
1337 q=prop.dist-propa.dla;
1340 propmode=0; // Line-Of-Sight Mode
1344 propmode=4; // Single Horizon
1346 else if (int(q)>0.0)
1347 propmode=8; // Double Horizon
1349 if (prop.dist<=propa.dlsa || prop.dist<=propa.dx)
1350 propmode+=1; // Diffraction Dominant
1352 else if (prop.dist>propa.dx)
1353 propmode+=2; // Troposcatter Dominant
1356 dbloss=avar(ztime, zloc, zconf, prop, propv)+fs; //avar(time,location,confidence)
1360 void point_to_pointDH (double elev[], double tht_m, double rht_m, double eps_dielect, double sgm_conductivity, double eno_ns_surfref, double frq_mhz, int radio_climate, int pol, double conf, double rel, double &dbloss, double &deltaH, int &errnum)
1362 // pol: 0-Horizontal, 1-Vertical
1363 // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1364 // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1365 // 7-Maritime Temperate, Over Sea
1366 // conf, rel: .01 to .99
1367 // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
1368 // errnum: 0- No Error.
1369 // 1- Warning: Some parameters are nearly out of range.
1370 // Results should be used with caution.
1371 // 2- Note: Default parameters have been substituted for impossible ones.
1372 // 3- Warning: A combination of parameters is out of range.
1373 // Results are probably invalid.
1374 // Other- Warning: Some parameters are out of range.
1375 // Results are probably invalid.
1383 double eno, enso, q;
1390 propv.klim=radio_climate;
1397 dkm=(elev[1]*elev[0])/1000.0;
1405 /* ja = 3.0 + 0.1 * elev[0]; */
1406 ja=(long)(3.0+0.1*elev[0]);
1410 for (i=ja-1; i<jb; ++i)
1418 qlrps(frq_mhz,zsys,q,pol,eps_dielect,sgm_conductivity,prop);
1419 qlrpfl(elev,propv.klim,propv.mdvar,prop,propa,propv);
1420 fs=32.45+20.0*log10(frq_mhz)+20.0*log10(prop.dist/1000.0);
1422 q=prop.dist-propa.dla;
1425 strcpy(strmode,"Line-Of-Sight Mode");
1429 strcpy(strmode,"Single Horizon");
1431 else if (int(q)>0.0)
1432 strcpy(strmode,"Double Horizon");
1434 if (prop.dist<=propa.dlsa || prop.dist <= propa.dx)
1435 strcat(strmode,", Diffraction Dominant");
1437 else if (prop.dist>propa.dx)
1438 strcat(strmode, ", Troposcatter Dominant");
1441 dbloss=avar(zr,0.0,zc,prop,propv)+fs; //avar(time,location,confidence)
1446 //********************************************************
1447 //* Area Mode Calculations *
1448 //********************************************************
1450 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_ns_surfref, double frq_mhz, int radio_climate, int pol, double pctTime, double pctLoc, double pctConf, double &dbloss, char *strmode, int &errnum)
1452 // pol: 0-Horizontal, 1-Vertical
1453 // TSiteCriteria, RSiteCriteria:
1454 // 0 - random, 1 - careful, 2 - very careful
1455 // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1456 // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1457 // 7-Maritime Temperate, Over Sea
1458 // ModVar: 0 - Single: pctConf is "Time/Situation/Location", pctTime, pctLoc not used
1459 // 1 - Individual: pctTime is "Situation/Location", pctConf is "Confidence", pctLoc not used
1460 // 2 - Mobile: pctTime is "Time/Locations (Reliability)", pctConf is "Confidence", pctLoc not used
1461 // 3 - Broadcast: pctTime is "Time", pctLoc is "Location", pctConf is "Confidence"
1462 // pctTime, pctLoc, pctConf: .01 to .99
1463 // errnum: 0- No Error.
1464 // 1- Warning: Some parameters are nearly out of range.
1465 // Results should be used with caution.
1466 // 2- Note: Default parameters have been substituted for impossible ones.
1467 // 3- Warning: A combination of parameters is out of range.
1468 // Results are probably invalid.
1469 // Other- Warning: Some parameters are out of range.
1470 // Results are probably invalid.
1471 // NOTE: strmode is not used at this time.
1476 double zt, zl, zc, xlb;
1479 double eps, eno, sgm;
1483 kst[0]=(int)TSiteCriteria;
1484 kst[1]=(int)RSiteCriteria;
1489 sgm=sgm_conductivity;
1494 /* propv.klim = (__int32) radio_climate; */
1495 propv.klim=(long)radio_climate;
1500 qlrps(frq_mhz, 0.0, eno, ipol, eps, sgm, prop);
1501 qlra(kst, propv.klim, ivar, prop, propv);
1506 lrprop(dist_km*1000.0, prop, propa);
1507 fs=32.45+20.0*log10(frq_mhz)+20.0*log10(prop.dist/1000.0);
1508 xlb=fs+avar(zt, zl, zc, prop, propv);
1517 double ITMAreadBLoss(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_ns_surfref, double frq_mhz, int radio_climate, int pol, double pctTime, double pctLoc, double pctConf)
1523 area(ModVar,deltaH,tht_m,rht_m,dist_km,TSiteCriteria,RSiteCriteria, eps_dielect,sgm_conductivity,eno_ns_surfref, frq_mhz,radio_climate,pol,pctTime,pctLoc, pctConf,dbloss,strmode,errnum);