Imported Upstream version 1.2.1
[debian/splat] / itm.cpp
diff --git a/itm.cpp b/itm.cpp
index 602aa0b6d607046f4cebde3ec1b50d7aec6c62a2..486dcb057e8cd3ac2ecd7a27bd158a3b64dec855 100644 (file)
--- a/itm.cpp
+++ b/itm.cpp
@@ -1,14 +1,16 @@
-/****************************************************************************
-*                                                                           *
-* This file was obtained from ftp://flattop.its.bldrdoc.gov/itm/ITMDLL.cpp  *
-* on Jan. 10, 2004 and is public domain.  It was modified by J. D. McDonald *
-* to remove Microsoft Windows dll-isms and to correct one case where a      *
-* compiler found an ambguity in overloaded calls.  It was further modified  *
-* by John A. Magliacane to remove unused variables, unneeded #includes,     *
-* and to replace pow() statements with explicit multiplications wherever    *
-* possible to increase execution speed and improve accuracy.                *
-*                                                                           *
-****************************************************************************/
+/*****************************************************************************\
+ *                                                                           *
+ *  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
 #include <assert.h>
 #include <string.h>
 
-#define THIRD (1.0/3.0)
+#define THIRD  (1.0/3.0)
 
 using namespace std;
 
 struct tcomplex
-{      double tcreal;
+{
+       double tcreal;
        double tcimag;
 };
 
 struct prop_type
-{      double aref;
+{
+       double aref;
        double dist;
        double hg[2];
        double wn;
@@ -51,14 +55,16 @@ struct prop_type
 };
 
 struct propv_type
-{      double sgc;
+{
+       double sgc;
        int lvar;
        int mdvar;
        int klim;
 };
 
 struct propa_type
-{      double dlsa;
+{
+       double dlsa;
        double dx;
        double ael;
        double ak1;
@@ -106,8 +112,8 @@ double mymax(const double &a, const double &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 */
+       // 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;
@@ -123,25 +129,27 @@ double aknfe(const double &v2)
                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<1.0e-5 || x*w*w*w > 5495.0)
+               /* 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;
        }
@@ -156,6 +164,7 @@ double fht(const double& x, const double& pk)
                        fhtv=(1.0-w)*fhtv+w*(17.372*log(x)-117.0);
                }
        }
+
        return fhtv;
 }
 
@@ -164,8 +173,8 @@ 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;
-       double h0fv, temp; 
        int it;
+       double h0fv;
 
        it=(int)et;
 
@@ -185,10 +194,8 @@ double h0f(double r, double et)
                q=et-it;
 
        /* x=pow(1.0/r,2.0); */
-
-       temp=1.0/r;
-       x=temp*temp;
-
+       x=(1.0/r);
+       x*=x;
        h0fv=4.343*log((a[it-1]*x+b[it-1])*x+1.0);
 
        if (q!=0.0)
@@ -200,9 +207,9 @@ double h0f(double r, double et)
 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};
+       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;
@@ -216,7 +223,7 @@ double ahd(double td)
        return a[i]+b[i]*td+c[i]*log(td);
 }
 
-double adiff(double d, prop_type &prop, propa_type &propa)
+double  adiff( double d, prop_type &prop, propa_type &propa)
 {
        complex<double> prop_zgnd(prop.zgndreal,prop.zgndimag);
        static double wd1, xd1, afo, qk, aht, xht;
@@ -253,7 +260,7 @@ double adiff(double d, prop_type &prop, propa_type &propa)
                adiffv=0.0;
        }
 
-       else
+       else
        {
                th=propa.tha+d*prop.gme;
                ds=d-propa.dla;
@@ -273,7 +280,7 @@ double adiff(double d, prop_type &prop, propa_type &propa)
        return adiffv;
 }
 
-double ascat( double d, prop_type &prop, propa_type &propa)
+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;
@@ -315,7 +322,6 @@ double ascat( double d, prop_type &prop, propa_type &propa)
                        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;
@@ -326,12 +332,7 @@ double ascat( double d, prop_type &prop, propa_type &propa)
                        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)); */
-
-                               temp=((1.0+1.4142/r1)*(1.0+1.4142/r2));
-                               h0=et*h0+(1.0-et)*4.343*log((temp*temp)*(r1+r2)/(r1+r2+2.8284));
-                       }
+                               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;
@@ -340,13 +341,15 @@ double ascat( double d, prop_type &prop, propa_type &propa)
                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;
