From 87d6b4cd2264dc6bac658a0885d6e942588da40f Mon Sep 17 00:00:00 2001 From: Cloud-RF Date: Sat, 6 Apr 2013 22:07:29 +0100 Subject: [PATCH] first commit --- itm.cpp | 1532 +++++++++++++++++++ main.cpp | 4287 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 5819 insertions(+) create mode 100644 itm.cpp create mode 100644 main.cpp diff --git a/itm.cpp b/itm.cpp new file mode 100644 index 0000000..486dcb0 --- /dev/null +++ b/itm.cpp @@ -0,0 +1,1532 @@ +/*****************************************************************************\ + * * + * The following code was derived from public domain ITM code available * + * at ftp://flattop.its.bldrdoc.gov/itm/ITMDLL.cpp that was released on * + * June 26, 2007. It was modified to remove Microsoft Windows "dll-isms", * + * redundant and unnecessary #includes, redundant and unnecessary { }'s, * + * to initialize uninitialized variables, type cast some variables, * + * re-format the code for easier reading, and to replace pow() function * + * calls with explicit multiplications wherever possible to increase * + * execution speed and improve computational accuracy. * + * * +\*****************************************************************************/ + + +// ************************************* +// C++ routines for this program are taken from +// a translation of the FORTRAN code written by +// U.S. Department of Commerce NTIA/ITS +// Institute for Telecommunication Sciences +// ***************** +// Irregular Terrain Model (ITM) (Longley-Rice) +// ************************************* + +#include +#include +#include +#include + +#define THIRD (1.0/3.0) + +using namespace std; + +struct tcomplex +{ + double tcreal; + double tcimag; +}; + +struct prop_type +{ + double aref; + double dist; + double hg[2]; + double wn; + double dh; + double ens; + double gme; + double zgndreal; + double zgndimag; + double he[2]; + double dl[2]; + double the[2]; + int kwx; + int mdp; +}; + +struct propv_type +{ + double sgc; + int lvar; + int mdvar; + int klim; +}; + +struct propa_type +{ + double dlsa; + double dx; + double ael; + double ak1; + double ak2; + double aed; + double emd; + double aes; + double ems; + double dls[2]; + double dla; + double tha; +}; + +int mymin(const int &i, const int &j) +{ + if (ij) + return i; + else + return j; +} + +double mymin(const double &a, const double &b) +{ + if (ab) + return a; + else + return b; +} + +double FORTRAN_DIM(const double &x, const double &y) +{ + // This performs the FORTRAN DIM function. + // result is x-y if x is greater than y; otherwise result is 0.0 + + if (x>y) + return x-y; + else + return 0.0; +} + +double aknfe(const double &v2) +{ + double a; + + if (v2<5.76) + a=6.02+9.11*sqrt(v2)-1.27*v2; + else + a=12.953+4.343*log(v2); + + return a; +} + +double fht(const double& x, const double& pk) +{ + double w, fhtv; + if (x<200.0) + { + w=-log(pk); + + /* if (pk < 1e-5 || x*pow(w,3.0) > 5495.0 ) */ + + if (pk < 1e-5 || (x*w*w*w) > 5495.0 ) + { + fhtv=-117.0; + + if (x>1.0) + fhtv=17.372*log(x)+fhtv; + } + + else + fhtv=2.5e-5*x*x/pk-8.686*w-15.0; + } + + else + { + fhtv=0.05751*x-4.343*log(x); + + if (x<2000.0) + { + w=0.0134*x*exp(-0.005*x); + fhtv=(1.0-w)*fhtv+w*(17.372*log(x)-117.0); + } + } + + return fhtv; +} + +double h0f(double r, double et) +{ + double a[5]={25.0, 80.0, 177.0, 395.0, 705.0}; + double b[5]={24.0, 45.0, 68.0, 80.0, 105.0}; + double q, x; + int it; + double h0fv; + + it=(int)et; + + if (it<=0) + { + it=1; + q=0.0; + } + + else if (it>=5) + { + it=5; + q=0.0; + } + + else + q=et-it; + + /* x=pow(1.0/r,2.0); */ + x=(1.0/r); + x*=x; + h0fv=4.343*log((a[it-1]*x+b[it-1])*x+1.0); + + if (q!=0.0) + h0fv=(1.0-q)*h0fv+q*4.343*log((a[it]*x+b[it])*x+1.0); + + return h0fv; +} + +double ahd(double td) +{ + int i; + double a[3] = { 133.4, 104.6, 71.8}; + double b[3] = {0.332e-3, 0.212e-3, 0.157e-3}; + double c[3] = { -4.343, -1.086, 2.171}; + + if (td<=10e3) + i=0; + + else if (td<=70e3) + i=1; + + else + i=2; + + return a[i]+b[i]*td+c[i]*log(td); +} + +double adiff( double d, prop_type &prop, propa_type &propa) +{ + complex prop_zgnd(prop.zgndreal,prop.zgndimag); + static double wd1, xd1, afo, qk, aht, xht; + double a, q, pk, ds, th, wa, ar, wd, adiffv; + + if (d==0) + { + q=prop.hg[0]*prop.hg[1]; + qk=prop.he[0]*prop.he[1]-q; + + if (prop.mdp<0.0) + q+=10.0; + + wd1=sqrt(1.0+qk/q); + xd1=propa.dla+propa.tha/prop.gme; + q=(1.0-0.8*exp(-propa.dlsa/50e3))*prop.dh; + q*=0.78*exp(-pow(q/16.0,0.25)); + afo=mymin(15.0,2.171*log(1.0+4.77e-4*prop.hg[0]*prop.hg[1]*prop.wn*q)); + qk=1.0/abs(prop_zgnd); + aht=20.0; + xht=0.0; + + for (int j=0; j<2; ++j) + { + /* a=0.5*pow(prop.dl[j],2.0)/prop.he[j]; */ + a=0.5*(prop.dl[j]*prop.dl[j])/prop.he[j]; + wa=pow(a*prop.wn,THIRD); + pk=qk/wa; + q=(1.607-pk)*151.0*wa*prop.dl[j]/a; + xht+=q; + aht+=fht(q,pk); + } + + adiffv=0.0; + } + + else + { + th=propa.tha+d*prop.gme; + ds=d-propa.dla; + /* q=0.0795775*prop.wn*ds*pow(th,2.0); */ + q=0.0795775*prop.wn*ds*th*th; + adiffv=aknfe(q*prop.dl[0]/(ds+prop.dl[0]))+aknfe(q*prop.dl[1]/(ds+prop.dl[1])); + a=ds/th; + wa=pow(a*prop.wn,THIRD); + pk=qk/wa; + q=(1.607-pk)*151.0*wa*th+xht; + ar=0.05751*q-4.343*log(q)-aht; + q=(wd1+xd1/d)*mymin(((1.0-0.8*exp(-d/50e3))*prop.dh*prop.wn),6283.2); + wd=25.1/(25.1+sqrt(q)); + adiffv=ar*wd+(1.0-wd)*adiffv+afo; + } + + return adiffv; +} + +double ascat( double d, prop_type &prop, propa_type &propa) +{ + static double ad, rr, etq, h0s; + double h0, r1, r2, z0, ss, et, ett, th, q; + double ascatv, temp; + + if (d==0.0) + { + ad=prop.dl[0]-prop.dl[1]; + rr=prop.he[1]/prop.he[0]; + + if (ad<0.0) + { + ad=-ad; + rr=1.0/rr; + } + + etq=(5.67e-6*prop.ens-2.32e-3)*prop.ens+0.031; + h0s=-15.0; + ascatv=0.0; + } + + else + { + if (h0s>15.0) + h0=h0s; + else + { + th=prop.the[0]+prop.the[1]+d*prop.gme; + r2=2.0*prop.wn*th; + r1=r2*prop.he[0]; + r2*=prop.he[1]; + + if (r1<0.2 && r2<0.2) + return 1001.0; // <==== early return + + ss=(d-ad)/(d+ad); + q=rr/ss; + ss=mymax(0.1,ss); + q=mymin(mymax(0.1,q),10.0); + z0=(d-ad)*(d+ad)*th*0.25/d; + /* et=(etq*exp(-pow(mymin(1.7,z0/8.0e3),6.0))+1.0)*z0/1.7556e3; */ + temp=mymin(1.7,z0/8.0e3); + temp=temp*temp*temp*temp*temp*temp; + et=(etq*exp(-temp)+1.0)*z0/1.7556e3; + + ett=mymax(et,1.0); + h0=(h0f(r1,ett)+h0f(r2,ett))*0.5; + h0+=mymin(h0,(1.38-log(ett))*log(ss)*log(q)*0.49); + h0=FORTRAN_DIM(h0,0.0); + + if (et<1.0) + 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)); + + if (h0>15.0 && h0s>=0.0) + h0=h0s; + } + + h0s=h0; + th=propa.tha+d*prop.gme; + /* 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; */ + + 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; + + } + + return ascatv; +} + +double qerfi (double q) +{ + double x, t, v; + double c0=2.515516698; + double c1=0.802853; + double c2=0.010328; + double d1=1.432788; + double d2=0.189269; + double d3=0.001308; + + x=0.5-q; + t=mymax(0.5-fabs(x),0.000001); + t=sqrt(-2.0*log(t)); + v=t-((c2*t+c1)*t+c0)/(((d3*t+d2)*t+d1)*t+1.0); + + if (x<0.0) + v=-v; + + return v; +} + +void qlrps( double fmhz, double zsys, double en0, int ipol, double eps, double sgm, prop_type &prop) + +{ + double gma=157e-9; + + prop.wn=fmhz/47.7; + prop.ens=en0; + + if (zsys!=0.0) + prop.ens*=exp(-zsys/9460.0); + + prop.gme=gma*(1.0-0.04665*exp(prop.ens/179.3)); + complex zq, prop_zgnd(prop.zgndreal,prop.zgndimag); + zq=complex (eps,376.62*sgm/prop.wn); + prop_zgnd=sqrt(zq-1.0); + + if (ipol!=0.0) + prop_zgnd=prop_zgnd/zq; + + prop.zgndreal=prop_zgnd.real(); + prop.zgndimag=prop_zgnd.imag(); +} + +double abq_alos(complex r) +{ + return r.real()*r.real()+r.imag()*r.imag(); +} + +double alos(double d, prop_type &prop, propa_type &propa) +{ + complex prop_zgnd(prop.zgndreal,prop.zgndimag); + static double wls; + complex r; + double s, sps, q; + double alosv; + + if (d==0.0) + { + wls=0.021/(0.021+prop.wn*prop.dh/mymax(10e3,propa.dlsa)); + alosv=0.0; + } + + else + { + q=(1.0-0.8*exp(-d/50e3))*prop.dh; + s=0.78*q*exp(-pow(q/16.0,0.25)); + q=prop.he[0]+prop.he[1]; + sps=q/sqrt(d*d+q*q); + r=(sps-prop_zgnd)/(sps+prop_zgnd)*exp(-mymin(10.0,prop.wn*s*sps)); + q=abq_alos(r); + + if (q<0.25 || q1.57) + q=3.14-2.4649/q; + + alosv=(-4.343*log(abq_alos(complex(cos(q),-sin(q))+r))-alosv)*wls+alosv; + } + + return alosv; +} + + +void qlra(int kst[], int klimx, int mdvarx, prop_type &prop, propv_type &propv) +{ + double q; + + for (int j=0; j<2; ++j) + { + if (kst[j]<=0) + prop.he[j]=prop.hg[j]; + else + { + q=4.0; + + if (kst[j]!=1) + q=9.0; + + if (prop.hg[j]<5.0) + q*=sin(0.3141593*prop.hg[j]); + + prop.he[j]=prop.hg[j]+(1.0+q)*exp(-mymin(20.0,2.0*prop.hg[j]/mymax(1e-3,prop.dh))); + } + + q=sqrt(2.0*prop.he[j]/prop.gme); + prop.dl[j]=q*exp(-0.07*sqrt(prop.dh/mymax(prop.he[j],5.0))); + prop.the[j]=(0.65*prop.dh*(q/prop.dl[j]-1.0)-2.0*prop.he[j])/q; + } + + prop.mdp=1; + propv.lvar=mymax(propv.lvar,3); + + if (mdvarx>=0) + { + propv.mdvar=mdvarx; + propv.lvar=mymax(propv.lvar,4); + } + + if (klimx>0) + { + propv.klim=klimx; + propv.lvar=5; + } +} + +void lrprop (double d, prop_type &prop, propa_type &propa) // PaulM_lrprop +{ + static bool wlos, wscat; + static double dmin, xae; + complex prop_zgnd(prop.zgndreal,prop.zgndimag); + double a0, a1, a2, a3, a4, a5, a6; + double d0, d1, d2, d3, d4, d5, d6; + bool wq; + double q; + int j; + + if (prop.mdp!=0) + { + for (j=0; j<2; j++) + propa.dls[j]=sqrt(2.0*prop.he[j]/prop.gme); + + propa.dlsa=propa.dls[0]+propa.dls[1]; + propa.dla=prop.dl[0]+prop.dl[1]; + propa.tha=mymax(prop.the[0]+prop.the[1],-propa.dla*prop.gme); + wlos=false; + wscat=false; + + if (prop.wn<0.838 || prop.wn>210.0) + prop.kwx=mymax(prop.kwx,1); + + for (j=0; j<2; j++) + if (prop.hg[j]<1.0 || prop.hg[j]>1000.0) + prop.kwx=mymax(prop.kwx,1); + + for (j=0; j<2; j++) + if (abs(prop.the[j]) >200e-3 || prop.dl[j]<0.1*propa.dls[j] || prop.dl[j]>3.0*propa.dls[j]) + prop.kwx=mymax(prop.kwx,3); + + 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) + + prop.kwx=4; + + for (j=0; j<2; j++) + + if (prop.hg[j]<0.5 || prop.hg[j]>3000.0) + { + prop.kwx=4; + } + + dmin=abs(prop.he[0]-prop.he[1])/200e-3; + q=adiff(0.0,prop,propa); + /* xae=pow(prop.wn*pow(prop.gme,2),-THIRD); */ + xae=pow(prop.wn*prop.gme*prop.gme,-THIRD); + d3=mymax(propa.dlsa,1.3787*xae+propa.dla); + d4=d3+2.7574*xae; + a3=adiff(d3,prop,propa); + a4=adiff(d4,prop,propa); + propa.emd=(a4-a3)/(d4-d3); + propa.aed=a3-propa.emd*d3; + } + + if (prop.mdp>=0) + { + prop.mdp=0; + prop.dist=d; + } + + if (prop.dist>0.0) + { + if (prop.dist>1000e3) + prop.kwx=mymax(prop.kwx,1); + + if (prop.dist2000e3) + prop.kwx=4; + } + + if (prop.dist=0.0) + { + d0=mymin(d0,0.5*propa.dla); + d1=d0+0.25*(propa.dla-d0); + } + + else + d1=mymax(-propa.aed/propa.emd,0.25*propa.dla); + + a1=alos(d1,prop,propa); + wq=false; + + if (d0=0.0 || propa.ak2>0.0; + + if (wq) + { + propa.ak1=(a2-a0-propa.ak2*q)/(d2-d0); + + if (propa.ak1<0.0) + { + propa.ak1=0.0; + propa.ak2=FORTRAN_DIM(a2,a0)/q; + + if (propa.ak2==0.0) + propa.ak1=propa.emd; + } + } + } + + if (!wq) + { + propa.ak1=FORTRAN_DIM(a2,a1)/(d2-d1); + propa.ak2=0.0; + + if (propa.ak1==0.0) + propa.ak1=propa.emd; + } + + propa.ael=a2-propa.ak1*d2-propa.ak2*log(d2); + wlos=true; + } + + if (prop.dist>0.0) + prop.aref=propa.ael+propa.ak1*prop.dist+propa.ak2*log(prop.dist); + } + + if (prop.dist<=0.0 || prop.dist>=propa.dlsa) + { + if (!wscat) + { + q=ascat(0.0,prop,propa); + d5=propa.dla+200e3; + d6=d5+200e3; + a6=ascat(d6,prop,propa); + a5=ascat(d5,prop,propa); + + if (a5<1000.0) + { + propa.ems=(a6-a5)/200e3; + 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))); + + propa.aes=(propa.emd-propa.ems)*propa.dx+propa.aed; + } + + else + { + propa.ems=propa.emd; + propa.aes=propa.aed; + propa.dx=10.e6; + } + + wscat=true; + } + + if (prop.dist>propa.dx) + prop.aref=propa.aes+propa.ems*prop.dist; + else + prop.aref=propa.aed+propa.emd*prop.dist; + } + + prop.aref=mymax(prop.aref,0.0); +} + +double curve (double const &c1, double const &c2, double const &x1, double const &x2, double const &x3, double const &de) +{ + /* return (c1+c2/(1.0+pow((de-x2)/x3,2.0)))*pow(de/x1,2.0)/(1.0+pow(de/x1,2.0)); */ + + double temp1, temp2; + + temp1=(de-x2)/x3; + temp2=de/x1; + + temp1*=temp1; + temp2*=temp2; + + return (c1+c2/(1.0+temp1))*temp2/(1.0+temp2); +} + +double avar(double zzt, double zzl, double zzc,prop_type &prop, propv_type &propv) +{ + static int kdv; + 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; + + double bv1[7]={-9.67,-0.62,1.26,-9.21,-0.62,-0.39,3.15}; + double bv2[7]={12.7,9.19,15.5,9.05,9.19,2.86,857.9}; + double xv1[7]={144.9e3,228.9e3,262.6e3,84.1e3,228.9e3,141.7e3,2222.e3}; + double xv2[7]={190.3e3,205.2e3,185.2e3,101.1e3,205.2e3,315.9e3,164.8e3}; + double xv3[7]={133.8e3,143.6e3,99.8e3,98.6e3,143.6e3,167.4e3,116.3e3}; + double bsm1[7]={2.13,2.66,6.11,1.98,2.68,6.86,8.51}; + double bsm2[7]={159.5,7.67,6.65,13.11,7.16,10.38,169.8}; + double xsm1[7]={762.2e3,100.4e3,138.2e3,139.1e3,93.7e3,187.8e3,609.8e3}; + double xsm2[7]={123.6e3,172.5e3,242.2e3,132.7e3,186.8e3,169.6e3,119.9e3}; + double xsm3[7]={94.5e3,136.4e3,178.6e3,193.5e3,133.5e3,108.9e3,106.6e3}; + double bsp1[7]={2.11,6.87,10.08,3.68,4.75,8.58,8.43}; + double bsp2[7]={102.3,15.53,9.60,159.3,8.12,13.97,8.19}; + double xsp1[7]={636.9e3,138.7e3,165.3e3,464.4e3,93.2e3,216.0e3,136.2e3}; + double xsp2[7]={134.8e3,143.7e3,225.7e3,93.1e3,135.9e3,152.0e3,188.5e3}; + double xsp3[7]={95.6e3,98.6e3,129.7e3,94.2e3,113.4e3,122.7e3,122.9e3}; + double bsd1[7]={1.224,0.801,1.380,1.000,1.224,1.518,1.518}; + double bzd1[7]={1.282,2.161,1.282,20.,1.282,1.282,1.282}; + double bfm1[7]={1.0,1.0,1.0,1.0,0.92,1.0,1.0}; + double bfm2[7]={0.0,0.0,0.0,0.0,0.25,0.0,0.0}; + double bfm3[7]={0.0,0.0,0.0,0.0,1.77,0.0,0.0}; + double bfp1[7]={1.0,0.93,1.0,0.93,0.93,1.0,1.0}; + double bfp2[7]={0.0,0.31,0.0,0.19,0.31,0.0,0.0}; + double bfp3[7]={0.0,2.00,0.0,1.79,2.00,0.0,0.0}; + static bool ws, w1; + double rt=7.8, rl=24.0, avarv, q, vs, zt, zl, zc; + double sgt, yr; + int temp_klim=propv.klim-1; + + if (propv.lvar>0) + { + switch (propv.lvar) + { + default: + if (propv.klim<=0 || propv.klim>7) + { + propv.klim = 5; + temp_klim = 4; + prop.kwx=mymax(prop.kwx,2); + } + + cv1 = bv1[temp_klim]; + cv2 = bv2[temp_klim]; + yv1 = xv1[temp_klim]; + yv2 = xv2[temp_klim]; + yv3 = xv3[temp_klim]; + csm1=bsm1[temp_klim]; + csm2=bsm2[temp_klim]; + ysm1=xsm1[temp_klim]; + ysm2=xsm2[temp_klim]; + ysm3=xsm3[temp_klim]; + csp1=bsp1[temp_klim]; + csp2=bsp2[temp_klim]; + ysp1=xsp1[temp_klim]; + ysp2=xsp2[temp_klim]; + ysp3=xsp3[temp_klim]; + csd1=bsd1[temp_klim]; + zd=bzd1[temp_klim]; + cfm1=bfm1[temp_klim]; + cfm2=bfm2[temp_klim]; + cfm3=bfm3[temp_klim]; + cfp1=bfp1[temp_klim]; + cfp2=bfp2[temp_klim]; + cfp3=bfp3[temp_klim]; + + case 4: + kdv=propv.mdvar; + ws=kdv>=20; + + if (ws) + kdv-=20; + + w1=kdv>=10; + + if (w1) + kdv-=10; + + if (kdv<0 || kdv>3) + { + kdv=0; + prop.kwx=mymax(prop.kwx,2); + } + + case 3: + q=log(0.133*prop.wn); + + /* gm=cfm1+cfm2/(pow(cfm3*q,2.0)+1.0); */ + /* gp=cfp1+cfp2/(pow(cfp3*q,2.0)+1.0); */ + + gm=cfm1+cfm2/((cfm3*q*cfm3*q)+1.0); + gp=cfp1+cfp2/((cfp3*q*cfp3*q)+1.0); + + case 2: + dexa=sqrt(18e6*prop.he[0])+sqrt(18e6*prop.he[1])+pow((575.7e12/prop.wn),THIRD); + + case 1: + if (prop.dist3.1 || fabs(zl)>3.1 || fabs(zc)>3.1) + prop.kwx=mymax(prop.kwx,1); + + if (zt<0.0) + sgt=sgtm; + + else if (zt<=zd) + sgt=sgtp; + + else + sgt=sgtd+tgtd/zt; + + /* vs=vs0+pow(sgt*zt,2.0)/(rt+zc*zc)+pow(sgl*zl,2.0)/(rl+zc*zc); */ + vs=vs0+(sgt*zt*sgt*zt)/(rt+zc*zc)+(sgl*zl*sgl*zl)/(rl+zc*zc); + + if (kdv==0) + { + yr=0.0; + propv.sgc=sqrt(sgt*sgt+sgl*sgl+vs); + } + + else if (kdv==1) + { + yr=sgt*zt; + propv.sgc=sqrt(sgl*sgl+vs); + } + + else if (kdv==2) + { + yr=sqrt(sgt*sgt+sgl*sgl)*zt; + propv.sgc=sqrt(vs); + } + + else + { + yr=sgt*zt+sgl*zl; + propv.sgc=sqrt(vs); + } + + avarv=prop.aref-vmd-yr-propv.sgc*zc; + + if (avarv<0.0) + avarv=avarv*(29.0-avarv)/(29.0-10.0*avarv); + + return avarv; +} + +void hzns(double pfl[], prop_type &prop) +{ + bool wq; + int np; + double xi, za, zb, qc, q, sb, sa; + + np=(int)pfl[0]; + xi=pfl[1]; + za=pfl[2]+prop.hg[0]; + zb=pfl[np+2]+prop.hg[1]; + qc=0.5*prop.gme; + q=qc*prop.dist; + prop.the[1]=(zb-za)/prop.dist; + prop.the[0]=prop.the[1]-q; + prop.the[1]=-prop.the[1]-q; + prop.dl[0]=prop.dist; + prop.dl[1]=prop.dist; + + if (np>=2) + { + sa=0.0; + sb=prop.dist; + wq=true; + + for (int i=1; i0.0) + { + prop.the[0]+=q/sa; + prop.dl[0]=sa; + wq=false; + } + + if (!wq) + { + q=pfl[i+2]-(qc*sb+prop.the[1])*sb-zb; + + if (q>0.0) + { + prop.the[1]+=q/sb; + prop.dl[1]=sb; + } + } + } + } +} + +void z1sq1 (double z[], const double &x1, const double &x2, double& z0, double& zn) + +{ + double xn, xa, xb, x, a, b; + int n, ja, jb; + + xn=z[0]; + xa=int(FORTRAN_DIM(x1/z[1],0.0)); + xb=xn-int(FORTRAN_DIM(xn,x2/z[1])); + + if (xb<=xa) + { + xa=FORTRAN_DIM(xa,1.0); + xb=xn-FORTRAN_DIM(xn,xb+1.0); + } + + ja=(int)xa; + jb=(int)xb; + n=jb-ja; + xa=xb-xa; + x=-0.5*xa; + xb+=x; + a=0.5*(z[ja+2]+z[jb+2]); + b=0.5*(z[ja+2]-z[jb+2])*x; + + for (int i=2; i<=n; ++i) + { + ++ja; + x+=1.0; + a+=z[ja+2]; + b+=z[ja+2]*x; + } + + a/=xa; + b=b*12.0/((xa*xa+2.0)*xa); + z0=a-b*xb; + zn=a+b*(xn-xb); +} + +double qtile (const int &nn, double a[], const int &ir) +{ + double q=0.0, r; /* Initializations by KD2BD */ + int m, n, i, j, j1=0, i0=0, k; /* Initializations by KD2BD */ + bool done=false; + bool goto10=true; + + m=0; + n=nn; + k=mymin(mymax(0,ir),n); + + while (!done) + { + if (goto10) + { + q=a[k]; + i0=m; + j1=n; + } + + i=i0; + + while (i<=n && a[i]>=q) + i++; + + if (i>n) + i=n; + j=j1; + + while (j>=m && a[j]<=q) + j--; + + if (jk) + { + a[k]=a[j]; + a[j]=q; + n=j-1; + goto10=true; + } + + else + done=true; + } + + return q; +} + +double qerf(const double &z) +{ + double b1=0.319381530, b2=-0.356563782, b3=1.781477937; + double b4=-1.821255987, b5=1.330274429; + double rp=4.317008, rrt2pi=0.398942280; + double t, x, qerfv; + + x=z; + t=fabs(x); + + if (t>=10.0) + qerfv=0.0; + else + { + t=rp/(t+rp); + qerfv=exp(-0.5*x*x)*rrt2pi*((((b5*t+b4)*t+b3)*t+b2)*t+b1)*t; + } + + if (x<0.0) + qerfv=1.0-qerfv; + + return qerfv; +} + +double d1thx(double pfl[], const double &x1, const double &x2) +{ + int np, ka, kb, n, k, j; + double d1thxv, sn, xa, xb; + double *s; + + np=(int)pfl[0]; + xa=x1/pfl[1]; + xb=x2/pfl[1]; + d1thxv=0.0; + + if (xb-xa<2.0) // exit out + return d1thxv; + + ka=(int)(0.1*(xb-xa+8.0)); + ka=mymin(mymax(4,ka),25); + n=10*ka-5; + kb=n-ka+1; + sn=n-1; + assert((s=new double[n+2])!=0); + s[0]=sn; + s[1]=1.0; + xb=(xb-xa)/sn; + k=(int)(xa+1.0); + xa-=(double)k; + + for (j=0; j0.0 && k1.5*prop.dist) + { + z1sq1(pfl,xl[0],xl[1],za,zb); + prop.he[0]=prop.hg[0]+FORTRAN_DIM(pfl[2],za); + prop.he[1]=prop.hg[1]+FORTRAN_DIM(pfl[np+2],zb); + + for (j=0; j<2; j++) + prop.dl[j]=sqrt(2.0*prop.he[j]/prop.gme)*exp(-0.07*sqrt(prop.dh/mymax(prop.he[j],5.0))); + + q=prop.dl[0]+prop.dl[1]; + + if (q<=prop.dist) + { + /* q=pow(prop.dist/q,2.0); */ + q=((prop.dist/q)*(prop.dist/q)); + + for (j=0; j<2; j++) + { + prop.he[j]*=q; + prop.dl[j]=sqrt(2.0*prop.he[j]/prop.gme)*exp(-0.07*sqrt(prop.dh/mymax(prop.he[j],5.0))); + } + } + + for (j=0; j<2; j++) + { + q=sqrt(2.0*prop.he[j]/prop.gme); + prop.the[j]=(0.65*prop.dh*(q/prop.dl[j]-1.0)-2.0*prop.he[j])/q; + } + } + + else + { + z1sq1(pfl,xl[0],0.9*prop.dl[0],za,q); + z1sq1(pfl,prop.dist-0.9*prop.dl[1],xl[1],q,zb); + prop.he[0]=prop.hg[0]+FORTRAN_DIM(pfl[2],za); + prop.he[1]=prop.hg[1]+FORTRAN_DIM(pfl[np+2],zb); + } + + prop.mdp=-1; + propv.lvar=mymax(propv.lvar,3); + + if (mdvarx>=0) + { + propv.mdvar=mdvarx; + propv.lvar=mymax(propv.lvar,4); + } + + if (klimx>0) + { + propv.klim=klimx; + propv.lvar=5; + } + + lrprop(0.0,prop,propa); +} + +double deg2rad(double d) +{ + return d*3.1415926535897/180.0; +} + +//******************************************************** +//* Point-To-Point Mode Calculations * +//******************************************************** + +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) +{ + // pol: 0-Horizontal, 1-Vertical + // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical, + // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land, + // 7-Maritime Temperate, Over Sea + // conf, rel: .01 to .99 + // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n] + // errnum: 0- No Error. + // 1- Warning: Some parameters are nearly out of range. + // Results should be used with caution. + // 2- Note: Default parameters have been substituted for impossible ones. + // 3- Warning: A combination of parameters is out of range. + // Results are probably invalid. + // Other- Warning: Some parameters are out of range. + // Results are probably invalid. + + prop_type prop; + propv_type propv; + propa_type propa; + + double zsys=0; + double zc, zr; + double eno, enso, q; + long ja, jb, i, np; + double dkm, xkm; + double fs; + + prop.hg[0]=tht_m; + prop.hg[1]=rht_m; + propv.klim=radio_climate; + prop.kwx=0; + propv.lvar=5; + prop.mdp=-1; + zc=qerfi(conf); + zr=qerfi(rel); + np=(long)elev[0]; + dkm=(elev[1]*elev[0])/1000.0; + xkm=elev[1]/1000.0; + eno=eno_ns_surfref; + enso=0.0; + q=enso; + + if (q<=0.0) + { + /* ja = 3.0 + 0.1 * elev[0]; */ + ja=(long)(3.0+0.1* elev[0]); + + jb=np-ja+6; + + for (i=ja-1; i0.0) + strcpy(strmode,"Double Horizon"); + + if (prop.dist<=propa.dlsa || prop.dist <= propa.dx) + strcat(strmode,", Diffraction Dominant"); + + else if (prop.dist>propa.dx) + strcat(strmode, ", Troposcatter Dominant"); + } + + dbloss=avar(zr,0.0,zc,prop,propv)+fs; + errnum=prop.kwx; +} + +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) +{ + // pol: 0-Horizontal, 1-Vertical + // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical, + // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land, + // 7-Maritime Temperate, Over Sea + // timepct, locpct, confpct: .01 to .99 + // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n] + // propmode: Value Mode + // -1 mode is undefined + // 0 Line of Sight + // 5 Single Horizon, Diffraction + // 6 Single Horizon, Troposcatter + // 9 Double Horizon, Diffraction + // 10 Double Horizon, Troposcatter + // errnum: 0- No Error. + // 1- Warning: Some parameters are nearly out of range. + // Results should be used with caution. + // 2- Note: Default parameters have been substituted for impossible ones. + // 3- Warning: A combination of parameters is out of range. + // Results are probably invalid. + // Other- Warning: Some parameters are out of range. + // Results are probably invalid. + + prop_type prop; + propv_type propv; + propa_type propa; + double zsys=0; + double ztime, zloc, zconf; + double eno, enso, q; + long ja, jb, i, np; + double dkm, xkm; + double fs; + + propmode=-1; // mode is undefined + prop.hg[0]=tht_m; + prop.hg[1]=rht_m; + propv.klim=radio_climate; + prop.kwx=0; + propv.lvar=5; + prop.mdp=-1; + ztime=qerfi(timepct); + zloc=qerfi(locpct); + zconf=qerfi(confpct); + + np=(long)elev[0]; + dkm=(elev[1]*elev[0])/1000.0; + xkm=elev[1]/1000.0; + eno=eno_ns_surfref; + enso=0.0; + q=enso; + + if (q<=0.0) + { + /* ja = 3.0 + 0.1 * elev[0]; */ + ja=(long)(3.0+0.1*elev[0]); + jb=np-ja+6; + + for (i=ja-1; i0.0) + propmode=8; // Double Horizon + + if (prop.dist<=propa.dlsa || prop.dist<=propa.dx) + propmode+=1; // Diffraction Dominant + + else if (prop.dist>propa.dx) + propmode+=2; // Troposcatter Dominant + } + + dbloss=avar(ztime, zloc, zconf, prop, propv)+fs; //avar(time,location,confidence) + errnum=prop.kwx; +} + +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) +{ + // pol: 0-Horizontal, 1-Vertical + // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical, + // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land, + // 7-Maritime Temperate, Over Sea + // conf, rel: .01 to .99 + // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n] + // errnum: 0- No Error. + // 1- Warning: Some parameters are nearly out of range. + // Results should be used with caution. + // 2- Note: Default parameters have been substituted for impossible ones. + // 3- Warning: A combination of parameters is out of range. + // Results are probably invalid. + // Other- Warning: Some parameters are out of range. + // Results are probably invalid. + + char strmode[100]; + prop_type prop; + propv_type propv; + propa_type propa; + double zsys=0; + double zc, zr; + double eno, enso, q; + long ja, jb, i, np; + double dkm, xkm; + double fs; + + prop.hg[0]=tht_m; + prop.hg[1]=rht_m; + propv.klim=radio_climate; + prop.kwx=0; + propv.lvar=5; + prop.mdp=-1; + zc=qerfi(conf); + zr=qerfi(rel); + np=(long)elev[0]; + dkm=(elev[1]*elev[0])/1000.0; + xkm=elev[1]/1000.0; + eno=eno_ns_surfref; + enso=0.0; + q=enso; + + if (q<=0.0) + { + /* ja = 3.0 + 0.1 * elev[0]; */ + ja=(long)(3.0+0.1*elev[0]); + + jb=np-ja+6; + + for (i=ja-1; i0.0) + strcpy(strmode,"Double Horizon"); + + if (prop.dist<=propa.dlsa || prop.dist <= propa.dx) + strcat(strmode,", Diffraction Dominant"); + + else if (prop.dist>propa.dx) + strcat(strmode, ", Troposcatter Dominant"); + } + + dbloss=avar(zr,0.0,zc,prop,propv)+fs; //avar(time,location,confidence) + errnum=prop.kwx; +} + + +//******************************************************** +//* Area Mode Calculations * +//******************************************************** + +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) +{ + // pol: 0-Horizontal, 1-Vertical + // TSiteCriteria, RSiteCriteria: + // 0 - random, 1 - careful, 2 - very careful + // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical, + // 4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land, + // 7-Maritime Temperate, Over Sea + // ModVar: 0 - Single: pctConf is "Time/Situation/Location", pctTime, pctLoc not used + // 1 - Individual: pctTime is "Situation/Location", pctConf is "Confidence", pctLoc not used + // 2 - Mobile: pctTime is "Time/Locations (Reliability)", pctConf is "Confidence", pctLoc not used + // 3 - Broadcast: pctTime is "Time", pctLoc is "Location", pctConf is "Confidence" + // pctTime, pctLoc, pctConf: .01 to .99 + // errnum: 0- No Error. + // 1- Warning: Some parameters are nearly out of range. + // Results should be used with caution. + // 2- Note: Default parameters have been substituted for impossible ones. + // 3- Warning: A combination of parameters is out of range. + // Results are probably invalid. + // Other- Warning: Some parameters are out of range. + // Results are probably invalid. + // NOTE: strmode is not used at this time. + + prop_type prop; + propv_type propv; + propa_type propa; + double zt, zl, zc, xlb; + double fs; + long ivar; + double eps, eno, sgm; + long ipol; + int kst[2]; + + kst[0]=(int)TSiteCriteria; + kst[1]=(int)RSiteCriteria; + zt=qerfi(pctTime); + zl=qerfi(pctLoc); + zc=qerfi(pctConf); + eps=eps_dielect; + sgm=sgm_conductivity; + eno=eno_ns_surfref; + prop.dh=deltaH; + prop.hg[0]=tht_m; + prop.hg[1]=rht_m; + /* propv.klim = (__int32) radio_climate; */ + propv.klim=(long)radio_climate; + prop.ens=eno; + prop.kwx=0; + ivar=(long)ModVar; + ipol=(long)pol; + qlrps(frq_mhz, 0.0, eno, ipol, eps, sgm, prop); + qlra(kst, propv.klim, ivar, prop, propv); + + if (propv.lvar<1) + propv.lvar=1; + + lrprop(dist_km*1000.0, prop, propa); + fs=32.45+20.0*log10(frq_mhz)+20.0*log10(prop.dist/1000.0); + xlb=fs+avar(zt, zl, zc, prop, propv); + dbloss=xlb; + + if (prop.kwx==0) + errnum=0; + else + errnum=prop.kwx; +} + +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) +{ + char strmode[200]; + int errnum; + double dbloss; + + 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); + + return dbloss; +} + +double ITMVersion() +{ + return 7.0; +} + diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..79b346f --- /dev/null +++ b/main.cpp @@ -0,0 +1,4287 @@ +/****************************************************************************\ +* Signal Server 1.3.3: Server optimised SPLAT! by Alex Farrant * +****************************************************************************** +* SPLAT! Project started in 1997 by John A. Magliacane, KD2BD * +* * +****************************************************************************** +* Please consult the SPLAT! documentation for a complete list of * +* individuals who have contributed to this project. * +****************************************************************************** +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU General Public License as published by the * +* Free Software Foundation; either version 2 of the License or any later * +* version. * +* * +* This program is distributed in the hope that it will useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License * +* for more details. * +* * +****************************************************************************** +* g++ -Wall -O3 -s -lm -fomit-frame-pointer itm.cpp main.cpp -o ss * +\****************************************************************************/ + +#include +#include +#include +#include +#include +#include + +#define GAMMA 2.5 +#define MAXPAGES 64 +#define ARRAYSIZE 76810 +#define IPPD 1200 + +#ifndef PI +#define PI 3.141592653589793 +#endif + +#ifndef TWOPI +#define TWOPI 6.283185307179586 +#endif + +#ifndef HALFPI +#define HALFPI 1.570796326794896 +#endif + +#define DEG2RAD 1.74532925199e-02 +#define EARTHRADIUS 20902230.97 +#define METERS_PER_MILE 1609.344 +#define METERS_PER_FOOT 0.3048 +#define KM_PER_MILE 1.609344 +#define FOUR_THIRDS 1.3333333333333 + +char string[255], sdf_path[255], udt_file[255], opened=0, gpsav=0, ss_name[16], + ss_version[6], dashes[80]; + +double earthradius, max_range=0.0, forced_erp, dpp, ppd, + fzone_clearance=0.6, forced_freq, clutter, lat,lon,txh,tercon,terdic, + north,east,south,west; + +int min_north=90, max_north=-90, min_west=360, max_west=-1, ippd, mpi, + max_elevation=-32768, min_elevation=32768, bzerror, contour_threshold, + pred,pblue,pgreen,ter,multiplier=256,debug=0,loops=64,jgets=0, MAXRAD; + +unsigned char got_elevation_pattern, got_azimuth_pattern, metric=0, dbm=0; + + + +struct site { double lat; + double lon; + float alt; + char name[50]; + char filename[255]; + } site; + +struct path { double lat[ARRAYSIZE]; + double lon[ARRAYSIZE]; + double elevation[ARRAYSIZE]; + double distance[ARRAYSIZE]; + int length; + } path; + +struct dem { int min_north; + int max_north; + int min_west; + int max_west; + int max_el; + int min_el; + short data[IPPD][IPPD]; + unsigned char mask[IPPD][IPPD]; + unsigned char signal[IPPD][IPPD]; + } dem[MAXPAGES]; + +struct LR { double eps_dielect; + double sgm_conductivity; + double eno_ns_surfref; + double frq_mhz; + double conf; + double rel; + double erp; + int radio_climate; + int pol; + float antenna_pattern[361][1001]; + } LR; + +struct region { unsigned char color[128][3]; + int level[128]; + int levels; + } region; + +double elev[ARRAYSIZE+10]; + + +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); + +double arccos(double x, double y) +{ + /* This function implements the arc cosine function, + returning a value between 0 and TWOPI. */ + + double result=0.0; + + if (y>0.0) + result=acos(x/y); + + if (y<0.0) + result=PI+acos(x/y); + + return result; +} + + +int ReduceAngle(double angle) +{ + /* This function normalizes the argument to + an integer angle between 0 and 180 degrees */ + + double temp; + + temp=acos(cos(angle*DEG2RAD)); + + return (int)rint(temp/DEG2RAD); +} + +double LonDiff(double lon1, double lon2) +{ + /* This function returns the short path longitudinal + difference between longitude1 and longitude2 + as an angle between -180.0 and +180.0 degrees. + If lon1 is west of lon2, the result is positive. + If lon1 is east of lon2, the result is negative. */ + + double diff; + + diff=lon1-lon2; + + if (diff<=-180.0) + diff+=360.0; + + if (diff>=180.0) + diff-=360.0; + + return diff; +} + +char *dec2dms(double decimal) +{ + /* Converts decimal degrees to degrees, minutes, seconds, + (DMS) and returns the result as a character string. */ + + char sign; + int degrees, minutes, seconds; + double a, b, c, d; + + if (decimal<0.0) + { + decimal=-decimal; + sign=-1; + } + + else + sign=1; + + a=floor(decimal); + b=60.0*(decimal-a); + c=floor(b); + d=60.0*(b-c); + + degrees=(int)a; + minutes=(int)c; + seconds=(int)d; + + if (seconds<0) + seconds=0; + + if (seconds>59) + seconds=59; + + string[0]=0; + snprintf(string,250,"%d%c %d\' %d\"", degrees*sign, 176, minutes, seconds); + return (string); +} + +int PutMask(double lat, double lon, int value) +{ + /* Lines, text, markings, and coverage areas are stored in a + mask that is combined with topology data when topographic + maps are generated by ss. This function sets and resets + bits in the mask based on the latitude and longitude of the + area pointed to. */ + + int x, y, indx; + char found; + + for (indx=0, found=0; indx=0 && x<=mpi && y>=0 && y<=mpi) + found=1; + else + indx++; + } + + if (found) + { + dem[indx].mask[x][y]=value; + return ((int)dem[indx].mask[x][y]); + } + + else + return -1; +} + +int OrMask(double lat, double lon, int value) +{ + /* Lines, text, markings, and coverage areas are stored in a + mask that is combined with topology data when topographic + maps are generated by ss. This function sets bits in + the mask based on the latitude and longitude of the area + pointed to. */ + + int x, y, indx; + char found; + + for (indx=0, found=0; indx=0 && x<=mpi && y>=0 && y<=mpi) + found=1; + else + indx++; + } + + if (found) + { + dem[indx].mask[x][y]|=value; + return ((int)dem[indx].mask[x][y]); + } + + else + return -1; +} + +int GetMask(double lat, double lon) +{ + /* This function returns the mask bits based on the latitude + and longitude given. */ + + return (OrMask(lat,lon,0)); +} + +int PutSignal(double lat, double lon, unsigned char signal) +{ + /* This function writes a signal level (0-255) + at the specified location for later recall. */ + + int x, y, indx; + char found; + + for (indx=0, found=0; indx=0 && x<=mpi && y>=0 && y<=mpi) + found=1; + else + indx++; + } + + if (found) + { + dem[indx].signal[x][y]=signal; + return (dem[indx].signal[x][y]); + } + + else + return 0; +} + +unsigned char GetSignal(double lat, double lon) +{ + /* This function reads the signal level (0-255) at the + specified location that was previously written by the + complimentary PutSignal() function. */ + + int x, y, indx; + char found; + + for (indx=0, found=0; indx=0 && x<=mpi && y>=0 && y<=mpi) + found=1; + else + indx++; + } + + if (found) + return (dem[indx].signal[x][y]); + else + return 0; +} + +double GetElevation(struct site location) +{ + /* This function returns the elevation (in feet) of any location + represented by the digital elevation model data in memory. + Function returns -5000.0 for locations not found in memory. */ + + char found; + int x, y, indx; + double elevation; + + for (indx=0, found=0; indx=0 && x<=mpi && y>=0 && y<=mpi) + found=1; + else + indx++; + } + + if (found) + elevation=3.28084*dem[indx].data[x][y]; + else + elevation=-5000.0; + + return elevation; +} + +int AddElevation(double lat, double lon, double height) +{ + /* This function adds a user-defined terrain feature + (in meters AGL) to the digital elevation model data + in memory. Does nothing and returns 0 for locations + not found in memory. */ + + char found; + int x, y, indx; + + + for (indx=0, found=0; indx=0 && x<=mpi && y>=0 && y<=mpi) + found=1; + else + indx++; + } + + if (found) + dem[indx].data[x][y]+=(short)rint(height); + + return found; +} + +double Distance(struct site site1, struct site site2) +{ + /* This function returns the great circle distance + in miles between any two site locations. */ + + double lat1, lon1, lat2, lon2, distance; + + lat1=site1.lat*DEG2RAD; + lon1=site1.lon*DEG2RAD; + lat2=site2.lat*DEG2RAD; + lon2=site2.lon*DEG2RAD; + + distance=3959.0*acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos((lon1)-(lon2))); + + return distance; +} + +double Azimuth(struct site source, struct site destination) +{ + /* This function returns the azimuth (in degrees) to the + destination as seen from the location of the source. */ + + double dest_lat, dest_lon, src_lat, src_lon, + beta, azimuth, diff, num, den, fraction; + + dest_lat=destination.lat*DEG2RAD; + dest_lon=destination.lon*DEG2RAD; + + src_lat=source.lat*DEG2RAD; + src_lon=source.lon*DEG2RAD; + + /* Calculate Surface Distance */ + + beta=acos(sin(src_lat)*sin(dest_lat)+cos(src_lat)*cos(dest_lat)*cos(src_lon-dest_lon)); + + /* Calculate Azimuth */ + + num=sin(dest_lat)-(sin(src_lat)*cos(beta)); + den=cos(src_lat)*sin(beta); + fraction=num/den; + + /* Trap potential problems in acos() due to rounding */ + + if (fraction>=1.0) + fraction=1.0; + + if (fraction<=-1.0) + fraction=-1.0; + + /* Calculate azimuth */ + + azimuth=acos(fraction); + + /* Reference it to True North */ + + diff=dest_lon-src_lon; + + if (diff<=-PI) + diff+=TWOPI; + + if (diff>=PI) + diff-=TWOPI; + + if (diff>0.0) + azimuth=TWOPI-azimuth; + + return (azimuth/DEG2RAD); +} + +double ElevationAngle(struct site source, struct site destination) +{ + /* This function returns the angle of elevation (in degrees) + of the destination as seen from the source location. + A positive result represents an angle of elevation (uptilt), + while a negative result represents an angle of depression + (downtilt), as referenced to a normal to the center of + the earth. */ + + register double a, b, dx; + + a=GetElevation(destination)+destination.alt+earthradius; + b=GetElevation(source)+source.alt+earthradius; + + dx=5280.0*Distance(source,destination); + + /* Apply the Law of Cosines */ + + return ((180.0*(acos(((b*b)+(dx*dx)-(a*a))/(2.0*b*dx)))/PI)-90.0); +} + +void ReadPath(struct site source, struct site destination) +{ + /* This function generates a sequence of latitude and + longitude positions between source and destination + locations along a great circle path, and stores + elevation and distance information for points + along that path in the "path" structure. */ + + int c; + double azimuth, distance, lat1, lon1, beta, den, num, + lat2, lon2, total_distance, dx, dy, path_length, + miles_per_sample, samples_per_radian=68755.0; + struct site tempsite; + + lat1=source.lat*DEG2RAD; + lon1=source.lon*DEG2RAD; + + lat2=destination.lat*DEG2RAD; + lon2=destination.lon*DEG2RAD; + + samples_per_radian=ppd*57.295833; + + azimuth=Azimuth(source,destination)*DEG2RAD; + + total_distance=Distance(source,destination); + + if (total_distance>(30.0/ppd)) + { + dx=samples_per_radian*acos(cos(lon1-lon2)); + dy=samples_per_radian*acos(cos(lat1-lat2)); + + path_length=sqrt((dx*dx)+(dy*dy)); + + miles_per_sample=total_distance/path_length; + } + + else + { + c=0; + dx=0.0; + dy=0.0; + path_length=0.0; + miles_per_sample=0.0; + total_distance=0.0; + + lat1=lat1/DEG2RAD; + lon1=lon1/DEG2RAD; + + path.lat[c]=lat1; + path.lon[c]=lon1; + path.elevation[c]=GetElevation(source); + path.distance[c]=0.0; + } + + for (distance=0.0, c=0; (total_distance!=0.0 && distance<=total_distance && cHALFPI-lat1)) + lon2=lon1+PI; + + else if (azimuth==HALFPI && (beta>HALFPI+lat1)) + lon2=lon1+PI; + + else if (fabs(num/den)>1.0) + lon2=lon1; + + else + { + if ((PI-azimuth)>=0.0) + lon2=lon1-arccos(num,den); + else + lon2=lon1+arccos(num,den); + } + + while (lon2<0.0) + lon2+=TWOPI; + + while (lon2>TWOPI) + lon2-=TWOPI; + + lat2=lat2/DEG2RAD; + lon2=lon2/DEG2RAD; + + path.lat[c]=lat2; + path.lon[c]=lon2; + tempsite.lat=lat2; + tempsite.lon=lon2; + path.elevation[c]=GetElevation(tempsite); + path.distance[c]=distance; + } + + /* Make sure exact destination point is recorded at path.length-1 */ + + if (c=cos_test_angle) + { + block=1; + first_obstruction_angle=((acos(cos_test_angle))/DEG2RAD)-90.0; + } + } + + if (block) + elevation=first_obstruction_angle; + + else + elevation=((acos(cos_xmtr_angle))/DEG2RAD)-90.0; + + path=temp; + + return elevation; +} + +double AverageTerrain(struct site source, double azimuthx, double start_distance, double end_distance) +{ + /* This function returns the average terrain calculated in + the direction of "azimuth" (degrees) between "start_distance" + and "end_distance" (miles) from the source location. If + the terrain is all water (non-critical error), -5000.0 is + returned. If not enough SDF data has been loaded into + memory to complete the survey (critical error), then + -9999.0 is returned. */ + + int c, samples, endpoint; + double beta, lat1, lon1, lat2, lon2, num, den, azimuth, terrain=0.0; + struct site destination; + + lat1=source.lat*DEG2RAD; + lon1=source.lon*DEG2RAD; + + /* Generate a path of elevations between the source + location and the remote location provided. */ + + beta=end_distance/3959.0; + + azimuth=DEG2RAD*azimuthx; + + lat2=asin(sin(lat1)*cos(beta)+cos(azimuth)*sin(beta)*cos(lat1)); + num=cos(beta)-(sin(lat1)*sin(lat2)); + den=cos(lat1)*cos(lat2); + + if (azimuth==0.0 && (beta>HALFPI-lat1)) + lon2=lon1+PI; + + else if (azimuth==HALFPI && (beta>HALFPI+lat1)) + lon2=lon1+PI; + + else if (fabs(num/den)>1.0) + lon2=lon1; + + else + { + if ((PI-azimuth)>=0.0) + lon2=lon1-arccos(num,den); + else + lon2=lon1+arccos(num,den); + } + + while (lon2<0.0) + lon2+=TWOPI; + + while (lon2>TWOPI) + lon2-=TWOPI; + + lat2=lat2/DEG2RAD; + lon2=lon2/DEG2RAD; + + destination.lat=lat2; + destination.lon=lon2; + + /* If SDF data is missing for the endpoint of + the radial, then the average terrain cannot + be accurately calculated. Return -9999.0 */ + + if (GetElevation(destination)<-4999.0) + return (-9999.0); + else + { + ReadPath(source,destination); + + endpoint=path.length; + + /* Shrink the length of the radial if the + outermost portion is not over U.S. land. */ + + for (c=endpoint-1; c>=0 && path.elevation[c]==0.0; c--); + + endpoint=c+1; + + for (c=0, samples=0; c=start_distance) + { + terrain+=(path.elevation[c]==0.0?path.elevation[c]:path.elevation[c]+clutter); + samples++; + } + } + + if (samples==0) + terrain=-5000.0; /* No land */ + else + terrain=(terrain/(double)samples); + + return terrain; + } +} + +double haat(struct site antenna) +{ + /* This function returns the antenna's Height Above Average + Terrain (HAAT) based on FCC Part 73.313(d). If a critical + error occurs, such as a lack of SDF data to complete the + survey, -5000.0 is returned. */ + + int azi, c; + char error=0; + double terrain, avg_terrain, haat, sum=0.0; + + /* Calculate the average terrain between 2 and 10 miles + from the antenna site at azimuths of 0, 45, 90, 135, + 180, 225, 270, and 315 degrees. */ + + for (c=0, azi=0; azi<=315 && error==0; azi+=45) + { + terrain=AverageTerrain(antenna, (double)azi, 2.0, 10.0); + + if (terrain<-9998.0) /* SDF data is missing */ + error=1; + + if (terrain>-4999.0) /* It's land, not water */ + { + sum+=terrain; /* Sum of averages */ + c++; + } + } + + if (error) + return -5000.0; + else + { + avg_terrain=(sum/(double)c); + haat=(antenna.alt+GetElevation(antenna))-avg_terrain; + return haat; + } +} +double ReadBearing(char *input) +{ + /* This function takes numeric input in the form of a character + string, and returns an equivalent bearing in degrees as a + decimal number (double). The input may either be expressed + in decimal format (40.139722) or degree, minute, second + format (40 08 23). This function also safely handles + extra spaces found either leading, trailing, or + embedded within the numbers expressed in the + input string. Decimal seconds are permitted. */ + + double seconds, bearing=0.0; + char string[20]; + int a, b, length, degrees, minutes; + + /* Copy "input" to "string", and ignore any extra + spaces that might be present in the process. */ + + string[0]=0; + length=strlen(input); + + for (a=0, b=0; a360.0 || bearing<-360.0) + bearing=0.0; + + return bearing; +} + +void LoadPAT(char *filename) +{ + /* This function reads and processes antenna pattern (.az + and .el) files that correspond in name to previously + loaded ss .lrp files. */ + + int a, b, w, x, y, z, last_index, next_index, span; + char string[255], azfile[255], elfile[255], *pointer=NULL, *s=NULL; + float az, xx, elevation, amplitude, rotation, valid1, valid2, + delta, azimuth[361], azimuth_pattern[361], el_pattern[10001], + elevation_pattern[361][1001], slant_angle[361], tilt, + mechanical_tilt=0.0, tilt_azimuth, tilt_increment, sum; + FILE *fd=NULL; + unsigned char read_count[10001]; + + for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++) + { + azfile[x]=filename[x]; + elfile[x]=filename[x]; + } + + azfile[x]='.'; + azfile[x+1]='a'; + azfile[x+2]='z'; + azfile[x+3]=0; + + elfile[x]='.'; + elfile[x+1]='e'; + elfile[x+2]='l'; + elfile[x+3]=0; + + rotation=0.0; + + got_azimuth_pattern=0; + got_elevation_pattern=0; + + /* Load .az antenna pattern file */ + + fd=fopen(azfile,"r"); + + if (fd!=NULL) + { + /* Clear azimuth pattern array */ + + for (x=0; x<=360; x++) + { + azimuth[x]=0.0; + read_count[x]=0; + } + + + /* Read azimuth pattern rotation + in degrees measured clockwise + from true North. */ + + s=fgets(string,254,fd); + pointer=strchr(string,';'); + + if (pointer!=NULL) + *pointer=0; + + sscanf(string,"%f",&rotation); + + + /* Read azimuth (degrees) and corresponding + normalized field radiation pattern amplitude + (0.0 to 1.0) until EOF is reached. */ + + s=fgets(string,254,fd); + pointer=strchr(string,';'); + + if (pointer!=NULL) + *pointer=0; + + sscanf(string,"%f %f",&az, &litude); + + do + { + x=(int)rintf(az); + + if (x>=0 && x<=360 && fd!=NULL) + { + azimuth[x]+=amplitude; + read_count[x]++; + } + + s=fgets(string,254,fd); + pointer=strchr(string,';'); + + if (pointer!=NULL) + *pointer=0; + + sscanf(string,"%f %f",&az, &litude); + + } while (feof(fd)==0); + + fclose(fd); + + + /* Handle 0=360 degree ambiguity */ + + if ((read_count[0]==0) && (read_count[360]!=0)) + { + read_count[0]=read_count[360]; + azimuth[0]=azimuth[360]; + } + + if ((read_count[0]!=0) && (read_count[360]==0)) + { + read_count[360]=read_count[0]; + azimuth[360]=azimuth[0]; + } + + /* Average pattern values in case more than + one was read for each degree of azimuth. */ + + for (x=0; x<=360; x++) + { + if (read_count[x]>1) + azimuth[x]/=(float)read_count[x]; + } + + /* Interpolate missing azimuths + to completely fill the array */ + + last_index=-1; + next_index=-1; + + for (x=0; x<=360; x++) + { + if (read_count[x]!=0) + { + if (last_index==-1) + last_index=x; + else + next_index=x; + } + + if (last_index!=-1 && next_index!=-1) + { + valid1=azimuth[last_index]; + valid2=azimuth[next_index]; + + span=next_index-last_index; + delta=(valid2-valid1)/(float)span; + + for (y=last_index+1; y=360) + y-=360; + + azimuth_pattern[y]=azimuth[x]; + } + + azimuth_pattern[360]=azimuth_pattern[0]; + + got_azimuth_pattern=255; + } + + /* Read and process .el file */ + + fd=fopen(elfile,"r"); + + if (fd!=NULL) + { + for (x=0; x<=10000; x++) + { + el_pattern[x]=0.0; + read_count[x]=0; + } + + /* Read mechanical tilt (degrees) and + tilt azimuth in degrees measured + clockwise from true North. */ + + s=fgets(string,254,fd); + pointer=strchr(string,';'); + + if (pointer!=NULL) + *pointer=0; + + sscanf(string,"%f %f",&mechanical_tilt, &tilt_azimuth); + + /* Read elevation (degrees) and corresponding + normalized field radiation pattern amplitude + (0.0 to 1.0) until EOF is reached. */ + + s=fgets(string,254,fd); + pointer=strchr(string,';'); + + if (pointer!=NULL) + *pointer=0; + + sscanf(string,"%f %f", &elevation, &litude); + + while (feof(fd)==0) + { + /* Read in normalized radiated field values + for every 0.01 degrees of elevation between + -10.0 and +90.0 degrees */ + + x=(int)rintf(100.0*(elevation+10.0)); + + if (x>=0 && x<=10000) + { + el_pattern[x]+=amplitude; + read_count[x]++; + } + + s=fgets(string,254,fd); + pointer=strchr(string,';'); + + if (pointer!=NULL) + *pointer=0; + + sscanf(string,"%f %f", &elevation, &litude); + } + + fclose(fd); + + /* Average the field values in case more than + one was read for each 0.01 degrees of elevation. */ + + for (x=0; x<=10000; x++) + { + if (read_count[x]>1) + el_pattern[x]/=(float)read_count[x]; + } + + /* Interpolate between missing elevations (if + any) to completely fill the array and provide + radiated field values for every 0.01 degrees of + elevation. */ + + last_index=-1; + next_index=-1; + + for (x=0; x<=10000; x++) + { + if (read_count[x]!=0) + { + if (last_index==-1) + last_index=x; + else + next_index=x; + } + + if (last_index!=-1 && next_index!=-1) + { + valid1=el_pattern[last_index]; + valid2=el_pattern[next_index]; + + span=next_index-last_index; + delta=(valid2-valid1)/(float)span; + + for (y=last_index+1; y=360) + y-=360; + + while (y<0) + y+=360; + + if (x<=180) + slant_angle[y]=-(tilt_increment*(90.0-xx)); + + if (x>180) + slant_angle[y]=-(tilt_increment*(xx-270.0)); + } + } + + slant_angle[360]=slant_angle[0]; /* 360 degree wrap-around */ + + for (w=0; w<=360; w++) + { + tilt=slant_angle[w]; + + /** Convert tilt angle to + an array index offset **/ + + y=(int)rintf(100.0*tilt); + + /* Copy shifted el_pattern[10001] field + values into elevation_pattern[361][1001] + at the corresponding azimuth, downsampling + (averaging) along the way in chunks of 10. */ + + for (x=y, z=0; z<=1000; x+=10, z++) + { + for (sum=0.0, a=0; a<10; a++) + { + b=a+x; + + if (b>=0 && b<=10000) + sum+=el_pattern[b]; + if (b<0) + sum+=el_pattern[0]; + if (b>10000) + sum+=el_pattern[10000]; + } + + elevation_pattern[w][z]=sum/10.0; + } + } + + got_elevation_pattern=255; + } + + for (x=0; x<=360; x++) + { + for (y=0; y<=1000; y++) + { + if (got_elevation_pattern) + elevation=elevation_pattern[x][y]; + else + elevation=1.0; + + if (got_azimuth_pattern) + az=azimuth_pattern[x]; + else + az=1.0; + + LR.antenna_pattern[x][y]=az*elevation; + } + } +} + +int LoadSDF_SDF(char *name) +{ + /* This function reads uncompressed ss Data Files (.sdf) + containing digital elevation model data into memory. + Elevation data, maximum and minimum elevations, and + quadrangle limits are stored in the first available + dem[] structure. */ + + int x, y, data, indx, minlat, minlon, maxlat, maxlon,j; + char found, free_page=0, line[20], jline[20], sdf_file[255], + path_plus_name[255], *s=NULL,*junk=NULL; + + + FILE *fd; + + for (x=0; name[x]!='.' && name[x]!=0 && x<250; x++) + sdf_file[x]=name[x]; + + sdf_file[x]=0; + + /* Parse filename for minimum latitude and longitude values */ + + sscanf(sdf_file,"%d:%d:%d:%d",&minlat,&maxlat,&minlon,&maxlon); + + sdf_file[x]='.'; + sdf_file[x+1]='s'; + sdf_file[x+2]='d'; + sdf_file[x+3]='f'; + sdf_file[x+4]=0; + + /* Is it already in memory? */ + + + for (indx=0, found=0; indx=0 && indxdem[indx].max_el) + dem[indx].max_el=data; + + if (datamax_elevation) + max_elevation=dem[indx].max_el; + + if (max_north==-90) + max_north=dem[indx].max_north; + + else if (dem[indx].max_north>max_north) + max_north=dem[indx].max_north; + + if (min_north==90) + min_north=dem[indx].min_north; + + else if (dem[indx].min_northmax_west) + max_west=dem[indx].max_west; + } + + else + { + if (dem[indx].max_westmin_west) + min_west=dem[indx].min_west; + } + } + + if (debug==1){ + fprintf(stdout," Done!\n"); + fflush(stdout); + } + + return 1; + } + + else + return -1; + } + + else + return 0; +} +char LoadSDF(char *name) +{ + /* This function loads the requested SDF file from the filesystem. + It first tries to invoke the LoadSDF_SDF() function to load an + uncompressed SDF file (since uncompressed files load slightly + faster). If that attempt fails, then it tries to load a + compressed SDF file by invoking the LoadSDF_BZ() function. + If that fails, then we can assume that no elevation data + exists for the region requested, and that the region + requested must be entirely over water. */ + + int x, y, indx, minlat, minlon, maxlat, maxlon; + char found, free_page=0; + int return_value=-1; + + /* Try to load an uncompressed SDF first. */ + + return_value=LoadSDF_SDF(name); + + /* If that fails, try loading a compressed SDF. */ + + //if (return_value==0 || return_value==-1) + // return_value=LoadSDF_BZ(name); + + /* If neither format can be found, then assume the area is water. */ + + if (return_value==0 || return_value==-1) + { + + /* Parse SDF name for minimum latitude and longitude values */ + + sscanf(name,"%d:%d:%d:%d",&minlat,&maxlat,&minlon,&maxlon); + + /* Is it already in memory? */ + + for (indx=0, found=0; indx=0 && indx0) + dem[indx].min_el=0; + } + + if (dem[indx].min_elmax_elevation) + max_elevation=dem[indx].max_el; + + if (max_north==-90) + max_north=dem[indx].max_north; + + else if (dem[indx].max_north>max_north) + max_north=dem[indx].max_north; + + if (min_north==90) + min_north=dem[indx].min_north; + + else if (dem[indx].min_northmax_west) + max_west=dem[indx].max_west; + } + + else + { + if (dem[indx].max_westmin_west) + min_west=dem[indx].min_west; + } + } + if(debug==1){ + fprintf(stdout," Done!\n"); + fflush(stdout); + } + return_value=1; + } + } + + return return_value; +} + +void PlotPath(struct site source, struct site destination, char mask_value) +{ + /* This function analyzes the path between the source and + destination locations. It determines which points along + the path have line-of-sight visibility to the source. + Points along with path having line-of-sight visibility + to the source at an AGL altitude equal to that of the + destination location are stored by setting bit 1 in the + mask[][] array, which are displayed in green when PPM + maps are later generated by ss. */ + + char block; + int x, y; + register double cos_xmtr_angle, cos_test_angle, test_alt; + double distance, rx_alt, tx_alt; + + ReadPath(source,destination); + + for (y=0; y=0 && block==0; x--) + { + distance=5280.0*(path.distance[y]-path.distance[x]); + test_alt=earthradius+(path.elevation[x]==0.0?path.elevation[x]:path.elevation[x]+clutter); + + cos_test_angle=((rx_alt*rx_alt)+(distance*distance)-(test_alt*test_alt))/(2.0*rx_alt*distance); + + /* Compare these two angles to determine if + an obstruction exists. Since we're comparing + the cosines of these angles rather than + the angles themselves, the following "if" + statement is reversed from what it would + be if the actual angles were compared. */ + + if (cos_xmtr_angle>=cos_test_angle) + block=1; + } + + if (block==0) + OrMask(path.lat[y],path.lon[y],mask_value); + } + } +} + +void PlotLRPath(struct site source, struct site destination, unsigned char mask_value, FILE *fd) +{ + /* This function plots the RF path loss between source and + destination points based on the Longley-Rice propagation + model, taking into account antenna pattern data, if available. */ + + int x, y, ifs, ofs, errnum; + char block=0, strmode[100]; + double loss, azimuth, pattern=0.0, + xmtr_alt, dest_alt, xmtr_alt2, dest_alt2, + cos_rcvr_angle, cos_test_angle=0.0, test_alt, + elevation=0.0, distance=0.0, four_thirds_earth, + field_strength=0.0, rxp, dBm; + struct site temp; + + ReadPath(source,destination); + + four_thirds_earth=FOUR_THIRDS*EARTHRADIUS; + + /* Copy elevations plus clutter along path into the elev[] array. */ + + for (x=1; x1.0) + cos_rcvr_angle=1.0; + + if (cos_rcvr_angle<-1.0) + cos_rcvr_angle=-1.0; + + if (got_elevation_pattern || fd!=NULL) + { + /* Determine the elevation angle to the first obstruction + along the path IF elevation pattern data is available + or an output (.ano) file has been designated. */ + + for (x=2, block=0; (x1.0) + cos_test_angle=1.0; + + if (cos_test_angle<-1.0) + cos_test_angle=-1.0; + + /* Compare these two angles to determine if + an obstruction exists. Since we're comparing + the cosines of these angles rather than + the angles themselves, the sense of the + following "if" statement is reversed from + what it would be if the angles themselves + were compared. */ + + if (cos_rcvr_angle>=cos_test_angle) + block=1; + } + + if (block) + elevation=((acos(cos_test_angle))/DEG2RAD)-90.0; + else + elevation=((acos(cos_rcvr_angle))/DEG2RAD)-90.0; + } + + /* Determine attenuation for each point along the + path using Longley-Rice's point_to_point mode + starting at y=2 (number_of_points = 1), the + shortest distance terrain can play a role in + path loss. */ + + elev[0]=y-1; /* (number of points - 1) */ + + /* Distance between elevation samples */ + + elev[1]=METERS_PER_MILE*(path.distance[y]-path.distance[y-1]); + + point_to_point(elev,source.alt*METERS_PER_FOOT, + destination.alt*METERS_PER_FOOT, LR.eps_dielect, + LR.sgm_conductivity, LR.eno_ns_surfref, LR.frq_mhz, + LR.radio_climate, LR.pol, LR.conf, LR.rel, loss, + strmode, errnum); + + temp.lat=path.lat[y]; + temp.lon=path.lon[y]; + + azimuth=(Azimuth(source,temp)); + + if (fd!=NULL) + fprintf(fd,"%.7f, %.7f, %.3f, %.3f, ",path.lat[y], path.lon[y], azimuth, elevation); + + /* If ERP==0, write path loss to alphanumeric + output file. Otherwise, write field strength + or received power level (below), as appropriate. */ + + if (fd!=NULL && LR.erp==0.0) + fprintf(fd,"%.2f",loss); + + /* Integrate the antenna's radiation + pattern into the overall path loss. */ + + x=(int)rint(10.0*(10.0-elevation)); + + if (x>=0 && x<=1000) + { + azimuth=rint(azimuth); + + pattern=(double)LR.antenna_pattern[(int)azimuth][x]; + + if (pattern!=0.0) + { + pattern=20.0*log10(pattern); + loss-=pattern; + } + } + + if (LR.erp!=0.0) + { + if (dbm) + { + /* dBm is based on EIRP (ERP + 2.14) */ + + rxp=LR.erp/(pow(10.0,(loss-2.14)/10.0)); + + dBm=10.0*(log10(rxp*1000.0)); + + if (fd!=NULL) + fprintf(fd,"%.3f",dBm); + + /* Scale roughly between 0 and 255 */ + + ifs=200+(int)rint(dBm); + + if (ifs<0) + ifs=0; + + if (ifs>255) + ifs=255; + + ofs=GetSignal(path.lat[y],path.lon[y]); + + if (ofs>ifs) + ifs=ofs; + + PutSignal(path.lat[y],path.lon[y],(unsigned char)ifs); + } + + else + { + field_strength=(139.4+(20.0*log10(LR.frq_mhz))-loss)+(10.0*log10(LR.erp/1000.0)); + + ifs=100+(int)rint(field_strength); + + if (ifs<0) + ifs=0; + + if (ifs>255) + ifs=255; + + ofs=GetSignal(path.lat[y],path.lon[y]); + + if (ofs>ifs) + ifs=ofs; + + PutSignal(path.lat[y],path.lon[y],(unsigned char)ifs); + + if (fd!=NULL) + fprintf(fd,"%.3f",field_strength); + } + } + + else + { + if (loss>255) + ifs=255; + else + ifs=(int)rint(loss); + + ofs=GetSignal(path.lat[y],path.lon[y]); + + if (ofs0.0 && debug) + fprintf(stdout," and %.2f %s of ground clutter",metric?clutter*METERS_PER_FOOT:clutter,metric?"meters":"feet"); + + fprintf(stdout,"...\n\n 0%c to 25%c ",37,37); + fflush(stdout); + + /* th=pixels/degree divided by 64 loops per + progress indicator symbol (.oOo) printed. */ + + th=ppd/loops; + + z=(int)(th*ReduceAngle(max_west-min_west)); + + minwest=dpp+(double)min_west; + maxnorth=(double)max_north-dpp; + + for (lon=minwest, x=0, y=0; (LonDiff(lon,(double)max_west)<=0.0); y++, lon=minwest+(dpp*(double)y)) + { + if (lon>=360.0) + lon-=360.0; + + edge.lat=max_north; + edge.lon=lon; + edge.alt=altitude; + + PlotPath(source,edge,mask_value); + count++; + + if (count==z) + { + fprintf(stdout,"%c",symbol[x]); + fflush(stdout); + count=0; + + if (x==3) + x=0; + else + x++; + } + } + + count=0; + fprintf(stdout,"\n25%c to 50%c ",37,37); + fflush(stdout); + + z=(int)(th*(double)(max_north-min_north)); + + for (lat=maxnorth, x=0, y=0; lat>=(double)min_north; y++, lat=maxnorth-(dpp*(double)y)) + { + edge.lat=lat; + edge.lon=min_west; + edge.alt=altitude; + + PlotPath(source,edge,mask_value); + count++; + + if (count==z) + { + //fprintf(stdout,"%c",symbol[x]); + //fflush(stdout); + count=0; + + if (x==3) + x=0; + else + x++; + } + } + + count=0; + fprintf(stdout,"\n50%c to 75%c ",37,37); + fflush(stdout); + + z=(int)(th*ReduceAngle(max_west-min_west)); + + for (lon=minwest, x=0, y=0; (LonDiff(lon,(double)max_west)<=0.0); y++, lon=minwest+(dpp*(double)y)) + { + if (lon>=360.0) + lon-=360.0; + + edge.lat=min_north; + edge.lon=lon; + edge.alt=altitude; + + PlotPath(source,edge,mask_value); + count++; + + if (count==z) + { + //fprintf(stdout,"%c",symbol[x]); + //fflush(stdout); + count=0; + + if (x==3) + x=0; + else + x++; + } + } + + count=0; + fprintf(stdout,"\n75%c to 100%c ",37,37); + fflush(stdout); + + z=(int)(th*(double)(max_north-min_north)); + + for (lat=(double)min_north, x=0, y=0; lat<(double)max_north; y++, lat=(double)min_north+(dpp*(double)y)) + { + edge.lat=lat; + edge.lon=max_west; + edge.alt=altitude; + + PlotPath(source,edge,mask_value); + count++; + + if (count==z) + { + //fprintf(stdout,"%c",symbol[x]); + //fflush(stdout); + count=0; + + if (x==3) + x=0; + else + x++; + } + } + + fprintf(stdout,"\nDone!\n"); + fflush(stdout); + + /* Assign next mask value */ + + switch (mask_value) + { + case 1: + mask_value=8; + break; + + case 8: + mask_value=16; + break; + + case 16: + mask_value=32; + } +} + +void PlotLRMap(struct site source, double altitude, char *plo_filename) +{ + /* This function performs a 360 degree sweep around the + transmitter site (source location), and plots the + Longley-Rice attenuation on the ss generated + topographic map based on a receiver located at + the specified altitude (in feet AGL). Results + are stored in memory, and written out in the form + of a topographic map when the DoPathLoss() or + DoSigStr() functions are later invoked. */ + + int y, z, count; + struct site edge; + double lat, lon, minwest, maxnorth, th; + unsigned char x, symbol[4]; + static unsigned char mask_value=1; + FILE *fd=NULL; + + minwest=dpp+(double)min_west; + maxnorth=(double)max_north-dpp; + + symbol[0]='.'; + symbol[1]='o'; + symbol[2]='O'; + symbol[3]='o'; + + count=0; + if (debug){ + fprintf(stdout,"\nComputing Longley-Rice "); + } + + if (LR.erp==0.0 && debug) + fprintf(stdout,"path loss"); + else + { + if(debug){ + if (dbm) + fprintf(stdout,"signal power level"); + else + fprintf(stdout,"field strength"); + } + } + if (debug){ + fprintf(stdout," contours of \"%s\"\nout to a radius of %.2f %s with Rx antenna(s) at %.2f %s AGL\n",source.name,metric?max_range*KM_PER_MILE:max_range,metric?"kilometers":"miles",metric?altitude*METERS_PER_FOOT:altitude,metric?"meters":"feet"); + } + + if (clutter>0.0 && debug) + fprintf(stdout,"\nand %.2f %s of ground clutter",metric?clutter*METERS_PER_FOOT:clutter,metric?"meters":"feet"); + if(debug){ + fprintf(stdout,"...\n\n 0%c to 25%c ",37,37); + fflush(stdout); + } + if (plo_filename[0]!=0) + fd=fopen(plo_filename,"wb"); + + if (fd!=NULL) + { + /* Write header information to output file */ + + fprintf(fd,"%d, %d\t; max_west, min_west\n%d, %d\t; max_north, min_north\n",max_west, min_west, max_north, min_north); + } + + /* th=pixels/degree divided by 64 loops per + progress indicator symbol (.oOo) printed. */ + + th=ppd/loops; + + z=(int)(th*ReduceAngle(max_west-min_west)); + + for (lon=minwest, x=0, y=0; (LonDiff(lon,(double)max_west)<=0.0); y++, lon=minwest+(dpp*(double)y)) + { + if (lon>=360.0) + lon-=360.0; + + edge.lat=max_north; + edge.lon=lon; + edge.alt=altitude; + + PlotLRPath(source,edge,mask_value,fd); + count++; + + if (count==z) + { + //fprintf(stdout,"%c",symbol[x]); + //fflush(stdout); + count=0; + + if (x==3) + x=0; + else + x++; + } + } + + count=0; + if(debug){ + fprintf(stdout,"\n25%c to 50%c ",37,37); + fflush(stdout); + } + z=(int)(th*(double)(max_north-min_north)); + + for (lat=maxnorth, x=0, y=0; lat>=(double)min_north; y++, lat=maxnorth-(dpp*(double)y)) + { + edge.lat=lat; + edge.lon=min_west; + edge.alt=altitude; + + PlotLRPath(source,edge,mask_value,fd); + count++; + + if (count==z) + { + //fprintf(stdout,"%c",symbol[x]); + //fflush(stdout); + count=0; + + if (x==3) + x=0; + else + x++; + } + } + + count=0; + if(debug){ + fprintf(stdout,"\n50%c to 75%c ",37,37); + fflush(stdout); + } + z=(int)(th*ReduceAngle(max_west-min_west)); + + for (lon=minwest, x=0, y=0; (LonDiff(lon,(double)max_west)<=0.0); y++, lon=minwest+(dpp*(double)y)) + { + if (lon>=360.0) + lon-=360.0; + + edge.lat=min_north; + edge.lon=lon; + edge.alt=altitude; + + PlotLRPath(source,edge,mask_value,fd); + count++; + if (count==z) + { + //fprintf(stdout,"%c",symbol[x]); + //fflush(stdout); + count=0; + + if (x==3) + x=0; + else + x++; + } + + } + + count=0; + if(debug){ + fprintf(stdout,"\n75%c to 100%c ",37,37); + fflush(stdout); + } + z=(int)(th*(double)(max_north-min_north)); + + for (lat=(double)min_north, x=0, y=0; lat<(double)max_north; y++, lat=(double)min_north+(dpp*(double)y)) + { + edge.lat=lat; + edge.lon=max_west; + edge.alt=altitude; + + PlotLRPath(source,edge,mask_value,fd); + count++; + + if (count==z) + { + //fprintf(stdout,"%c",symbol[x]); + //fflush(stdout); + count=0; + + if (x==3) + x=0; + else + x++; + } + } + + if (fd!=NULL) + fclose(fd); + + + if (mask_value<30) + mask_value++; +} + +void LoadSignalColors(struct site xmtr) +{ + int x, y, ok, val[4]; + char filename[255], string[80], *pointer=NULL, *s=NULL; + FILE *fd=NULL; + + for (x=0; xmtr.filename[x]!='.' && xmtr.filename[x]!=0 && x<250; x++) + filename[x]=xmtr.filename[x]; + + filename[x]='.'; + filename[x+1]='s'; + filename[x+2]='c'; + filename[x+3]='f'; + filename[x+4]=0; + + /* Default values */ + + region.level[0]=128; + region.color[0][0]=255; + region.color[0][1]=0; + region.color[0][2]=0; + + region.level[1]=118; + region.color[1][0]=255; + region.color[1][1]=165; + region.color[1][2]=0; + + region.level[2]=108; + region.color[2][0]=255; + region.color[2][1]=206; + region.color[2][2]=0; + + region.level[3]=98; + region.color[3][0]=255; + region.color[3][1]=255; + region.color[3][2]=0; + + region.level[4]=88; + region.color[4][0]=184; + region.color[4][1]=255; + region.color[4][2]=0; + + region.level[5]=78; + region.color[5][0]=0; + region.color[5][1]=255; + region.color[5][2]=0; + + region.level[6]=68; + region.color[6][0]=0; + region.color[6][1]=208; + region.color[6][2]=0; + + region.level[7]=58; + region.color[7][0]=0; + region.color[7][1]=196; + region.color[7][2]=196; + + region.level[8]=48; + region.color[8][0]=0; + region.color[8][1]=148; + region.color[8][2]=255; + + region.level[9]=38; + region.color[9][0]=80; + region.color[9][1]=80; + region.color[9][2]=255; + + region.level[10]=28; + region.color[10][0]=0; + region.color[10][1]=38; + region.color[10][2]=255; + + region.level[11]=18; + region.color[11][0]=142; + region.color[11][1]=63; + region.color[11][2]=255; + + region.level[12]=8; + region.color[12][0]=140; + region.color[12][1]=0; + region.color[12][2]=128; + + region.levels=13; + + fd=fopen(filename,"r"); + + if (fd==NULL) + fd=fopen(filename,"r"); + + if (fd==NULL) + { + fd=fopen(filename,"w"); + + + for (x=0; x255) + val[y]=255; + + if (val[y]<0) + val[y]=0; + } + + region.level[x]=val[0]; + region.color[x][0]=val[1]; + region.color[x][1]=val[2]; + region.color[x][2]=val[3]; + x++; + } + + s=fgets(string,80,fd); + } + + fclose(fd); + region.levels=x; + } +} + +void LoadLossColors(struct site xmtr) +{ + int x, y, ok, val[4]; + char filename[255], string[80], *pointer=NULL, *s=NULL; + FILE *fd=NULL; + + for (x=0; xmtr.filename[x]!='.' && xmtr.filename[x]!=0 && x<250; x++) + filename[x]=xmtr.filename[x]; + + filename[x]='.'; + filename[x+1]='l'; + filename[x+2]='c'; + filename[x+3]='f'; + filename[x+4]=0; + + /* Default values */ + + region.level[0]=80; + region.color[0][0]=255; + region.color[0][1]=0; + region.color[0][2]=0; + + region.level[1]=90; + region.color[1][0]=255; + region.color[1][1]=128; + region.color[1][2]=0; + + region.level[2]=100; + region.color[2][0]=255; + region.color[2][1]=165; + region.color[2][2]=0; + + region.level[3]=110; + region.color[3][0]=255; + region.color[3][1]=206; + region.color[3][2]=0; + + region.level[4]=120; + region.color[4][0]=255; + region.color[4][1]=255; + region.color[4][2]=0; + + region.level[5]=130; + region.color[5][0]=184; + region.color[5][1]=255; + region.color[5][2]=0; + + region.level[6]=140; + region.color[6][0]=0; + region.color[6][1]=255; + region.color[6][2]=0; + + region.level[7]=150; + region.color[7][0]=0; + region.color[7][1]=208; + region.color[7][2]=0; + + region.level[8]=160; + region.color[8][0]=0; + region.color[8][1]=196; + region.color[8][2]=196; + + region.level[9]=170; + region.color[9][0]=0; + region.color[9][1]=148; + region.color[9][2]=255; + + region.level[10]=180; + region.color[10][0]=80; + region.color[10][1]=80; + region.color[10][2]=255; + + region.level[11]=190; + region.color[11][0]=0; + region.color[11][1]=38; + region.color[11][2]=255; + + region.level[12]=200; + region.color[12][0]=142; + region.color[12][1]=63; + region.color[12][2]=255; + + region.level[13]=210; + region.color[13][0]=196; + region.color[13][1]=54; + region.color[13][2]=255; + + region.level[14]=220; + region.color[14][0]=255; + region.color[14][1]=0; + region.color[14][2]=255; + + region.level[15]=230; + region.color[15][0]=255; + region.color[15][1]=194; + region.color[15][2]=204; + + region.levels=16; + + fd=fopen(filename,"r"); + + if (fd==NULL) + fd=fopen(filename,"r"); + + if (fd==NULL) + { + fd=fopen(filename,"w"); + + + + for (x=0; x255) + val[y]=255; + + if (val[y]<0) + val[y]=0; + } + + region.level[x]=val[0]; + region.color[x][0]=val[1]; + region.color[x][1]=val[2]; + region.color[x][2]=val[3]; + x++; + } + + s=fgets(string,80,fd); + } + + fclose(fd); + region.levels=x; + } +} + +void LoadDBMColors(struct site xmtr) +{ + int x, y, ok, val[4]; + char filename[255], string[80], *pointer=NULL, *s=NULL; + FILE *fd=NULL; + + for (x=0; xmtr.filename[x]!='.' && xmtr.filename[x]!=0 && x<250; x++) + filename[x]=xmtr.filename[x]; + + filename[x]='.'; + filename[x+1]='d'; + filename[x+2]='c'; + filename[x+3]='f'; + filename[x+4]=0; + + /* Default values */ + + region.level[0]=0; + region.color[0][0]=255; + region.color[0][1]=0; + region.color[0][2]=0; + + region.level[1]=-10; + region.color[1][0]=255; + region.color[1][1]=128; + region.color[1][2]=0; + + region.level[2]=-20; + region.color[2][0]=255; + region.color[2][1]=165; + region.color[2][2]=0; + + region.level[3]=-30; + region.color[3][0]=255; + region.color[3][1]=206; + region.color[3][2]=0; + + region.level[4]=-40; + region.color[4][0]=255; + region.color[4][1]=255; + region.color[4][2]=0; + + region.level[5]=-50; + region.color[5][0]=184; + region.color[5][1]=255; + region.color[5][2]=0; + + region.level[6]=-60; + region.color[6][0]=0; + region.color[6][1]=255; + region.color[6][2]=0; + + region.level[7]=-70; + region.color[7][0]=0; + region.color[7][1]=208; + region.color[7][2]=0; + + region.level[8]=-80; + region.color[8][0]=0; + region.color[8][1]=196; + region.color[8][2]=196; + + region.level[9]=-90; + region.color[9][0]=0; + region.color[9][1]=148; + region.color[9][2]=255; + + region.level[10]=-100; + region.color[10][0]=80; + region.color[10][1]=80; + region.color[10][2]=255; + + region.level[11]=-110; + region.color[11][0]=0; + region.color[11][1]=38; + region.color[11][2]=255; + + region.level[12]=-120; + region.color[12][0]=142; + region.color[12][1]=63; + region.color[12][2]=255; + + region.level[13]=-130; + region.color[13][0]=196; + region.color[13][1]=54; + region.color[13][2]=255; + + region.level[14]=-140; + region.color[14][0]=255; + region.color[14][1]=0; + region.color[14][2]=255; + + region.level[15]=-150; + region.color[15][0]=255; + region.color[15][1]=194; + region.color[15][2]=204; + + region.levels=16; + + fd=fopen(filename,"r"); + + if (fd==NULL) + fd=fopen(filename,"r"); + + if (fd==NULL) + { + fd=fopen(filename,"w"); + + + for (x=0; x+40) + val[0]=+40; + + region.level[x]=val[0]; + + for (y=1; y<4; y++) + { + if (val[y]>255) + val[y]=255; + + if (val[y]<0) + val[y]=0; + } + + region.color[x][0]=val[1]; + region.color[x][1]=val[2]; + region.color[x][2]=val[3]; + x++; + } + + s=fgets(string,80,fd); + } + + fclose(fd); + region.levels=x; + } +} + + +void DoPathLoss(char *filename, unsigned char geo, unsigned char kml, unsigned char ngs, struct site *xmtr, unsigned char txsites) +{ + /* This function generates a topographic map in Portable Pix Map + (PPM) format based on the content of flags held in the mask[][] + array (only). The image created is rotated counter-clockwise + 90 degrees from its representation in dem[][] so that north + points up and east points right in the image generated. */ + + char mapfile[255], geofile[255], kmlfile[255]; + unsigned width, height, red, green, blue, terrain=0; + unsigned char found, mask, cityorcounty; + int indx, x, y, z, x0, y0, loss, match; + double lat, lon, conversion, one_over_gamma,minwest; + FILE *fd; + + one_over_gamma=1.0/GAMMA; + conversion=255.0/pow((double)(max_elevation-min_elevation),one_over_gamma); + + width=(unsigned)(ippd*ReduceAngle(max_west-min_west)); + height=(unsigned)(ippd*ReduceAngle(max_north-min_north)); + + LoadLossColors(xmtr[0]); + + if (filename[0]==0) + { + strncpy(filename, xmtr[0].filename,254); + filename[strlen(filename)-4]=0; /* Remove .qth */ + } + + y=strlen(filename); + + if (y>4) + { + if (filename[y-1]=='m' && filename[y-2]=='p' && filename[y-3]=='p' && filename[y-4]=='.') + y-=4; + } + + for (x=0; x360.0) + minwest-=360.0; + + north=(double)max_north-dpp; + + if (kml || geo) + south=(double)min_north; /* No bottom legend */ + else + south=(double)min_north-(30.0/ppd); /* 30 pixels for bottom legend */ + + east=(minwest<180.0?-minwest:360.0-min_west); + west=(double)(max_west<180?-max_west:360-max_west); + + + // WriteKML() + //writeKML(xmtr,filename); + + fd=fopen(mapfile,"wb"); + + fprintf(fd,"P6\n%u %u\n255\n",width,(kml?height:height+30)); + if(debug){ + fprintf(stdout,"\nWriting \"%s\" (%ux%u pixmap image)... ",mapfile,width,(kml?height:height+30)); + fflush(stdout); + } + for (y=0, lat=north; y<(int)height; y++, lat=north-(dpp*(double)y)) + { + for (x=0, lon=max_west; x<(int)width; x++, lon=max_west-(dpp*(double)x)) + { + if (lon<0.0) + lon+=360.0; + + for (indx=0, found=0; indx=0 && x0<=mpi && y0>=0 && y0<=mpi) + found=1; + else + indx++; + } + + + + if (found) + { + mask=dem[indx].mask[x0][y0]; + loss=(dem[indx].signal[x0][y0]); + cityorcounty=0; + + match=255; + + red=0; + green=0; + blue=0; + + if (loss<=region.level[0]) + match=0; + else + { + for (z=1; (z=region.level[z-1] && loss=180 && green<=75 && blue<=75 && loss==0) + fprintf(fd,"%c%c%c",255^red,255^green,255^blue); + else + fprintf(fd,"%c%c%c",255,0,0); + + cityorcounty=1; + } + + else if (mask&4) + { + /* County Boundaries: Black */ + + fprintf(fd,"%c%c%c",0,0,0); + + cityorcounty=1; + } + + if (cityorcounty==0) + { + if (loss==0 || (contour_threshold!=0 && loss>abs(contour_threshold))) + { + if (ngs) /* No terrain */ + fprintf(fd,"%c%c%c",255,255,255); + else + { + /* Display land or sea elevation */ + + if (dem[indx].data[x0][y0]==0) + fprintf(fd,"%c%c%c",0,0,170); + else + { + terrain=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion); + fprintf(fd,"%c%c%c",terrain,terrain,terrain); + } + } + } + + else + { + /* Plot path loss in color */ + + if (red!=0 || green!=0 || blue!=0) + fprintf(fd,"%c%c%c",red,green,blue); + + else /* terrain / sea-level */ + { + if (dem[indx].data[x0][y0]==0) + fprintf(fd,"%c%c%c",0,0,170); + else + { + /* Elevation: Greyscale */ + terrain=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion); + fprintf(fd,"%c%c%c",terrain,terrain,terrain); + } + } + } + } + } + + else + { + /* We should never get here, but if */ + /* we do, display the region as black */ + + fprintf(fd,"%c%c%c",0,0,0); + } + } + } + + + + fclose(fd); + if(debug){ + fprintf(stdout,"Done!\n"); + fflush(stdout); + } +} + +void DoSigStr(char *filename, unsigned char geo, unsigned char kml, unsigned char ngs, struct site *xmtr, unsigned char txsites) +{ + /* This function generates a topographic map in Portable Pix Map + (PPM) format based on the signal strength values held in the + signal[][] array. The image created is rotated counter-clockwise + 90 degrees from its representation in dem[][] so that north + points up and east points right in the image generated. */ + + char mapfile[255], geofile[255], kmlfile[255]; + unsigned width, height, terrain, red, green, blue; + unsigned char found, mask, cityorcounty; + int indx, x, y, z=1, x0, y0, signal,match; + double conversion, one_over_gamma, lat, lon, minwest; + FILE *fd; + + one_over_gamma=1.0/GAMMA; + conversion=255.0/pow((double)(max_elevation-min_elevation),one_over_gamma); + + width=(unsigned)(ippd*ReduceAngle(max_west-min_west)); + height=(unsigned)(ippd*ReduceAngle(max_north-min_north)); + + LoadSignalColors(xmtr[0]); + + if (filename[0]==0) + { + strncpy(filename, xmtr[0].filename,254); + filename[strlen(filename)-4]=0; /* Remove .qth */ + } + + y=strlen(filename); + + if (y>4) + { + if (filename[y-1]=='m' && filename[y-2]=='p' && filename[y-3]=='p' && filename[y-4]=='.') + y-=4; + } + + for (x=0; x360.0) + minwest-=360.0; + + north=(double)max_north-dpp; + + + south=(double)min_north; /* No bottom legend */ + + east=(minwest<180.0?-minwest:360.0-min_west); + west=(double)(max_west<180?-max_west:360-max_west); + + // WriteKML() + //writeKML(xmtr,filename); + + fd=fopen(mapfile,"wb"); + + fprintf(fd,"P6\n%u %u\n255\n",width,(kml?height:height+30)); + if(debug){ + fprintf(stdout,"\nWriting \"%s\" (%ux%u pixmap image)... ",mapfile,width,(kml?height:height+30)); + fflush(stdout); + } + for (y=0, lat=north; y<(int)height; y++, lat=north-(dpp*(double)y)) + { + for (x=0, lon=max_west; x<(int)width; x++, lon=max_west-(dpp*(double)x)) + { + if (lon<0.0) + lon+=360.0; + + for (indx=0, found=0; indx=0 && x0<=mpi && y0>=0 && y0<=mpi) + found=1; + else + indx++; + } + + if (found) + { + mask=dem[indx].mask[x0][y0]; + signal=(dem[indx].signal[x0][y0])-100; + cityorcounty=0; + + match=255; + + red=0; + green=0; + blue=0; + + if (signal>=region.level[0]) + match=0; + else + { + for (z=1; (z=region.level[z]) + match=z; + } + } + + if (match=180 && green<=75 && blue<=75) + fprintf(fd,"%c%c%c",255^red,255^green,255^blue); + else + fprintf(fd,"%c%c%c",255,0,0); + + cityorcounty=1; + } + + else if (mask&4) + { + /* County Boundaries: Black */ + + fprintf(fd,"%c%c%c",0,0,0); + + cityorcounty=1; + } + + if (cityorcounty==0) + { + if (contour_threshold!=0 && signal4) + { + if (filename[y-1]=='m' && filename[y-2]=='p' && filename[y-3]=='p' && filename[y-4]=='.') + y-=4; + } + + for (x=0; x360.0) + minwest-=360.0; + + north=(double)max_north-dpp; + + + south=(double)min_north; /* No bottom legend */ + + + east=(minwest<180.0?-minwest:360.0-min_west); + west=(double)(max_west<180?-max_west:360-max_west); + + fd=fopen(mapfile,"wb"); + + fprintf(fd,"P6\n%u %u\n255\n",width,(kml?height:height+30)); + if(debug){ + fprintf(stdout,"\nWriting \"%s\" (%ux%u pixmap image)... ",mapfile,width,(kml?height:height+30)); + fflush(stdout); + } + // WriteKML() + //writeKML(xmtr,filename); + + for (y=0, lat=north; y<(int)height; y++, lat=north-(dpp*(double)y)) + { + for (x=0, lon=max_west; x<(int)width; x++, lon=max_west-(dpp*(double)x)) + { + if (lon<0.0) + lon+=360.0; + + for (indx=0, found=0; indx=0 && x0<=mpi && y0>=0 && y0<=mpi) + found=1; + else + indx++; + } + + if (found) + { + mask=dem[indx].mask[x0][y0]; + dBm=(dem[indx].signal[x0][y0])-200; + cityorcounty=0; + + match=255; + + red=0; + green=0; + blue=0; + + if (dBm>=region.level[0]) + match=0; + else + { + for (z=1; (z=region.level[z]) + match=z; + } + } + + if (match=180 && green<=75 && blue<=75 && dBm!=0) + fprintf(fd,"%c%c%c",255^red,255^green,255^blue); + else + fprintf(fd,"%c%c%c",255,0,0); + + cityorcounty=1; + } + + else if (mask&4) + { + /* County Boundaries: Black */ + + fprintf(fd,"%c%c%c",0,0,0); + + cityorcounty=1; + } + + if (cityorcounty==0) + { + if (contour_threshold!=0 && dBm=360) + ymin-=360; + + ymax=ymin+1; + + while (ymax<0) + ymax+=360; + + while (ymax>=360) + ymax-=360; + + if (ippd==3600) + snprintf(string,19,"%d:%d:%d:%d-hd",x, x+1, ymin, ymax); + else + snprintf(string,16,"%d:%d:%d:%d",x, x+1, ymin, ymax); + + + LoadSDF(string); + } + } + + else + { + for (y=0; y<=width; y++) + for (x=min_lat; x<=max_lat; x++) + { + ymin=max_lon+y; + + while (ymin<0) + ymin+=360; + + while (ymin>=360) + ymin-=360; + + ymax=ymin+1; + + while (ymax<0) + ymax+=360; + + while (ymax>=360) + ymax-=360; + + if (ippd==3600) + snprintf(string,19,"%d:%d:%d:%d-hd",x, x+1, ymin, ymax); + else + snprintf(string,16,"%d:%d:%d:%d",x, x+1, ymin, ymax); + LoadSDF(string); + } + } +} + + + +void LoadUDT(char *filename) +{ + /* This function reads a file containing User-Defined Terrain + features for their addition to the digital elevation model + data used by SPLAT!. Elevations in the UDT file are evaluated + and then copied into a temporary file under /tmp. Then the + contents of the temp file are scanned, and if found to be unique, + are added to the ground elevations described by the digital + elevation data already loaded into memory. */ + + int i, x, y, z, ypix, xpix, tempxpix, tempypix, fd=0, n=0, pixelfound=0; + char input[80], str[3][80], tempname[15], *pointer=NULL, *s=NULL; + double latitude, longitude, height, tempheight; + FILE *fd1=NULL, *fd2=NULL; + + strcpy(tempname,"/tmp/XXXXXX\0"); + + fd1=fopen(filename,"r"); + + if (fd1!=NULL) + { + fd=mkstemp(tempname); + fd2=fopen(tempname,"w"); + + s=fgets(input,78,fd1); + + pointer=strchr(input,';'); + + if (pointer!=NULL) + *pointer=0; + + + while (feof(fd1)==0) + { + /* Parse line for latitude, longitude, height */ + + for (x=0, y=0, z=0; x<78 && input[x]!=0 && z<3; x++) + { + if (input[x]!=',' && y<78) + { + str[z][y]=input[x]; + y++; + } + + else + { + str[z][y]=0; + z++; + y=0; + } + } + + latitude=ReadBearing(str[0]); + longitude=ReadBearing(str[1]); + + if (longitude<0.0) + longitude+=360; + + /* Remove and/or from antenna height string */ + + for (i=0; str[2][i]!=13 && str[2][i]!=10 && str[2][i]!=0; i++); + + str[2][i]=0; + + /* The terrain feature may be expressed in either + feet or meters. If the letter 'M' or 'm' is + discovered in the string, then this is an + indication that the value given is expressed + in meters. Otherwise the height is interpreted + as being expressed in feet. */ + + for (i=0; str[2][i]!='M' && str[2][i]!='m' && str[2][i]!=0 && i<48; i++); + + if (str[2][i]=='M' || str[2][i]=='m') + { + str[2][i]=0; + height=rint(atof(str[2])); + } + + else + { + str[2][i]=0; + height=rint(METERS_PER_FOOT*atof(str[2])); + } + + if (height>0.0) + fprintf(fd2,"%d, %d, %f\n",(int)rint(latitude/dpp), (int)rint(longitude/dpp), height); + + + s=fgets(input,78,fd1); + + pointer=strchr(input,';'); + + if (pointer!=NULL) + *pointer=0; + } + + fclose(fd1); + fclose(fd2); + close(fd); + + + fd1=fopen(tempname,"r"); + fd2=fopen(tempname,"r"); + + y=0; + + n=fscanf(fd1,"%d, %d, %lf", &xpix, &ypix, &height); + + do + { + x=0; + z=0; + + n=fscanf(fd2,"%d, %d, %lf", &tempxpix, &tempypix, &tempheight); + + do + { + if (x>y && xpix==tempxpix && ypix==tempypix) + { + z=1; /* Dupe! */ + + if (tempheight>height) + height=tempheight; + } + + else + { + n=fscanf(fd2,"%d, %d, %lf", &tempxpix, &tempypix, &tempheight); + x++; + } + + } while (feof(fd2)==0 && z==0); + + if (z==0) /* No duplicate found */ + + //fprintf(stdout,"%lf, %lf \n",xpix*dpp, ypix*dpp); + fflush(stdout); + pixelfound = AddElevation(xpix*dpp, ypix*dpp, height); + //fprintf(stdout,"%d \n",pixelfound); + fflush(stdout); + + n=fscanf(fd1,"%d, %d, %lf", &xpix, &ypix, &height); + y++; + + rewind(fd2); + + } while (feof(fd1)==0); + + fclose(fd1); + fclose(fd2); + unlink(tempname); + } + + //else + //fprintf(stderr,"\n*** ERROR: \"%s\": not found!",filename); + + //fprintf(stdout,"\n"); +} + + +int main(int argc, char *argv[]) +{ + int x, y, z=0, min_lat, min_lon, max_lat, max_lon, + rxlat, rxlon, txlat, txlon, west_min, west_max, + north_min, north_max; + + unsigned char LRmap=0, map=0,txsites=0, + topomap=0, geo=0, kml=0, area_mode=0, max_txsites, ngs=0; + + char mapfile[255], elevation_file[255], longley_file[255], terrain_file[255], + string[255], rxfile[255],txfile[255], udt_file[255], rxsite=0, ani_filename[255], + ano_filename[255]; + + double altitude=0.0, altitudeLR=0.0, tx_range=0.0, + rx_range=0.0, deg_range=0.0, deg_limit=0.0, + deg_range_lon; + + struct site tx_site[32], rx_site; + + + + strncpy(ss_version,"1.3.3\0",6); + strncpy(ss_name,"Signal Server\0",14); + + if (argc==1) + { + fprintf(stdout,"\n\t\t -- %s %s options --\n\n",ss_name, ss_version); + fprintf(stdout," -d Directory containing .sdf tiles\n"); + fprintf(stdout," -lat Tx Latitude (decimal degrees)\n"); + fprintf(stdout," -lon Tx Longitude (decimal degrees) Positive 0-360 \n"); + fprintf(stdout," -txh Tx Height (above ground)\n"); + fprintf(stdout," -f Tx Frequency (MHz)\n"); + fprintf(stdout," -erp Tx Effective Radiated Power (Watts)\n"); + fprintf(stdout," -rxh Rx Height(s) (optional. Default=0.1)\n"); + fprintf(stdout," -rt Rx Threshold (dB / dBm / dBuV/m)\n"); + fprintf(stdout," -hp Horizontal Polarisation (default=vertical)\n"); + fprintf(stdout," -gc Ground clutter (feet/meters)\n"); + fprintf(stdout," -udt User defined terrain filename\n"); + fprintf(stdout," -dbm Plot Rxd signal power instead of field strength\n"); + fprintf(stdout," -m Metric units of measurement\n"); + fprintf(stdout," -te Terrain code 1-6 (optional)\n"); + fprintf(stdout," -terdic Terrain dielectric value 2-80 (optional)\n"); + fprintf(stdout," -tercon Terrain conductivity 0.01-0.0001 (optional)\n"); + fprintf(stdout," -cl Climate code 1-6 (optional)\n"); + fprintf(stdout," -o Filename. Required. \n"); + fprintf(stdout," -R Radius (miles/kilometers)\n"); + fprintf(stdout," -res Resolution. Default=1200, Med=600, Low=300)\n"); + fprintf(stdout," -t Terrain background\n"); + fprintf(stdout," -dbg Debug\n\n"); + + + + fflush(stdout); + + return 1; + } + + y=argc-1; + + kml=0; + geo=0; + dbm=0; + gpsav=0; + metric=0; + rxfile[0]=0; + txfile[0]=0; + string[0]=0; + mapfile[0]=0; + clutter=0.0; + forced_erp=-1.0; + forced_freq=0.0; + elevation_file[0]=0; + terrain_file[0]=0; + sdf_path[0]=0; + udt_file[0]=0; + path.length=0; + max_txsites=30; + fzone_clearance=0.6; + contour_threshold=0; + rx_site.lat=91.0; + rx_site.lon=361.0; + longley_file[0]=0; + ano_filename[0]=0; + ani_filename[0]=0; + earthradius=EARTHRADIUS; + max_range=1.0; + + + lat=0; + lon=0; + txh=0; + ngs=1; // no terrain background + kml=1; + + map=1; + LRmap=1; + area_mode=1; + ippd=IPPD; // default resolution + + sscanf("0.1","%lf",&altitudeLR); + + // Defaults + LR.eps_dielect=15.0; // Farmland + LR.sgm_conductivity=0.005; // Farmland + LR.eno_ns_surfref=301.0; + LR.frq_mhz=19.0; // Deliberately too low + LR.radio_climate=5; // continental + LR.pol=1; // vert + LR.conf=0.50; + LR.rel=0.50; + LR.erp=0.0; // will default to Path Loss + + tx_site[0].lat=91.0; + tx_site[0].lon=361.0; + + + for (x=0; x 90 || tx_site[0].lat < -90) + { + fprintf(stdout,"ERROR: Either the lat was missing or out of range!"); + exit(0); + + } + if (tx_site[0].lon > 360 || tx_site[0].lon < 0) + { + fprintf(stdout,"ERROR: Either the lon was missing or out of range!"); + exit(0); + + } + + if (LR.frq_mhz < 20 || LR.frq_mhz > 20000) + { + fprintf(stdout,"ERROR: Either the Frequency was missing or out of range!"); + exit(0); + } + if (LR.erp>2000000) + { + fprintf(stdout,"ERROR: Power was out of range!"); + exit(0); + + } + if (LR.eps_dielect > 80 || LR.eps_dielect < 0.1) + { + fprintf(stdout,"ERROR: Ground Dielectric value out of range!"); + exit(0); + + } + if (LR.sgm_conductivity > 0.01 || LR.sgm_conductivity < 0.000001) + { + fprintf(stdout,"ERROR: Ground conductivity out of range!"); + exit(0); + + } + + if (tx_site[0].alt < 0 || tx_site[0].alt > 60000) + { + fprintf(stdout,"ERROR: Tx altitude above ground was too high!"); + exit(0); + } + if (altitudeLR < 0 || altitudeLR > 60000) + { + fprintf(stdout,"ERROR: Rx altitude above ground was too high!"); + exit(0); + } + + + if (ippd < 300 || ippd > 1200){ + fprintf(stdout,"ERROR: resolution out of range!"); + exit(0); + } + + if(contour_threshold < -200 || contour_threshold > 200) + { + fprintf(stdout,"ERROR: Receiver threshold out of range (-200 / +200)"); + exit(0); + } + + + /* ERROR DETECTION COMPLETE */ + + + + if (metric) + { + altitudeLR/=METERS_PER_FOOT; /* RXH meters --> feet */ + max_range/=KM_PER_MILE; /* RAD kilometers --> miles */ + //altitude/=METERS_PER_FOOT; + tx_site[0].alt/=METERS_PER_FOOT; /* TXH meters --> feet */ + clutter/=METERS_PER_FOOT; /* CLH meters --> feet */ + + } + + /* Ensure a trailing '/' is present in sdf_path */ + + if (sdf_path[0]) + { + x=strlen(sdf_path); + + if (sdf_path[x-1]!='/' && x!=0) + { + sdf_path[x]='/'; + sdf_path[x+1]=0; + } + } + + x=0; + y=0; + + min_lat=70; + max_lat=-70; + + min_lon=(int)floor(tx_site[0].lon); + max_lon=(int)floor(tx_site[0].lon); + + + txlat=(int)floor(tx_site[0].lat); + txlon=(int)floor(tx_site[0].lon); + + if (txlatmax_lat) + max_lat=txlat; + + if (LonDiff(txlon,min_lon)<0.0) + min_lon=txlon; + + if (LonDiff(txlon,max_lon)>=0.0) + max_lon=txlon; + + + if (rxsite) + { + rxlat=(int)floor(rx_site.lat); + rxlon=(int)floor(rx_site.lon); + + if (rxlatmax_lat) + max_lat=rxlat; + + if (LonDiff(rxlon,min_lon)<0.0) + min_lon=rxlon; + + if (LonDiff(rxlon,max_lon)>=0.0) + max_lon=rxlon; + } + + /* Load the required SDF files */ + + LoadTopoData(max_lon, min_lon, max_lat, min_lat); + + if (area_mode || topomap) + { + for (z=0; zdeg_limit) + deg_range=deg_limit; + + if (deg_range_lon>deg_limit) + deg_range_lon=deg_limit; + + north_min=(int)floor(tx_site[z].lat-deg_range); + north_max=(int)floor(tx_site[z].lat+deg_range); + + west_min=(int)floor(tx_site[z].lon-deg_range_lon); + + while (west_min<0) + west_min+=360; + + while (west_min>=360) + west_min-=360; + + west_max=(int)floor(tx_site[z].lon+deg_range_lon); + + while (west_max<0) + west_max+=360; + + while (west_max>=360) + west_max-=360; + + if (north_minmax_lat) + max_lat=north_max; + + if (LonDiff(west_min,min_lon)<0.0) + min_lon=west_min; + + if (LonDiff(west_max,max_lon)>=0.0) + max_lon=west_max; + } + + /* Load any additional SDF files, if required */ + + LoadTopoData(max_lon, min_lon, max_lat, min_lat); + } + + // UDT clutter + LoadUDT (udt_file); + +// Calculate it + if (debug){ + fprintf(stdout,"\nUsing %.5f\n",LR.eps_dielect); + } + + if (area_mode && topomap==0) + { + + PlotLRMap(tx_site[0],altitudeLR,ano_filename); + + } + + if (debug){ + fprintf(stdout,"\nCalc complete\n"); + } + if (map || topomap) + { + + + if (LR.erp==0.0) + DoPathLoss(mapfile,geo,kml,ngs,tx_site,txsites); + else + if (dbm) + DoRxdPwr(mapfile,geo,kml,ngs,tx_site,txsites); + else + DoSigStr(mapfile,geo,kml,ngs,tx_site,txsites); + } + + fprintf(stdout,"|%.5f",north); + fprintf(stdout,"|%.5f",east); + fprintf(stdout,"|%.5f",south); + fprintf(stdout,"|%.5f|",west); + + + fflush(stdout); + + printf("\n"); + + return 0; +}