1 /****************************************************************************
3 * This file was obtained from ftp://flattop.its.bldrdoc.gov/itm/ITMDLL.cpp *
4 * on Jan. 10, 2004 and is public domain. It was modified by J. D. McDonald *
5 * to remove Microsoft Windows dll-isms and to correct one case where a *
6 * compiler found an ambguity in overloaded calls. It was further modified *
7 * by John A. Magliacane to remove unused variables, unneeded #includes, *
8 * and to replace pow() statements with explicit multiplications wherever *
9 * possible to increase execution speed and improve accuracy. *
11 ****************************************************************************/
13 // *************************************
14 // C++ routines for this program are taken from
15 // a translation of the FORTRAN code written by
16 // U.S. Department of Commerce NTIA/ITS
17 // Institute for Telecommunication Sciences
19 // Irregular Terrain Model (ITM) (Longley-Rice)
20 // *************************************
27 #define THIRD (1.0/3.0)
75 int mymin(const int &i, const int &j)
83 int mymax(const int &i, const int &j)
91 double mymin(const double &a, const double &b)
99 double mymax(const double &a, const double &b)
107 double FORTRAN_DIM(const double &x, const double &y)
109 /* This performs the FORTRAN DIM function. Result is x-y
110 if x is greater than y; otherwise result is 0.0 */
118 double aknfe(const double &v2)
123 a=6.02+9.11*sqrt(v2)-1.27*v2;
125 a=12.953+4.343*log(v2);
129 double fht(const double& x, const double& pk)
137 /* if (pk<1e-5 || x*pow(w,3.0) > 5495.0) */
138 if (pk<1.0e-5 || x*w*w*w > 5495.0)
143 fhtv=17.372*log(x)+fhtv;
146 fhtv=2.5e-5*x*x/pk-8.686*w-15.0;
151 fhtv=0.05751*x-4.343*log(x);
155 w=0.0134*x*exp(-0.005*x);
156 fhtv=(1.0-w)*fhtv+w*(17.372*log(x)-117.0);
162 double h0f(double r, double et)
164 double a[5]={25.0, 80.0, 177.0, 395.0, 705.0};
165 double b[5]={24.0, 45.0, 68.0, 80.0, 105.0};
187 /* x=pow(1.0/r,2.0); */
192 h0fv=4.343*log((a[it-1]*x+b[it-1])*x+1.0);
195 h0fv=(1.0-q)*h0fv+q*4.343*log((a[it]*x+b[it])*x+1.0);
200 double ahd(double td)
203 double a[3]={ 133.4, 104.6, 71.8};
204 double b[3]={0.332e-3, 0.212e-3, 0.157e-3};
205 double c[3]={ -4.343, -1.086, 2.171};
216 return a[i]+b[i]*td+c[i]*log(td);
219 double adiff(double d, prop_type &prop, propa_type &propa)
221 complex<double> prop_zgnd(prop.zgndreal,prop.zgndimag);
222 static double wd1, xd1, afo, qk, aht, xht;
223 double a, q, pk, ds, th, wa, ar, wd, adiffv;
227 q=prop.hg[0]*prop.hg[1];
228 qk=prop.he[0]*prop.he[1]-q;
234 xd1=propa.dla+propa.tha/prop.gme;
235 q=(1.0-0.8*exp(-propa.dlsa/50e3))*prop.dh;
236 q*=0.78*exp(-pow(q/16.0,0.25));
237 afo=mymin(15.0,2.171*log(1.0+4.77e-4*prop.hg[0]*prop.hg[1]*prop.wn*q));
238 qk=1.0/abs(prop_zgnd);
242 for (int j=0; j<2; ++j)
244 /* a=0.5*pow(prop.dl[j],2.0)/prop.he[j]; */
245 a=0.5*(prop.dl[j]*prop.dl[j])/prop.he[j];
246 wa=pow(a*prop.wn,THIRD);
248 q=(1.607-pk)*151.0*wa*prop.dl[j]/a;
258 th=propa.tha+d*prop.gme;
260 /* q=0.0795775*prop.wn*ds*pow(th,2.0); */
261 q=0.0795775*prop.wn*ds*th*th;
262 adiffv=aknfe(q*prop.dl[0]/(ds+prop.dl[0]))+aknfe(q*prop.dl[1]/(ds+prop.dl[1]));
264 wa=pow(a*prop.wn,THIRD);
266 q=(1.607-pk)*151.0*wa*th+xht;
267 ar=0.05751*q-4.343*log(q)-aht;
268 q=(wd1+xd1/d)*mymin(((1.0-0.8*exp(-d/50e3))*prop.dh*prop.wn),6283.2);
269 wd=25.1/(25.1+sqrt(q));
270 adiffv=ar*wd+(1.0-wd)*adiffv+afo;
276 double ascat( double d, prop_type &prop, propa_type &propa)
278 static double ad, rr, etq, h0s;
279 double h0, r1, r2, z0, ss, et, ett, th, q;
284 ad=prop.dl[0]-prop.dl[1];
285 rr=prop.he[1]/prop.he[0];
293 etq=(5.67e-6*prop.ens-2.32e-3)*prop.ens+0.031;
304 th=prop.the[0]+prop.the[1]+d*prop.gme;
309 if (r1<0.2 && r2<0.2)
310 return 1001.0; // <==== early return
315 q=mymin(mymax(0.1,q),10.0);
316 z0=(d-ad)*(d+ad)*th*0.25/d;
317 /* et=(etq*exp(-pow(mymin(1.7,z0/8.0e3),6.0))+1.0)*z0/1.7556e3; */
319 temp=mymin(1.7,z0/8.0e3);
320 temp=temp*temp*temp*temp*temp*temp;
321 et=(etq*exp(-temp)+1.0)*z0/1.7556e3;
324 h0=(h0f(r1,ett)+h0f(r2,ett))*0.5;
325 h0+=mymin(h0,(1.38-log(ett))*log(ss)*log(q)*0.49);
326 h0=FORTRAN_DIM(h0,0.0);
330 /* 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)); */
332 temp=((1.0+1.4142/r1)*(1.0+1.4142/r2));
333 h0=et*h0+(1.0-et)*4.343*log((temp*temp)*(r1+r2)/(r1+r2+2.8284));
336 if (h0>15.0 && h0s>=0.0)
341 th=propa.tha+d*prop.gme;
342 /* 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; */
343 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;
349 double qerfi(double q)
352 double c0=2.515516698;
360 t=mymax(0.5-fabs(x),0.000001);
362 v=t-((c2*t+c1)*t+c0)/(((d3*t+d2)*t+d1)*t+1.0);
370 void qlrps(double fmhz, double zsys, double en0, int ipol, double eps, double sgm, prop_type &prop)
378 prop.ens*=exp(-zsys/9460.0);
380 prop.gme=gma*(1.0-0.04665*exp(prop.ens/179.3));
381 complex<double> zq, prop_zgnd(prop.zgndreal,prop.zgndimag);
382 zq=complex<double> (eps,376.62*sgm/prop.wn);
383 prop_zgnd=sqrt(zq-1.0);
386 prop_zgnd=prop_zgnd/zq;
388 prop.zgndreal=prop_zgnd.real();
389 prop.zgndimag=prop_zgnd.imag();
392 double abq_alos (complex<double> r)
394 return r.real()*r.real()+r.imag()*r.imag();
397 double alos(double d, prop_type &prop, propa_type &propa)
399 complex<double> prop_zgnd(prop.zgndreal,prop.zgndimag);
407 wls=0.021/(0.021+prop.wn*prop.dh/mymax(10e3,propa.dlsa));
413 q=(1.0-0.8*exp(-d/50e3))*prop.dh;
414 s=0.78*q*exp(-pow(q/16.0,0.25));
415 q=prop.he[0]+prop.he[1];
417 r=(sps-prop_zgnd)/(sps+prop_zgnd)*exp(-mymin(10.0,prop.wn*s*sps));
423 alosv=propa.emd*d+propa.aed;
424 q=prop.wn*prop.he[0]*prop.he[1]*2.0/d;
429 alosv=(-4.343*log(abq_alos(complex<double>(cos(q),-sin(q))+r))-alosv)*wls+alosv;
436 void qlra(int kst[], int klimx, int mdvarx, prop_type &prop, propv_type &propv)
440 for (int j=0; j<2; ++j)
443 prop.he[j]=prop.hg[j];
452 q*=sin(0.3141593*prop.hg[j]);
454 prop.he[j]=prop.hg[j]+(1.0+q)*exp(-mymin(20.0,2.0*prop.hg[j]/mymax(1e-3,prop.dh)));
457 q=sqrt(2.0*prop.he[j]/prop.gme);
458 prop.dl[j]=q*exp(-0.07*sqrt(prop.dh/mymax(prop.he[j],5.0)));
459 prop.the[j]=(0.65*prop.dh*(q/prop.dl[j]-1.0)-2.0*prop.he[j])/q;
463 propv.lvar=mymax(propv.lvar,3);
468 propv.lvar=mymax(propv.lvar,4);
478 void freds_lrprop (double d, prop_type &prop, propa_type &propa)
482 static bool wlos, wscat;
483 static double dmin, xae;
484 complex<double> prop_zgnd(prop.zgndreal,prop.zgndimag);
485 double a0, a1, a2, a3, a4, a5, a6;
486 double d0, d1, d2, d3, d4, d5, d6;
493 propa.dls[j]=sqrt(2.0*prop.he[j]/prop.gme);
495 propa.dlsa=propa.dls[0]+propa.dls[1];
496 propa.dla=prop.dl[0]+prop.dl[1];
497 propa.tha=mymax(prop.the[0]+prop.the[1],-propa.dla*prop.gme);
501 if (prop.wn<0.838 || prop.wn>210.0)
502 prop.kwx=mymax(prop.kwx,1);
505 if (prop.hg[j]<1.0 || prop.hg[j]>1000.0)
506 prop.kwx=mymax(prop.kwx,1);
509 if (abs(prop.the[j]) >200e-3 || prop.dl[j]<0.1*propa.dls[j] || prop.dl[j]>3.0*propa.dls[j])
510 prop.kwx=mymax(prop.kwx,3);
512 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)
516 if (prop.hg[j]<0.5 || prop.hg[j]>3000.0)
519 dmin=abs(prop.he[0]-prop.he[1])/200e-3;
520 q=adiff(0.0,prop,propa);
521 xae=pow(prop.wn*prop.gme*prop.gme,-THIRD);
522 d3=mymax(propa.dlsa,1.3787*xae+propa.dla);
524 a3=adiff(d3,prop,propa);
525 a4=adiff(d4,prop,propa);
526 propa.emd=(a4-a3)/(d4-d3);
527 propa.aed=a3-propa.emd*d3;
538 if ((prop.dist>0.0) || (prop.mdp<0.0))
540 if (prop.dist>1000e3)
541 prop.kwx=mymax(prop.kwx,1);
544 prop.kwx=mymax(prop.kwx,3);
546 if (prop.dist<1e3 || prop.dist>2000e3)
557 if (prop.dist>1000e3)
558 prop.kwx=mymax(prop.kwx,1);
561 prop.kwx=mymax(prop.kwx,3);
563 if (prop.dist<1e3 || prop.dist>2000e3)
568 if (prop.dist<propa.dlsa)
572 /* Coefficients for the line-of-sight range */
574 q=alos(0.0,prop,propa);
576 a2=propa.aed+d2*propa.emd;
577 d0=1.908*prop.wn*prop.he[0]*prop.he[1];
581 d0=mymin(d0,0.5*propa.dla);
582 d1=d0+0.25*(propa.dla-d0);
586 d1=mymax(-propa.aed/propa.emd,0.25*propa.dla);
588 a1=alos(d1,prop,propa);
592 a0=alos(d0,prop,propa);
594 propa.ak2=mymax(0.0,((d2-d0)*(a1-a0)-(d1-d0)*(a2-a0))/((d2-d0)*log(d1-d0)-(d1-d0)*q));
596 if (propa.ak2<=0.0 && propa.aed<0.0)
599 propa.ak1=(a2-a1)/(d2-d1);
607 propa.ak1=(a2-a0-propa.ak2*q)/(d2-d0);
612 propa.ak2=FORTRAN_DIM(a2,a0)/q;
618 propa.ael=a2-propa.ak1*d2-propa.ak2*log(d2);
624 prop.aref=propa.ael+propa.ak1*prop.dist+propa.ak2*log(prop.dist);
631 q=ascat(0.0,prop,propa);
634 a6=ascat(d6,prop,propa);
635 a5=ascat(d5,prop,propa);
646 propa.ems=(a6-a5)/200e3;
647 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)));
648 propa.aes=(propa.emd-propa.ems)*propa.dx+propa.aed;
654 if (prop.dist<=propa.dx)
655 prop.aref=propa.aed+propa.emd*prop.dist;
657 prop.aref=propa.aes+propa.ems*prop.dist;
660 prop.aref=FORTRAN_DIM(prop.aref,0.0);
663 void lrprop (double d, prop_type &prop, propa_type &propa)
666 static bool wlos, wscat;
667 static double dmin, xae;
668 complex<double> prop_zgnd(prop.zgndreal,prop.zgndimag);
669 double a0, a1, a2, a3, a4, a5, a6;
670 double d0, d1, d2, d3, d4, d5, d6;
678 propa.dls[j]=sqrt(2.0*prop.he[j]/prop.gme);
680 propa.dlsa=propa.dls[0]+propa.dls[1];
681 propa.dla=prop.dl[0]+prop.dl[1];
682 propa.tha=mymax(prop.the[0]+prop.the[1],-propa.dla*prop.gme);
686 if (prop.wn<0.838 || prop.wn>210.0)
687 prop.kwx=mymax(prop.kwx,1);
690 if (prop.hg[j]<1.0 || prop.hg[j]>1000.0)
691 prop.kwx=mymax(prop.kwx,1);
694 if (abs(prop.the[j]) >200e-3 || prop.dl[j]<0.1*propa.dls[j] || prop.dl[j]>3.0*propa.dls[j] )
695 prop.kwx=mymax(prop.kwx,3);
697 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)
701 if (prop.hg[j]<0.5 || prop.hg[j]>3000.0)
704 dmin=abs(prop.he[0]-prop.he[1])/200e-3;
705 q=adiff(0.0,prop,propa);
706 /* xae=pow(prop.wn*pow(prop.gme,2.),-THIRD); -- JDM made argument 2 a double */
707 xae=pow(prop.wn*(prop.gme*prop.gme),-THIRD); /* No 2nd pow() */
708 d3=mymax(propa.dlsa,1.3787*xae+propa.dla);
710 a3=adiff(d3,prop,propa);
711 a4=adiff(d4,prop,propa);
712 propa.emd=(a4-a3)/(d4-d3);
713 propa.aed=a3-propa.emd*d3;
724 if (prop.dist>1000e3)
725 prop.kwx=mymax(prop.kwx,1);
728 prop.kwx=mymax(prop.kwx,3);
730 if (prop.dist<1e3 || prop.dist>2000e3)
734 if (prop.dist<propa.dlsa)
738 q=alos(0.0,prop,propa);
740 a2=propa.aed+d2*propa.emd;
741 d0=1.908*prop.wn*prop.he[0]*prop.he[1];
745 d0=mymin(d0,0.5*propa.dla);
746 d1=d0+0.25*(propa.dla-d0);
750 d1=mymax(-propa.aed/propa.emd,0.25*propa.dla);
752 a1=alos(d1,prop,propa);
757 a0=alos(d0,prop,propa);
759 propa.ak2=mymax(0.0,((d2-d0)*(a1-a0)-(d1-d0)*(a2-a0))/((d2-d0)*log(d1/d0)-(d1-d0)*q));
760 wq=propa.aed>=0.0 || propa.ak2>0.0;
764 propa.ak1=(a2-a0-propa.ak2*q)/(d2-d0);
769 propa.ak2=FORTRAN_DIM(a2,a0)/q;
779 propa.ak1=(a2-a1)/(d2-d1);
788 propa.ak1=(a2-a1)/(d2-d1);
795 propa.ael=a2-propa.ak1*d2-propa.ak2*log(d2);
800 prop.aref=propa.ael+propa.ak1*prop.dist+propa.ak2*log(prop.dist);
804 if (prop.dist<=0.0 || prop.dist>=propa.dlsa)
808 q=ascat(0.0,prop,propa);
811 a6=ascat(d6,prop,propa);
812 a5=ascat(d5,prop,propa);
816 propa.ems=(a6-a5)/200e3;
817 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)));
818 propa.aes=(propa.emd-propa.ems)*propa.dx+propa.aed;
831 if (prop.dist>propa.dx)
832 prop.aref=propa.aes+propa.ems*prop.dist;
834 prop.aref=propa.aed+propa.emd*prop.dist;
837 prop.aref=mymax(prop.aref,0.0);
840 double curve (double const &c1, double const &c2, double const &x1,
841 double const &x2, double const &x3, double const &de)
843 /* return (c1+c2/(1.0+pow((de-x2)/x3,2.0)))*pow(de/x1,2.0)/(1.0+pow(de/x1,2.0)); */
852 return (c1+c2/(1.0+temp1))*temp2/(1.0+temp2);
855 double avar(double zzt, double zzl, double zzc, prop_type &prop, propv_type &propv)
858 static double dexa, de, vmd, vs0, sgl, sgtm, sgtp, sgtd, tgtd,
859 gm, gp, cv1, cv2, yv1, yv2, yv3, csm1, csm2, ysm1, ysm2,
860 ysm3, csp1, csp2, ysp1, ysp2, ysp3, csd1, zd, cfm1, cfm2,
861 cfm3, cfp1, cfp2, cfp3;
863 double bv1[7]={-9.67,-0.62,1.26,-9.21,-0.62,-0.39,3.15};
864 double bv2[7]={12.7,9.19,15.5,9.05,9.19,2.86,857.9};
865 double xv1[7]={144.9e3,228.9e3,262.6e3,84.1e3,228.9e3,141.7e3,2222.e3};
866 double xv2[7]={190.3e3,205.2e3,185.2e3,101.1e3,205.2e3,315.9e3,164.8e3};
867 double xv3[7]={133.8e3,143.6e3,99.8e3,98.6e3,143.6e3,167.4e3,116.3e3};
868 double bsm1[7]={2.13,2.66,6.11,1.98,2.68,6.86,8.51};
869 double bsm2[7]={159.5,7.67,6.65,13.11,7.16,10.38,169.8};
870 double xsm1[7]={762.2e3,100.4e3,138.2e3,139.1e3,93.7e3,187.8e3,609.8e3};
871 double xsm2[7]={123.6e3,172.5e3,242.2e3,132.7e3,186.8e3,169.6e3,119.9e3};
872 double xsm3[7]={94.5e3,136.4e3,178.6e3,193.5e3,133.5e3,108.9e3,106.6e3};
873 double bsp1[7]={2.11,6.87,10.08,3.68,4.75,8.58,8.43};
874 double bsp2[7]={102.3,15.53,9.60,159.3,8.12,13.97,8.19};
875 double xsp1[7]={636.9e3,138.7e3,165.3e3,464.4e3,93.2e3,216.0e3,136.2e3};
876 double xsp2[7]={134.8e3,143.7e3,225.7e3,93.1e3,135.9e3,152.0e3,188.5e3};
877 double xsp3[7]={95.6e3,98.6e3,129.7e3,94.2e3,113.4e3,122.7e3,122.9e3};
878 double bsd1[7]={1.224,0.801,1.380,1.000,1.224,1.518,1.518};
879 double bzd1[7]={1.282,2.161,1.282,20.,1.282,1.282,1.282};
880 double bfm1[7]={1.0,1.0,1.0,1.0,0.92,1.0,1.0};
881 double bfm2[7]={0.0,0.0,0.0,0.0,0.25,0.0,0.0};
882 double bfm3[7]={0.0,0.0,0.0,0.0,1.77,0.0,0.0};
883 double bfp1[7]={1.0,0.93,1.0,0.93,0.93,1.0,1.0};
884 double bfp2[7]={0.0,0.31,0.0,0.19,0.31,0.0,0.0};
885 double bfp3[7]={0.0,2.00,0.0,1.79,2.00,0.0,0.0};
887 double rt=7.8, rl=24.0, avarv, q, vs, zt, zl, zc;
888 double sgt, yr, temp1, temp2;
889 int temp_klim=propv.klim-1;
896 if (propv.klim<=0 || propv.klim>7)
900 prop.kwx=mymax(prop.kwx,2);
908 csm1=bsm1[temp_klim];
909 csm2=bsm2[temp_klim];
910 ysm1=xsm1[temp_klim];
911 ysm2=xsm2[temp_klim];
912 ysm3=xsm3[temp_klim];
913 csp1=bsp1[temp_klim];
914 csp2=bsp2[temp_klim];
915 ysp1=xsp1[temp_klim];
916 ysp2=xsp2[temp_klim];
917 ysp3=xsp3[temp_klim];
918 csd1=bsd1[temp_klim];
920 cfm1=bfm1[temp_klim];
921 cfm2=bfm2[temp_klim];
922 cfm3=bfm3[temp_klim];
923 cfp1=bfp1[temp_klim];
924 cfp2=bfp2[temp_klim];
925 cfp3=bfp3[temp_klim];
942 prop.kwx=mymax(prop.kwx,2);
946 q=log(0.133*prop.wn);
948 /* gm=cfm1+cfm2/(pow(cfm3*q,2.0)+1.0); */
949 /* gp=cfp1+cfp2/(pow(cfp3*q,2.0)+1.0); */
951 gm=cfm1+cfm2/((cfm3*q*cfm3*q)+1.0);
952 gp=cfp1+cfp2/((cfp3*q*cfp3*q)+1.0);
955 dexa=sqrt(18e6*prop.he[0])+sqrt(18e6*prop.he[1])+pow((575.7e12/prop.wn),THIRD);
959 de=130e3*prop.dist/dexa;
961 de=130e3+prop.dist-dexa;
964 vmd=curve(cv1,cv2,yv1,yv2,yv3,de);
965 sgtm=curve(csm1,csm2,ysm1,ysm2,ysm3,de)*gm;
966 sgtp=curve(csp1,csp2,ysp1,ysp2,ysp3,de)*gp;
974 q=(1.0-0.8*exp(-prop.dist/50e3))*prop.dh*prop.wn;
982 /* vs0=pow(5.0+3.0*exp(-de/100e3),2.0); */
983 temp1=(5.0+3.0*exp(-de/100e3));
1010 if (fabs(zt)>3.1 || fabs(zl)>3.1 || fabs(zc)>3.1)
1011 prop.kwx=mymax(prop.kwx,1);
1022 /* vs=vs0+pow(sgt*zt,2.0)/(rt+zc*zc)+pow(sgl*zl,2.0)/(rl+zc*zc); */
1027 vs=vs0+(temp1*temp1)/(rt+zc*zc)+(temp2*temp2)/(rl+zc*zc);
1032 propv.sgc=sqrt(sgt*sgt+sgl*sgl+vs);
1038 propv.sgc=sqrt(sgl*sgl+vs);
1043 yr=sqrt(sgt*sgt+sgl*sgl)*zt;
1053 avarv=prop.aref-vmd-yr-propv.sgc*zc;
1056 avarv=avarv*(29.0-avarv)/(29.0-10.0*avarv);
1061 void hzns (double pfl[], prop_type &prop)
1065 double xi, za, zb, qc, q, sb, sa;
1069 za=pfl[2]+prop.hg[0];
1070 zb=pfl[np+2]+prop.hg[1];
1073 prop.the[1]=(zb-za)/prop.dist;
1074 prop.the[0]=prop.the[1]-q;
1075 prop.the[1]=-prop.the[1]-q;
1076 prop.dl[0]=prop.dist;
1077 prop.dl[1]=prop.dist;
1085 for (int i=1; i<np; i++)
1089 q=pfl[i+2]-(qc*sa+prop.the[0])*sa-za;
1100 q=pfl[i+2]-(qc*sb+prop.the[1])*sb-zb;
1112 void z1sq1 (double z[], const double &x1, const double &x2, double& z0, double& zn)
1114 double xn, xa, xb, x, a, b;
1118 xa=int(FORTRAN_DIM(x1/z[1],0.0));
1119 xb=xn-int(FORTRAN_DIM(xn,x2/z[1]));
1123 xa=FORTRAN_DIM(xa,1.0);
1124 xb=xn-FORTRAN_DIM(xn,xb+1.0);
1133 a=0.5*(z[ja+2]+z[jb+2]);
1134 b=0.5*(z[ja+2]-z[jb+2])*x;
1136 for (int i=2; i<=n; ++i)
1145 b=b*12.0/((xa*xa+2.0)*xa);
1150 double qtile (const int &nn, double a[], const int &ir)
1152 double q=0.0, r; /* q initialization -- KD2BD */
1153 int m, n, i, j, j1=0, i0=0, k; /* more initializations -- KD2BD */
1159 k=mymin(mymax(0,ir),n);
1172 while (i<=n && a[i]>=q)
1180 while (j>=m && a[j]<=q)
1219 double qerf(const double &z)
1221 double b1=0.319381530, b2=-0.356563782, b3=1.781477937;
1222 double b4=-1.821255987, b5=1.330274429;
1223 double rp=4.317008, rrt2pi=0.398942280;
1234 qerfv=exp(-0.5*x*x)*rrt2pi*((((b5*t+b4)*t+b3)*t+b2)*t+b1)*t;
1243 double d1thx(double pfl[], const double &x1, const double &x2)
1245 int np, ka, kb, n, k, j;
1246 double d1thxv, sn, xa, xb;
1254 if (xb-xa<2.0) // exit out
1257 ka=(int)(0.1*(xb-xa+8.0));
1258 ka=mymin(mymax(4,ka),25);
1262 assert((s=new double[n+2])!=0);
1271 while (xa>0.0 && k<np)
1277 s[j+2]=pfl[k+2]+(pfl[k+2]-pfl[k+1])*xa;
1281 z1sq1(s,0.0,sn,xa,xb);
1290 d1thxv=qtile(n-1,s+2,ka-1)-qtile(n-1,s+2,kb-1);
1291 d1thxv/=1.0-0.8*exp(-(x2-x1)/50.0e3);
1297 void qlrpfl(double pfl[], int klimx, int mdvarx, prop_type &prop, propa_type &propa, propv_type &propv)
1300 double xl[2], q, za, zb, temp;
1302 prop.dist=pfl[0]*pfl[1];
1307 xl[j]=mymin(15.0*prop.hg[j],0.1*prop.dl[j]);
1309 xl[1]=prop.dist-xl[1];
1310 prop.dh=d1thx(pfl,xl[0],xl[1]);
1312 if (prop.dl[0]+prop.dl[1]>1.5*prop.dist)
1314 z1sq1(pfl,xl[0],xl[1],za,zb);
1315 prop.he[0]=prop.hg[0]+FORTRAN_DIM(pfl[2],za);
1316 prop.he[1]=prop.hg[1]+FORTRAN_DIM(pfl[np+2],zb);
1319 prop.dl[j]=sqrt(2.0*prop.he[j]/prop.gme)*exp(-0.07*sqrt(prop.dh/mymax(prop.he[j],5.0)));
1321 q=prop.dl[0]+prop.dl[1];
1325 /* q=pow(prop.dist/q,2.0); */
1332 prop.dl[j]=sqrt(2.0*prop.he[j]/prop.gme)*exp(-0.07*sqrt(prop.dh/mymax(prop.he[j],5.0)));
1338 q=sqrt(2.0*prop.he[j]/prop.gme);
1339 prop.the[j]=(0.65*prop.dh*(q/prop.dl[j]-1.0)-2.0*prop.he[j])/q;
1345 z1sq1(pfl,xl[0],0.9*prop.dl[0],za,q);
1346 z1sq1(pfl,prop.dist-0.9*prop.dl[1],xl[1],q,zb);
1347 prop.he[0]=prop.hg[0]+FORTRAN_DIM(pfl[2],za);
1348 prop.he[1]=prop.hg[1]+FORTRAN_DIM(pfl[np+2],zb);
1352 propv.lvar=mymax(propv.lvar,3);
1357 propv.lvar=mymax(propv.lvar,4);
1366 lrprop(0.0,prop,propa);
1369 double deg2rad(double d)
1371 return d*3.1415926535897/180.0;
1374 //********************************************************
1375 //* Point-To-Point Mode Calculations *
1376 //********************************************************
1378 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)
1380 /******************************************************************************
1383 0-Horizontal, 1-Vertical
1386 1-Equatorial, 2-Continental Subtropical,
1387 3-Maritime Tropical, 4-Desert, 5-Continental Temperate,
1388 6-Maritime Temperate, Over Land, 7-Maritime Temperate,
1391 conf, rel: .01 to .99
1393 elev[]: [num points - 1], [delta dist(meters)],
1394 [height(meters) point 1], ..., [height(meters) point n]
1396 errnum: 0- No Error.
1397 1- Warning: Some parameters are nearly out of range.
1398 Results should be used with caution.
1399 2- Note: Default parameters have been substituted for
1401 3- Warning: A combination of parameters is out of range.
1402 Results are probably invalid.
1403 Other- Warning: Some parameters are out of range.
1404 Results are probably invalid.
1406 *****************************************************************************/
1413 double eno, enso, q;
1420 propv.klim=radio_climate;
1427 dkm=(elev[1]*elev[0])/1000.0;
1435 ja=(long)(3.0+0.1*elev[0]); /* KD2BD added (long) */
1438 for (i=ja-1; i<jb; ++i)
1446 qlrps(frq_mhz,zsys,q,pol,eps_dielect,sgm_conductivity,prop);
1447 qlrpfl(elev,propv.klim,propv.mdvar,prop,propa,propv);
1448 fs=32.45+20.0*log10(frq_mhz)+20.0*log10(prop.dist/1000.0);
1449 q=prop.dist-propa.dla;
1452 strcpy(strmode,"Line-Of-Sight Mode");
1456 strcpy(strmode,"Single Horizon");
1458 else if (int(q)>0.0)
1459 strcpy(strmode,"Double Horizon");
1461 if (prop.dist<=propa.dlsa || prop.dist <= propa.dx)
1462 strcat(strmode,", Diffraction Dominant");
1464 else if (prop.dist>propa.dx)
1465 strcat(strmode, ", Troposcatter Dominant");
1468 dbloss=avar(zr,0.0,zc,prop,propv)+fs;
1472 //********************************************************
1473 //* Area Mode Calculations *
1474 //********************************************************
1476 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)
1478 // pol: 0-Horizontal, 1-Vertical
1479 // TSiteCriteria, RSiteCriteria:
1480 // 0 - random, 1 - careful, 2 - very careful
1482 // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1483 // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1484 // 7-Maritime Temperate, Over Sea
1485 // ModVar: 0 - Single: pctConf is "Time/Situation/Location", pctTime, pctLoc not used
1486 // 1 - Individual: pctTime is "Situation/Location", pctConf is "Confidence", pctLoc not used
1487 // 2 - Mobile: pctTime is "Time/Locations (Reliability)", pctConf is "Confidence", pctLoc not used
1488 // 3 - Broadcast: pctTime is "Time", pctLoc is "Location", pctConf is "Confidence"
1489 // pctTime, pctLoc, pctConf: .01 to .99
1490 // errnum: 0- No Error.
1491 // 1- Warning: Some parameters are nearly out of range.
1492 // Results should be used with caution.
1493 // 2- Note: Default parameters have been substituted for impossible ones.
1494 // 3- Warning: A combination of parameters is out of range.
1495 // Results are probably invalid.
1496 // Other- Warning: Some parameters are out of range.
1497 // Results are probably invalid.
1498 // NOTE: strmode is not used at this time.
1503 double zt, zl, zc, xlb;
1506 double eps, eno, sgm;
1510 kst[0]=(int)TSiteCriteria;
1511 kst[1]=(int)RSiteCriteria;
1512 zt=qerfi(pctTime/100.0);
1513 zl=qerfi(pctLoc/100.0);
1514 zc=qerfi(pctConf/100.0);
1516 sgm=sgm_conductivity;
1521 /* propv.klim = (__int32) radio_climate; -KD2BD replaced below */
1522 propv.klim=(long)radio_climate;
1527 qlrps(frq_mhz, 0.0, eno, ipol, eps, sgm, prop);
1528 qlra(kst, propv.klim, ivar, prop, propv);
1533 lrprop(dist_km*1000.0, prop, propa);
1534 fs=32.45+20.0*log10(frq_mhz)+20.0*log10(prop.dist/1000.0);
1535 xlb=fs+avar(zt, zl, zc, prop, propv);