+
+               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 qerfi (double q)
 {
        double x, t, v;
        double c0=2.515516698;
@@ -367,7 +370,8 @@ double qerfi(double q)
        return v;
 }
 
-void qlrps(double fmhz, double zsys, double en0, int ipol, double eps, double sgm, prop_type &prop)
+void qlrps( double fmhz, double zsys, double en0, int ipol, double eps, double sgm, prop_type &prop)
+
 {
        double gma=157e-9;
 
@@ -389,7 +393,7 @@ void qlrps(double fmhz, double zsys, double en0, int ipol, double eps, double sg
        prop.zgndimag=prop_zgnd.imag();
 }
 
-double abq_alos (complex<double> r)
+double abq_alos(complex<double> r)
 {
        return r.real()*r.real()+r.imag()*r.imag();
 }
@@ -427,8 +431,8 @@ double alos(double d, prop_type &prop, propa_type &propa)
                        q=3.14-2.4649/q;
 
                alosv=(-4.343*log(abq_alos(complex<double>(cos(q),-sin(q))+r))-alosv)*wls+alosv;
+        }
 
-       }
        return alosv;
 }
 
@@ -440,7 +444,7 @@ void qlra(int kst[], int klimx, int mdvarx, prop_type &prop, propv_type &propv)
        for (int j=0; j<2; ++j)
        {
                if (kst[j]<=0)
-                        prop.he[j]=prop.hg[j];
+                       prop.he[j]=prop.hg[j];
                else
                {
                        q=4.0;
@@ -461,7 +465,7 @@ void qlra(int kst[], int klimx, int mdvarx, prop_type &prop, propv_type &propv)
 
        prop.mdp=1;
        propv.lvar=mymax(propv.lvar,3);
-  
+
        if (mdvarx>=0)
        {
                propv.mdvar=mdvarx;
@@ -475,194 +479,8 @@ void qlra(int kst[], int klimx, int mdvarx, prop_type &prop, propv_type &propv)
        }
 }
 
-void freds_lrprop (double d, prop_type &prop, propa_type &propa)
+void lrprop (double d, prop_type &prop, propa_type &propa)  // PaulM_lrprop
 {
-       /* freds_lrprop */
-
-       static bool wlos, wscat;
-       static double dmin, xae;
-       complex<double> prop_zgnd(prop.zgndreal,prop.zgndimag);
-       double a0, a1, a2, a3, a4, a5, a6;
-       double d0, d1, d2, d3, d4, d5, d6;
-       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*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.dist=0;
-
-               else if (prop.mdp>0)
-               {
-                       prop.mdp=0;
-                       prop.dist=0;
-               }
-
-               if ((prop.dist>0.0) || (prop.mdp<0.0))
-               {
-                       if (prop.dist>1000e3)
-                               prop.kwx=mymax(prop.kwx,1);
-
-                       if (prop.dist<dmin)
-                               prop.kwx=mymax(prop.kwx,3);
-
-                       if (prop.dist<1e3 || prop.dist>2000e3)
-                               prop.kwx=4;
-               }
-       }
-
-       else
-       {
-               prop.dist=d;
-
-               if (prop.dist>0.0)
-               {
-                       if (prop.dist>1000e3)
-                               prop.kwx=mymax(prop.kwx,1);
-
-                       if (prop.dist<dmin)
-                               prop.kwx=mymax(prop.kwx,3);
-
-                       if (prop.dist<1e3 || prop.dist>2000e3)
-                               prop.kwx=4;
-               }
-       }
-
-       if (prop.dist<propa.dlsa)
-       {
-               if (!wlos)
-               {
-                       /* Coefficients for the line-of-sight range */
-
-                       q=alos(0.0,prop,propa);
-                       d2=propa.dlsa;
-                       a2=propa.aed+d2*propa.emd;
-                       d0=1.908*prop.wn*prop.he[0]*prop.he[1];
-
-                       if (propa.aed>=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);
-        
-                       if (d0<d1)
-                       {
-                               a0=alos(d0,prop,propa);
-                               q=log(d2/d0);
-                               propa.ak2=mymax(0.0,((d2-d0)*(a1-a0)-(d1-d0)*(a2-a0))/((d2-d0)*log(d1-d0)-(d1-d0)*q));
-
-                               if (propa.ak2<=0.0 && propa.aed<0.0)
-                               {
-                                       propa.ak2=0.0;
-                                       propa.ak1=(a2-a1)/(d2-d1);
-
-                                       if (propa.ak1<=0.0)
-                                               propa.ak2=propa.emd;
-                               }
-
-                               else
-                               {
-                                       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;
-                                       }
-                               }
-                               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);
-       }
-
-       else
-       {
-               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=propa.emd;
-                               propa.aes=propa.aed;
-                               propa.dx=10e6;
-                       }
-
-                       else
-                       {
-                               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;
-                       }
-
-                       wscat=true;
-               }
-
-               if (prop.dist<=propa.dx)
-                       prop.aref=propa.aed+propa.emd*prop.dist;
-               else
-                       prop.aref=propa.aes+propa.ems*prop.dist;
-       }
-
-       prop.aref=FORTRAN_DIM(prop.aref,0.0);
-}
-
-void lrprop (double d, prop_type &prop, propa_type &propa)
-{
-       /* PaulM_lrprop */
        static bool wlos, wscat;
        static double dmin, xae;
        complex<double> prop_zgnd(prop.zgndreal,prop.zgndimag);
@@ -673,7 +491,7 @@ void lrprop (double d, prop_type &prop, propa_type &propa)
        int j;
 
        if (prop.mdp!=0)
-       {
+       {
                for (j=0; j<2; j++)
                        propa.dls[j]=sqrt(2.0*prop.he[j]/prop.gme);
 
@@ -691,20 +509,24 @@ void lrprop (double d, prop_type &prop, propa_type &propa)
                                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] )
+                       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;
+
+               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); -- JDM made argument 2 a double */
-               xae=pow(prop.wn*(prop.gme*prop.gme),-THIRD);  /* No 2nd pow() */
+               /* 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);
@@ -766,29 +588,20 @@ void lrprop (double d, prop_type &prop, propa_type &propa)
                                        if (propa.ak1<0.0)
                                        {
                                                propa.ak1=0.0;
-                                               propa.ak2=FORTRAN_DIM(a2,a0)/q;
+                                               propa.ak2=FORTRAN_DIM(a2,a0)/q;
 
                                                if (propa.ak2==0.0)
                                                        propa.ak1=propa.emd;
                                        }
                                }
-
-                               else
-                               {
-                                       propa.ak2=0.0;
-                                       propa.ak1=(a2-a1)/(d2-d1);
-
-                                       if (propa.ak1<=0.0)
-                                               propa.ak1=propa.emd;
-                               }
                        }
-       
-                       else
+
+                       if (!wq)
                        {
-                               propa.ak1=(a2-a1)/(d2-d1);
+                               propa.ak1=FORTRAN_DIM(a2,a1)/(d2-d1);
                                propa.ak2=0.0;
 
-                               if (propa.ak1<=0.0)
+                               if (propa.ak1==0.0)
                                        propa.ak1=propa.emd;
                        }
 
@@ -798,12 +611,11 @@ void lrprop (double d, prop_type &prop, propa_type &propa)
 
                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)
+               if (!wscat)
                { 
                        q=ascat(0.0,prop,propa);
                        d5=propa.dla+200e3;
@@ -815,6 +627,7 @@ void lrprop (double d, prop_type &prop, propa_type &propa)
                        {
                                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;
                        }
 
@@ -837,10 +650,10 @@ void lrprop (double d, prop_type &prop, propa_type &propa)
        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)
+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;
@@ -852,13 +665,10 @@ double curve (double const &c1, double const &c2, double const &x1,
        return (c1+c2/(1.0+temp1))*temp2/(1.0+temp2);
 }
 
-double avar(double zzt, double zzl, double zzc, prop_type &prop, propv_type &propv)
+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;
+       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};
@@ -885,7 +695,7 @@ double avar(double zzt, double zzl, double zzc, prop_type &prop, propv_type &pro
        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, temp1, temp2;
+       double sgt, yr;
        int temp_klim=propv.klim-1;
 
        if (propv.lvar>0)
@@ -895,16 +705,16 @@ double avar(double zzt, double zzl, double zzc, prop_type &prop, propv_type &pro
                        default:
                        if (propv.klim<=0 || propv.klim>7)
                        {
-                               propv.klim=5;
-                               temp_klim=4;
+                               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];
+                       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];
@@ -980,9 +790,8 @@ double avar(double zzt, double zzl, double zzc, prop_type &prop, propv_type &pro
                else
                {
                        /* vs0=pow(5.0+3.0*exp(-de/100e3),2.0); */
-                       temp1=(5.0+3.0*exp(-de/100e3));
-                       vs0=temp1*temp1;
-
+                       vs0=(5.0+3.0*exp(-de/100e3));
+                       vs0*=vs0;
                }
 
                propv.lvar=0;
@@ -1020,11 +829,7 @@ double avar(double zzt, double zzl, double zzc, prop_type &prop, propv_type &pro
                sgt=sgtd+tgtd/zt;
 
        /* vs=vs0+pow(sgt*zt,2.0)/(rt+zc*zc)+pow(sgl*zl,2.0)/(rl+zc*zc); */
-
-       temp1=sgt*zt;
-       temp2=sgl*zl;
-
-       vs=vs0+(temp1*temp1)/(rt+zc*zc)+(temp2*temp2)/(rl+zc*zc);
+       vs=vs0+(sgt*zt*sgt*zt)/(rt+zc*zc)+(sgl*zl*sgl*zl)/(rl+zc*zc);
 
        if (kdv==0)
        {
@@ -1058,7 +863,7 @@ double avar(double zzt, double zzl, double zzc, prop_type &prop, propv_type &pro
        return avarv;
 }
 
-void hzns (double pfl[], prop_type &prop)
+void hzns(double pfl[], prop_type &prop)
 {
        bool wq;
        int np;
@@ -1110,6 +915,7 @@ void hzns (double pfl[], prop_type &prop)
 }
   
 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;
@@ -1149,8 +955,8 @@ void z1sq1 (double z[], const double &x1, const double &x2, double& z0, double&
 
 double qtile (const int &nn, double a[], const int &ir)
 {
-       double q=0.0, r; /* q initialization -- KD2BD */
-       int m, n, i, j, j1=0, i0=0, k;  /* more initializations -- KD2BD */
+       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;
 
@@ -1174,7 +980,6 @@ double qtile (const int &nn, double a[], const int &ir)
 
                if (i>n)
                        i=n;
-
                j=j1;
 
                while (j>=m && a[j]<=q)
@@ -1210,7 +1015,7 @@ double qtile (const int &nn, double a[], const int &ir)
                }
 
                else
-               done=true;
+                       done=true;
        }
 
        return q;
@@ -1294,10 +1099,10 @@ double d1thx(double pfl[], const double &x1, const double &x2)
        return d1thxv;
 }
 
-void qlrpfl(double pfl[], int klimx, int mdvarx, prop_type &prop, propa_type &propa, propv_type &propv)
+void qlrpfl (double pfl[], int klimx, int mdvarx, prop_type &prop, propa_type &propa, propv_type &propv)
 {
        int np, j;
-       double xl[2], q, za, zb, temp;
+       double xl[2], q, za, zb;
 
        prop.dist=pfl[0]*pfl[1];
        np=(int)pfl[0];
@@ -1323,8 +1128,7 @@ void qlrpfl(double pfl[], int klimx, int mdvarx, prop_type &prop, propa_type &pr
                if (q<=prop.dist)
                {
                        /* q=pow(prop.dist/q,2.0); */
-                       temp=prop.dist/q;
-                       q=temp*temp;
+                       q=((prop.dist/q)*(prop.dist/q));
 
                        for (j=0; j<2; j++)
                        {
@@ -1337,7 +1141,7 @@ void qlrpfl(double pfl[], int klimx, int mdvarx, prop_type &prop, propa_type &pr
                {
                        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
@@ -1376,35 +1180,201 @@ double deg2rad(double d)
 //********************************************************
 
 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;
 
-       pol:
-               0-Horizontal, 1-Vertical
+       double zsys=0;
+       double zc, zr;
+       double eno, enso, q;
+       long ja, jb, i, np;
+       double dkm, xkm;
+       double fs;
 
-       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
+       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]);
 
-       conf, rel: .01 to .99
+               jb=np-ja+6;
 
-       elev[]: [num points - 1], [delta dist(meters)],
-               [height(meters) point 1], ..., [height(meters) point n]
+               for (i=ja-1; i<jb; ++i)
+                       zsys+=elev[i];
 
-       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.
+               zsys/=(jb-ja+1);
+               q=eno;
+       }
 
-*****************************************************************************/
+       propv.mdvar=12;
+       qlrps(frq_mhz,zsys,q,pol,eps_dielect,sgm_conductivity,prop);
+       qlrpfl(elev,propv.klim,propv.mdvar,prop,propa,propv);
+       fs=32.45+20.0*log10(frq_mhz)+20.0*log10(prop.dist/1000.0);
+       q=prop.dist-propa.dla;
+
+       if (int(q)<0.0)
+               strcpy(strmode,"Line-Of-Sight Mode");
+
+       else
+       {
+               if (int(q)==0.0)
+                       strcpy(strmode,"Single Horizon");
+
+               else if (int(q)>0.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; i<jb; ++i)
+                       zsys+=elev[i];
+
+               zsys/=(jb-ja+1);
+               q=eno;
+       }
+
+       propv.mdvar=12;
+       qlrps(frq_mhz,zsys,q,pol,eps_dielect,sgm_conductivity,prop);
+       qlrpfl(elev,propv.klim,propv.mdvar,prop,propa,propv);
+       fs=32.45+20.0*log10(frq_mhz)+20.0*log10(prop.dist/1000.0);
+       deltaH=prop.dh;
+       q=prop.dist-propa.dla;
+
+       if (int(q)<0.0)
+               propmode=0;  // Line-Of-Sight Mode
+       else
+       {
+               if (int(q)==0.0)
+                       propmode=4;  // Single Horizon
+
+               else if (int(q)>0.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;
@@ -1432,7 +1402,9 @@ void point_to_point(double elev[], double tht_m, double rht_m, double eps_dielec
 
        if (q<=0.0)
        {
-               ja=(long)(3.0+0.1*elev[0]);  /* KD2BD added (long) */
+               /* ja = 3.0 + 0.1 * elev[0]; */
+               ja=(long)(3.0+0.1*elev[0]);
+
                jb=np-ja+6;
 
                for (i=ja-1; i<jb; ++i)
@@ -1446,6 +1418,7 @@ void point_to_point(double elev[], double tht_m, double rht_m, double eps_dielec
        qlrps(frq_mhz,zsys,q,pol,eps_dielect,sgm_conductivity,prop);
        qlrpfl(elev,propv.klim,propv.mdvar,prop,propa,propv);
        fs=32.45+20.0*log10(frq_mhz)+20.0*log10(prop.dist/1000.0);
+       deltaH=prop.dh;
        q=prop.dist-propa.dla;
 
        if (int(q)<0.0)
@@ -1465,10 +1438,11 @@ void point_to_point(double elev[], double tht_m, double rht_m, double eps_dielec
                        strcat(strmode, ", Troposcatter Dominant");
        }
 
-       dbloss=avar(zr,0.0,zc,prop,propv)+fs;
+       dbloss=avar(zr,0.0,zc,prop,propv)+fs;      //avar(time,location,confidence)
        errnum=prop.kwx;
 }
 
+
 //********************************************************
 //* Area Mode Calculations                               *
 //********************************************************
@@ -1478,14 +1452,13 @@ void area(long ModVar, double deltaH, double tht_m, double rht_m, double dist_km
        // 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"
+       //         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.
@@ -1509,16 +1482,16 @@ void area(long ModVar, double deltaH, double tht_m, double rht_m, double dist_km
 
        kst[0]=(int)TSiteCriteria;
        kst[1]=(int)RSiteCriteria;
-       zt=qerfi(pctTime/100.0);
-       zl=qerfi(pctLoc/100.0);
-       zc=qerfi(pctConf/100.0);
+       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; -KD2BD replaced below */
+       /* propv.klim = (__int32) radio_climate; */
        propv.klim=(long)radio_climate;
        prop.ens=eno;
        prop.kwx=0;
@@ -1534,9 +1507,26 @@ void area(long ModVar, double deltaH, double tht_m, double rht_m, double dist_km
        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;
+}
+