git.gag.com Git - debian/splat/blob - itm.cpp

   1 /*****************************************************************************\
   2  *                                                                           *
   3  *  The following code was derived from public domain ITM code available     *
   4  *  at ftp://flattop.its.bldrdoc.gov/itm/ITMDLL.cpp that was released on     *
   5  *  June 26, 2007.  It was modified to remove Microsoft Windows "dll-isms",  *
   6  *  redundant and unnecessary #includes, redundant and unnecessary { }'s,    *
   7  *  to initialize uninitialized variables, type cast some variables,         *
   8  *  re-format the code for easier reading, and to replace pow() function     *
   9  *  calls with explicit multiplications wherever possible to increase        *
  10  *  execution speed and improve computational accuracy.                      *
  11  *                                                                           *
  12 \*****************************************************************************/
  13
  14
  15 // *************************************
  16 // C++ routines for this program are taken from
  17 // a translation of the FORTRAN code written by
  18 // U.S. Department of Commerce NTIA/ITS
  19 // Institute for Telecommunication Sciences
  20 // *****************
  21 // Irregular Terrain Model (ITM) (Longley-Rice)
  22 // *************************************
  23
  24 #include <math.h>
  25 #include <complex>
  26 #include <assert.h>
  27 #include <string.h>
  28
  29 #define THIRD  (1.0/3.0)
  30
  31 using namespace std;
  32
  33 struct tcomplex
  34 {
  35         double tcreal;
  36         double tcimag;
  37 };
  38
  39 struct prop_type
  40 {
  41         double aref;
  42         double dist;
  43         double hg[2];
  44         double wn;
  45         double dh;
  46         double ens;
  47         double gme;
  48         double zgndreal;
  49         double zgndimag;
  50         double he[2];
  51         double dl[2];
  52         double the[2];
  53         int kwx;
  54         int mdp;
  55 };
  56
  57 struct propv_type
  58 {
  59         double sgc;
  60         int lvar;
  61         int mdvar;
  62         int klim;
  63 };
  64
  65 struct propa_type
  66 {
  67         double dlsa;
  68         double dx;
  69         double ael;
  70         double ak1;
  71         double ak2;
  72         double aed;
  73         double emd;
  74         double aes;
  75         double ems;
  76         double dls[2];
  77         double dla;
  78         double tha;
  79 };
  80
  81 int mymin(const int &i, const int &j)
  82 {
  83         if (i<j)
  84                 return i;
  85         else
  86                 return j;
  87 }
  88
  89 int mymax(const int &i, const int &j)
  90 {
  91         if (i>j)
  92                 return i;
  93         else
  94                 return j;
  95 }
  96
  97 double mymin(const double &a, const double &b)
  98 {
  99         if (a<b)
 100                 return a;
 101         else
 102                 return b;
 103 }
 104
 105 double mymax(const double &a, const double &b)
 106 {
 107         if (a>b)
 108                 return a;
 109         else
 110                 return b;
 111 }
 112
 113 double FORTRAN_DIM(const double &x, const double &y)
 114 {
 115         // This performs the FORTRAN DIM function.
 116         // result is x-y if x is greater than y; otherwise result is 0.0
 117
 118         if (x>y)
 119                 return x-y;
 120         else
 121                 return 0.0;
 122 }
 123
 124 double aknfe(const double &v2)
 125 {
 126         double a;
 127
 128         if (v2<5.76)
 129                 a=6.02+9.11*sqrt(v2)-1.27*v2;
 130         else
 131                 a=12.953+4.343*log(v2);
 132
 133         return a;
 134 }
 135
 136 double fht(const double& x, const double& pk)
 137 {
 138         double w, fhtv;
 139         if (x<200.0)
 140         {
 141                 w=-log(pk);
 142
 143                 /* if (pk < 1e-5 || x*pow(w,3.0) > 5495.0 ) */
 144
 145                 if (pk < 1e-5 || (x*w*w*w) > 5495.0 )
 146                 {
 147                         fhtv=-117.0;
 148
 149                         if (x>1.0)
 150                                 fhtv=17.372*log(x)+fhtv;
 151                 }
 152
 153                 else
 154                         fhtv=2.5e-5*x*x/pk-8.686*w-15.0;
 155         }
 156
 157         else
 158         {
 159                 fhtv=0.05751*x-4.343*log(x);
 160
 161                 if (x<2000.0)
 162                 {
 163                         w=0.0134*x*exp(-0.005*x);
 164                         fhtv=(1.0-w)*fhtv+w*(17.372*log(x)-117.0);
 165                 }
 166         }
 167
 168         return fhtv;
 169 }
 170
 171 double h0f(double r, double et)
 172 {
 173         double a[5]={25.0, 80.0, 177.0, 395.0, 705.0};
 174         double b[5]={24.0, 45.0,  68.0,  80.0, 105.0};
 175         double q, x;
 176         int it;
 177         double h0fv;
 178
 179         it=(int)et;
 180
 181         if (it<=0)
 182         {
 183                 it=1;
 184                 q=0.0;
 185         }
 186
 187         else if (it>=5)
 188         {
 189                 it=5;
 190                 q=0.0;
 191         }
 192
 193         else
 194                 q=et-it;
 195
 196         /* x=pow(1.0/r,2.0); */
 197         x=(1.0/r);
 198         x*=x;
 199         h0fv=4.343*log((a[it-1]*x+b[it-1])*x+1.0);
 200
 201         if (q!=0.0)
 202                 h0fv=(1.0-q)*h0fv+q*4.343*log((a[it]*x+b[it])*x+1.0);
 203
 204         return h0fv;
 205 }
 206
 207 double ahd(double td)
 208 {
 209         int i;
 210         double a[3] = {   133.4,    104.6,     71.8};
 211         double b[3] = {0.332e-3, 0.212e-3, 0.157e-3};
 212         double c[3] = {  -4.343,   -1.086,    2.171};
 213
 214         if (td<=10e3)
 215                 i=0;
 216
 217         else if (td<=70e3)
 218                 i=1;
 219
 220         else
 221                 i=2;
 222
 223         return a[i]+b[i]*td+c[i]*log(td);
 224 }
 225
 226 double  adiff( double d, prop_type &prop, propa_type &propa)
 227 {
 228         complex<double> prop_zgnd(prop.zgndreal,prop.zgndimag);
 229         static double wd1, xd1, afo, qk, aht, xht;
 230         double a, q, pk, ds, th, wa, ar, wd, adiffv;
 231
 232         if (d==0)
 233         {
 234                 q=prop.hg[0]*prop.hg[1];
 235                 qk=prop.he[0]*prop.he[1]-q;
 236
 237                 if (prop.mdp<0.0)
 238                         q+=10.0;
 239
 240                 wd1=sqrt(1.0+qk/q);
 241                 xd1=propa.dla+propa.tha/prop.gme;
 242                 q=(1.0-0.8*exp(-propa.dlsa/50e3))*prop.dh;
 243                 q*=0.78*exp(-pow(q/16.0,0.25));
 244                 afo=mymin(15.0,2.171*log(1.0+4.77e-4*prop.hg[0]*prop.hg[1]*prop.wn*q));
 245                 qk=1.0/abs(prop_zgnd);
 246                 aht=20.0;
 247                 xht=0.0;
 248
 249                 for (int j=0; j<2; ++j)
 250                 {
 251                         /* a=0.5*pow(prop.dl[j],2.0)/prop.he[j]; */
 252                         a=0.5*(prop.dl[j]*prop.dl[j])/prop.he[j];
 253                         wa=pow(a*prop.wn,THIRD);
 254                         pk=qk/wa;
 255                         q=(1.607-pk)*151.0*wa*prop.dl[j]/a;
 256                         xht+=q;
 257                         aht+=fht(q,pk);
 258                 }
 259
 260                 adiffv=0.0;
 261         }
 262
 263         else
 264         {
 265                 th=propa.tha+d*prop.gme;
 266                 ds=d-propa.dla;
 267                 /* q=0.0795775*prop.wn*ds*pow(th,2.0); */
 268                 q=0.0795775*prop.wn*ds*th*th;
 269                 adiffv=aknfe(q*prop.dl[0]/(ds+prop.dl[0]))+aknfe(q*prop.dl[1]/(ds+prop.dl[1]));
 270                 a=ds/th;
 271                 wa=pow(a*prop.wn,THIRD);
 272                 pk=qk/wa;
 273                 q=(1.607-pk)*151.0*wa*th+xht;
 274                 ar=0.05751*q-4.343*log(q)-aht;
 275                 q=(wd1+xd1/d)*mymin(((1.0-0.8*exp(-d/50e3))*prop.dh*prop.wn),6283.2);
 276                 wd=25.1/(25.1+sqrt(q));
 277                 adiffv=ar*wd+(1.0-wd)*adiffv+afo;
 278         }
 279
 280         return adiffv;
 281 }
 282
 283 double  ascat( double d, prop_type &prop, propa_type &propa)
 284 {
 285         static double ad, rr, etq, h0s;
 286         double h0, r1, r2, z0, ss, et, ett, th, q;
 287         double ascatv, temp;
 288
 289         if (d==0.0)
 290         {
 291                 ad=prop.dl[0]-prop.dl[1];
 292                 rr=prop.he[1]/prop.he[0];
 293
 294                 if (ad<0.0)
 295                 {
 296                         ad=-ad;
 297                         rr=1.0/rr;
 298                 }
 299
 300                 etq=(5.67e-6*prop.ens-2.32e-3)*prop.ens+0.031;
 301                 h0s=-15.0;
 302                 ascatv=0.0;
 303         }
 304
 305         else
 306         {
 307                 if (h0s>15.0)
 308                         h0=h0s;
 309                 else
 310                 {
 311                         th=prop.the[0]+prop.the[1]+d*prop.gme;
 312                         r2=2.0*prop.wn*th;
 313                         r1=r2*prop.he[0];
 314                         r2*=prop.he[1];
 315
 316                         if (r1<0.2 && r2<0.2)
 317                                 return 1001.0;  // <==== early return
 318
 319                         ss=(d-ad)/(d+ad);
 320                         q=rr/ss;
 321                         ss=mymax(0.1,ss);
 322                         q=mymin(mymax(0.1,q),10.0);
 323                         z0=(d-ad)*(d+ad)*th*0.25/d;
 324                         /* et=(etq*exp(-pow(mymin(1.7,z0/8.0e3),6.0))+1.0)*z0/1.7556e3; */
 325                         temp=mymin(1.7,z0/8.0e3);
 326                         temp=temp*temp*temp*temp*temp*temp;
 327                         et=(etq*exp(-temp)+1.0)*z0/1.7556e3;
 328
 329                         ett=mymax(et,1.0);
 330                         h0=(h0f(r1,ett)+h0f(r2,ett))*0.5;
 331                         h0+=mymin(h0,(1.38-log(ett))*log(ss)*log(q)*0.49);
 332                         h0=FORTRAN_DIM(h0,0.0);
 333
 334                         if (et<1.0)
 335                                 h0=et*h0+(1.0-et)*4.343*log(pow((1.0+1.4142/r1) *(1.0+1.4142/r2),2.0)*(r1+r2)/(r1+r2+2.8284));
 336
 337                         if (h0>15.0 && h0s>=0.0)
 338                                 h0=h0s;
 339                 }
 340
 341                 h0s=h0;
 342                 th=propa.tha+d*prop.gme;
 343                 /* ascatv=ahd(th*d)+4.343*log(47.7*prop.wn*pow(th,4.0))-0.1*(prop.ens-301.0)*exp(-th*d/40e3)+h0; */
 344
 345                 ascatv=ahd(th*d)+4.343*log(47.7*prop.wn*th*th*th*th)-0.1*(prop.ens-301.0)*exp(-th*d/40e3)+h0;
 346
 347         }
 348
 349         return ascatv;
 350 }
 351
 352 double qerfi (double q)
 353 {
 354         double x, t, v;
 355         double c0=2.515516698;
 356         double c1=0.802853;
 357         double c2=0.010328;
 358         double d1=1.432788;
 359         double d2=0.189269;
 360         double d3=0.001308;
 361
 362         x=0.5-q;
 363         t=mymax(0.5-fabs(x),0.000001);
 364         t=sqrt(-2.0*log(t));
 365         v=t-((c2*t+c1)*t+c0)/(((d3*t+d2)*t+d1)*t+1.0);
 366
 367         if (x<0.0)
 368                 v=-v;
 369
 370         return v;
 371 }
 372
 373 void qlrps( double fmhz, double zsys, double en0, int ipol, double eps, double sgm, prop_type &prop)
 374
 375 {
 376         double gma=157e-9;
 377
 378         prop.wn=fmhz/47.7;
 379         prop.ens=en0;
 380
 381         if (zsys!=0.0)
 382                 prop.ens*=exp(-zsys/9460.0);
 383
 384         prop.gme=gma*(1.0-0.04665*exp(prop.ens/179.3));
 385         complex<double> zq, prop_zgnd(prop.zgndreal,prop.zgndimag);
 386         zq=complex<double> (eps,376.62*sgm/prop.wn);
 387         prop_zgnd=sqrt(zq-1.0);
 388
 389         if (ipol!=0.0)
 390                 prop_zgnd=prop_zgnd/zq;
 391
 392         prop.zgndreal=prop_zgnd.real();
 393         prop.zgndimag=prop_zgnd.imag();
 394 }
 395
 396 double abq_alos(complex<double> r)
 397 {
 398         return r.real()*r.real()+r.imag()*r.imag();
 399 }
 400
 401 double alos(double d, prop_type &prop, propa_type &propa)
 402 {
 403         complex<double> prop_zgnd(prop.zgndreal,prop.zgndimag);
 404         static double wls;
 405         complex<double> r;
 406         double s, sps, q;
 407         double alosv;
 408
 409         if (d==0.0)
 410         {
 411                 wls=0.021/(0.021+prop.wn*prop.dh/mymax(10e3,propa.dlsa));
 412                 alosv=0.0;
 413         }
 414
 415         else
 416         {
 417                 q=(1.0-0.8*exp(-d/50e3))*prop.dh;
 418                 s=0.78*q*exp(-pow(q/16.0,0.25));
 419                 q=prop.he[0]+prop.he[1];
 420                 sps=q/sqrt(d*d+q*q);
 421                 r=(sps-prop_zgnd)/(sps+prop_zgnd)*exp(-mymin(10.0,prop.wn*s*sps));
 422                 q=abq_alos(r);
 423
 424                 if (q<0.25 || q<sps)
 425                         r=r*sqrt(sps/q);
 426
 427                 alosv=propa.emd*d+propa.aed;
 428                 q=prop.wn*prop.he[0]*prop.he[1]*2.0/d;
 429
 430                 if (q>1.57)
 431                         q=3.14-2.4649/q;
 432
 433                 alosv=(-4.343*log(abq_alos(complex<double>(cos(q),-sin(q))+r))-alosv)*wls+alosv;
 434          }
 435
 436         return alosv;
 437 }
 438
 439
 440 void qlra(int kst[], int klimx, int mdvarx, prop_type &prop, propv_type &propv)
 441 {
 442         double q;
 443
 444         for (int j=0; j<2; ++j)
 445         {
 446                 if (kst[j]<=0)
 447                         prop.he[j]=prop.hg[j];
 448                 else
 449                 {
 450                         q=4.0;
 451
 452                         if (kst[j]!=1)
 453                                 q=9.0;
 454
 455                         if (prop.hg[j]<5.0)
 456                                 q*=sin(0.3141593*prop.hg[j]);
 457
 458                         prop.he[j]=prop.hg[j]+(1.0+q)*exp(-mymin(20.0,2.0*prop.hg[j]/mymax(1e-3,prop.dh)));
 459                 }
 460
 461                 q=sqrt(2.0*prop.he[j]/prop.gme);
 462                 prop.dl[j]=q*exp(-0.07*sqrt(prop.dh/mymax(prop.he[j],5.0)));
 463                 prop.the[j]=(0.65*prop.dh*(q/prop.dl[j]-1.0)-2.0*prop.he[j])/q;
 464         }
 465
 466         prop.mdp=1;
 467         propv.lvar=mymax(propv.lvar,3);
 468
 469         if (mdvarx>=0)
 470         {
 471                 propv.mdvar=mdvarx;
 472                 propv.lvar=mymax(propv.lvar,4);
 473         }
 474
 475         if (klimx>0)
 476         {
 477                 propv.klim=klimx;
 478                 propv.lvar=5;
 479         }
 480 }
 481
 482 void lrprop (double d, prop_type &prop, propa_type &propa)  // PaulM_lrprop
 483 {
 484         static bool wlos, wscat;
 485         static double dmin, xae;
 486         complex<double> prop_zgnd(prop.zgndreal,prop.zgndimag);
 487         double a0, a1, a2, a3, a4, a5, a6;
 488         double d0, d1, d2, d3, d4, d5, d6;
 489         bool wq;
 490         double q;
 491         int j;
 492
 493         if (prop.mdp!=0)
 494         {
 495                 for (j=0; j<2; j++)
 496                         propa.dls[j]=sqrt(2.0*prop.he[j]/prop.gme);
 497
 498                 propa.dlsa=propa.dls[0]+propa.dls[1];
 499                 propa.dla=prop.dl[0]+prop.dl[1];
 500                 propa.tha=mymax(prop.the[0]+prop.the[1],-propa.dla*prop.gme);
 501                 wlos=false;
 502                 wscat=false;
 503
 504                 if (prop.wn<0.838 || prop.wn>210.0)
 505                         prop.kwx=mymax(prop.kwx,1);
 506
 507                 for (j=0; j<2; j++)
 508                         if (prop.hg[j]<1.0 || prop.hg[j]>1000.0)
 509                                 prop.kwx=mymax(prop.kwx,1);
 510
 511                 for (j=0; j<2; j++)
 512                         if (abs(prop.the[j]) >200e-3 || prop.dl[j]<0.1*propa.dls[j] || prop.dl[j]>3.0*propa.dls[j])
 513                                 prop.kwx=mymax(prop.kwx,3);
 514
 515                 if (prop.ens < 250.0 || prop.ens > 400.0 || prop.gme < 75e-9 || prop.gme > 250e-9 || prop_zgnd.real() <= abs(prop_zgnd.imag()) || prop.wn < 0.419 || prop.wn > 420.0)
 516
 517                         prop.kwx=4;
 518
 519                 for (j=0; j<2; j++)
 520
 521                 if (prop.hg[j]<0.5 || prop.hg[j]>3000.0)
 522                 {
 523                         prop.kwx=4;
 524                 }
 525
 526                 dmin=abs(prop.he[0]-prop.he[1])/200e-3;
 527                 q=adiff(0.0,prop,propa);
 528                 /* xae=pow(prop.wn*pow(prop.gme,2),-THIRD); */
 529                 xae=pow(prop.wn*prop.gme*prop.gme,-THIRD);
 530                 d3=mymax(propa.dlsa,1.3787*xae+propa.dla);
 531                 d4=d3+2.7574*xae;
 532                 a3=adiff(d3,prop,propa);
 533                 a4=adiff(d4,prop,propa);
 534                 propa.emd=(a4-a3)/(d4-d3);
 535                 propa.aed=a3-propa.emd*d3;
 536         }
 537
 538         if (prop.mdp>=0)
 539         {
 540                 prop.mdp=0;
 541                 prop.dist=d;
 542         }
 543
 544         if (prop.dist>0.0)
 545         {
 546                 if (prop.dist>1000e3)
 547                         prop.kwx=mymax(prop.kwx,1);
 548
 549                 if (prop.dist<dmin)
 550                         prop.kwx=mymax(prop.kwx,3);
 551
 552                 if (prop.dist<1e3 || prop.dist>2000e3)
 553                         prop.kwx=4;
 554         }
 555
 556         if (prop.dist<propa.dlsa)
 557         {
 558                 if (!wlos)
 559                 {
 560                         q=alos(0.0,prop,propa);
 561                         d2=propa.dlsa;
 562                         a2=propa.aed+d2*propa.emd;
 563                         d0=1.908*prop.wn*prop.he[0]*prop.he[1];
 564
 565                         if (propa.aed>=0.0)
 566                         {
 567                                 d0=mymin(d0,0.5*propa.dla);
 568                                 d1=d0+0.25*(propa.dla-d0);
 569                         }
 570
 571                         else
 572                                 d1=mymax(-propa.aed/propa.emd,0.25*propa.dla);
 573
 574                         a1=alos(d1,prop,propa);
 575                         wq=false;
 576
 577                         if (d0<d1)
 578                         {
 579                                 a0=alos(d0,prop,propa);
 580                                 q=log(d2/d0);
 581                                 propa.ak2=mymax(0.0,((d2-d0)*(a1-a0)-(d1-d0)*(a2-a0))/((d2-d0)*log(d1/d0)-(d1-d0)*q));
 582                                 wq=propa.aed>=0.0 || propa.ak2>0.0;
 583
 584                                 if (wq)
 585                                 {
 586                                         propa.ak1=(a2-a0-propa.ak2*q)/(d2-d0);
 587
 588                                         if (propa.ak1<0.0)
 589                                         {
 590                                                 propa.ak1=0.0;
 591                                                 propa.ak2=FORTRAN_DIM(a2,a0)/q;
 592
 593                                                 if (propa.ak2==0.0)
 594                                                         propa.ak1=propa.emd;
 595                                         }
 596                                 }
 597                         }
 598
 599                         if (!wq)
 600                         {
 601                                 propa.ak1=FORTRAN_DIM(a2,a1)/(d2-d1);
 602                                 propa.ak2=0.0;
 603
 604                                 if (propa.ak1==0.0)
 605                                         propa.ak1=propa.emd;
 606                         }
 607
 608                         propa.ael=a2-propa.ak1*d2-propa.ak2*log(d2);
 609                         wlos=true;
 610                 }
 611
 612                 if (prop.dist>0.0)
 613                         prop.aref=propa.ael+propa.ak1*prop.dist+propa.ak2*log(prop.dist);
 614         }
 615
 616         if (prop.dist<=0.0 || prop.dist>=propa.dlsa)
 617         {
 618                 if (!wscat)
 619                 {
 620                         q=ascat(0.0,prop,propa);
 621                         d5=propa.dla+200e3;
 622                         d6=d5+200e3;
 623                         a6=ascat(d6,prop,propa);
 624                         a5=ascat(d5,prop,propa);
 625
 626                         if (a5<1000.0)
 627                         {
 628                                 propa.ems=(a6-a5)/200e3;
 629                                 propa.dx=mymax(propa.dlsa,mymax(propa.dla+0.3*xae*log(47.7*prop.wn),(a5-propa.aed-propa.ems*d5)/(propa.emd-propa.ems)));
 630
 631                                 propa.aes=(propa.emd-propa.ems)*propa.dx+propa.aed;
 632                         }
 633
 634                         else
 635                         {
 636                                 propa.ems=propa.emd;
 637                                 propa.aes=propa.aed;
 638                                 propa.dx=10.e6;
 639                         }
 640
 641                         wscat=true;
 642                 }
 643
 644                 if (prop.dist>propa.dx)
 645                         prop.aref=propa.aes+propa.ems*prop.dist;
 646                 else
 647                         prop.aref=propa.aed+propa.emd*prop.dist;
 648         }
 649
 650         prop.aref=mymax(prop.aref,0.0);
 651 }
 652
 653 double curve (double const &c1, double const &c2, double const &x1, double const &x2, double const &x3, double const &de)
 654 {
 655         /* return (c1+c2/(1.0+pow((de-x2)/x3,2.0)))*pow(de/x1,2.0)/(1.0+pow(de/x1,2.0)); */
 656
 657         double temp1, temp2;
 658
 659         temp1=(de-x2)/x3;
 660         temp2=de/x1;
 661
 662         temp1*=temp1;
 663         temp2*=temp2;
 664
 665         return (c1+c2/(1.0+temp1))*temp2/(1.0+temp2);
 666 }
 667
 668 double avar(double zzt, double zzl, double zzc,prop_type &prop, propv_type &propv)
 669 {
 670         static int kdv;
 671         static double dexa, de, vmd, vs0, sgl, sgtm, sgtp, sgtd, tgtd, gm, gp, cv1, cv2, yv1, yv2, yv3, csm1, csm2, ysm1, ysm2, ysm3, csp1, csp2, ysp1, ysp2, ysp3, csd1, zd, cfm1, cfm2, cfm3, cfp1, cfp2, cfp3;
 672
 673         double bv1[7]={-9.67,-0.62,1.26,-9.21,-0.62,-0.39,3.15};
 674         double bv2[7]={12.7,9.19,15.5,9.05,9.19,2.86,857.9};
 675         double xv1[7]={144.9e3,228.9e3,262.6e3,84.1e3,228.9e3,141.7e3,2222.e3};
 676         double xv2[7]={190.3e3,205.2e3,185.2e3,101.1e3,205.2e3,315.9e3,164.8e3};
 677         double xv3[7]={133.8e3,143.6e3,99.8e3,98.6e3,143.6e3,167.4e3,116.3e3};
 678         double bsm1[7]={2.13,2.66,6.11,1.98,2.68,6.86,8.51};
 679         double bsm2[7]={159.5,7.67,6.65,13.11,7.16,10.38,169.8};
 680         double xsm1[7]={762.2e3,100.4e3,138.2e3,139.1e3,93.7e3,187.8e3,609.8e3};
 681         double xsm2[7]={123.6e3,172.5e3,242.2e3,132.7e3,186.8e3,169.6e3,119.9e3};
 682         double xsm3[7]={94.5e3,136.4e3,178.6e3,193.5e3,133.5e3,108.9e3,106.6e3};
 683         double bsp1[7]={2.11,6.87,10.08,3.68,4.75,8.58,8.43};
 684         double bsp2[7]={102.3,15.53,9.60,159.3,8.12,13.97,8.19};
 685         double xsp1[7]={636.9e3,138.7e3,165.3e3,464.4e3,93.2e3,216.0e3,136.2e3};
 686         double xsp2[7]={134.8e3,143.7e3,225.7e3,93.1e3,135.9e3,152.0e3,188.5e3};
 687         double xsp3[7]={95.6e3,98.6e3,129.7e3,94.2e3,113.4e3,122.7e3,122.9e3};
 688         double bsd1[7]={1.224,0.801,1.380,1.000,1.224,1.518,1.518};
 689         double bzd1[7]={1.282,2.161,1.282,20.,1.282,1.282,1.282};
 690         double bfm1[7]={1.0,1.0,1.0,1.0,0.92,1.0,1.0};
 691         double bfm2[7]={0.0,0.0,0.0,0.0,0.25,0.0,0.0};
 692         double bfm3[7]={0.0,0.0,0.0,0.0,1.77,0.0,0.0};
 693         double bfp1[7]={1.0,0.93,1.0,0.93,0.93,1.0,1.0};
 694         double bfp2[7]={0.0,0.31,0.0,0.19,0.31,0.0,0.0};
 695         double bfp3[7]={0.0,2.00,0.0,1.79,2.00,0.0,0.0};
 696         static bool ws, w1;
 697         double rt=7.8, rl=24.0, avarv, q, vs, zt, zl, zc;
 698         double sgt, yr;
 699         int temp_klim=propv.klim-1;
 700
 701         if (propv.lvar>0)
 702         {
 703                 switch (propv.lvar)
 704                 {
 705                         default:
 706                         if (propv.klim<=0 || propv.klim>7)
 707                         {
 708                                 propv.klim = 5;
 709                                 temp_klim = 4;
 710                                 prop.kwx=mymax(prop.kwx,2);
 711                         }
 712
 713                         cv1 = bv1[temp_klim];
 714                         cv2 = bv2[temp_klim];
 715                         yv1 = xv1[temp_klim];
 716                         yv2 = xv2[temp_klim];
 717                         yv3 = xv3[temp_klim];
 718                         csm1=bsm1[temp_klim];
 719                         csm2=bsm2[temp_klim];
 720                         ysm1=xsm1[temp_klim];
 721                         ysm2=xsm2[temp_klim];
 722                         ysm3=xsm3[temp_klim];
 723                         csp1=bsp1[temp_klim];
 724                         csp2=bsp2[temp_klim];
 725                         ysp1=xsp1[temp_klim];
 726                         ysp2=xsp2[temp_klim];
 727                         ysp3=xsp3[temp_klim];
 728                         csd1=bsd1[temp_klim];
 729                         zd=bzd1[temp_klim];
 730                         cfm1=bfm1[temp_klim];
 731                         cfm2=bfm2[temp_klim];
 732                         cfm3=bfm3[temp_klim];
 733                         cfp1=bfp1[temp_klim];
 734                         cfp2=bfp2[temp_klim];
 735                         cfp3=bfp3[temp_klim];
 736
 737                         case 4:
 738                         kdv=propv.mdvar;
 739                         ws=kdv>=20;
 740
 741                         if (ws)
 742                                 kdv-=20;
 743
 744                         w1=kdv>=10;
 745
 746                         if (w1)
 747                                 kdv-=10;
 748
 749                         if (kdv<0 || kdv>3)
 750                         {
 751                                 kdv=0;
 752                                 prop.kwx=mymax(prop.kwx,2);
 753                         }
 754
 755                         case 3:
 756                         q=log(0.133*prop.wn);
 757
 758                         /* gm=cfm1+cfm2/(pow(cfm3*q,2.0)+1.0); */
 759                         /* gp=cfp1+cfp2/(pow(cfp3*q,2.0)+1.0); */
 760
 761                         gm=cfm1+cfm2/((cfm3*q*cfm3*q)+1.0);
 762                         gp=cfp1+cfp2/((cfp3*q*cfp3*q)+1.0);
 763
 764                         case 2:
 765                         dexa=sqrt(18e6*prop.he[0])+sqrt(18e6*prop.he[1])+pow((575.7e12/prop.wn),THIRD);
 766
 767                         case 1:
 768                         if (prop.dist<dexa)
 769                                 de=130e3*prop.dist/dexa;
 770                         else
 771                                 de=130e3+prop.dist-dexa;
 772                 }
 773
 774                 vmd=curve(cv1,cv2,yv1,yv2,yv3,de);
 775                 sgtm=curve(csm1,csm2,ysm1,ysm2,ysm3,de)*gm;
 776                 sgtp=curve(csp1,csp2,ysp1,ysp2,ysp3,de)*gp;
 777                 sgtd=sgtp*csd1;
 778                 tgtd=(sgtp-sgtd)*zd;
 779
 780                 if (w1)
 781                         sgl=0.0;
 782                 else
 783                 {
 784                         q=(1.0-0.8*exp(-prop.dist/50e3))*prop.dh*prop.wn;
 785                         sgl=10.0*q/(q+13.0);
 786                 }
 787
 788                 if (ws)
 789                         vs0=0.0;
 790                 else
 791                 {
 792                         /* vs0=pow(5.0+3.0*exp(-de/100e3),2.0); */
 793                         vs0=(5.0+3.0*exp(-de/100e3));
 794                         vs0*=vs0;
 795                 }
 796
 797                 propv.lvar=0;
 798         }
 799
 800         zt=zzt;
 801         zl=zzl;
 802         zc=zzc;
 803
 804         switch (kdv)
 805         {
 806                 case 0:
 807                 zt=zc;
 808                 zl=zc;
 809                 break;
 810
 811                 case 1:
 812                 zl=zc;
 813                 break;
 814
 815                 case 2:
 816                 zl=zt;
 817         }
 818
 819         if (fabs(zt)>3.1 || fabs(zl)>3.1 || fabs(zc)>3.1)
 820                 prop.kwx=mymax(prop.kwx,1);
 821
 822         if (zt<0.0)
 823                 sgt=sgtm;
 824
 825         else if (zt<=zd)
 826                 sgt=sgtp;
 827
 828         else
 829                 sgt=sgtd+tgtd/zt;
 830
 831         /* vs=vs0+pow(sgt*zt,2.0)/(rt+zc*zc)+pow(sgl*zl,2.0)/(rl+zc*zc); */
 832         vs=vs0+(sgt*zt*sgt*zt)/(rt+zc*zc)+(sgl*zl*sgl*zl)/(rl+zc*zc);
 833
 834         if (kdv==0)
 835         {
 836                 yr=0.0;
 837                 propv.sgc=sqrt(sgt*sgt+sgl*sgl+vs);
 838         }
 839
 840         else if (kdv==1)
 841         {
 842                 yr=sgt*zt;
 843                 propv.sgc=sqrt(sgl*sgl+vs);
 844         }
 845
 846         else if (kdv==2)
 847         {
 848                 yr=sqrt(sgt*sgt+sgl*sgl)*zt;
 849                 propv.sgc=sqrt(vs);
 850         }
 851
 852         else
 853         {
 854                 yr=sgt*zt+sgl*zl;
 855                 propv.sgc=sqrt(vs);
 856         }
 857
 858         avarv=prop.aref-vmd-yr-propv.sgc*zc;
 859
 860         if (avarv<0.0)
 861                 avarv=avarv*(29.0-avarv)/(29.0-10.0*avarv);
 862
 863         return avarv;
 864 }
 865
 866 void hzns(double pfl[], prop_type &prop)
 867 {
 868         bool wq;
 869         int np;
 870         double xi, za, zb, qc, q, sb, sa;
 871
 872         np=(int)pfl[0];
 873         xi=pfl[1];
 874         za=pfl[2]+prop.hg[0];
 875         zb=pfl[np+2]+prop.hg[1];
 876         qc=0.5*prop.gme;
 877         q=qc*prop.dist;
 878         prop.the[1]=(zb-za)/prop.dist;
 879         prop.the[0]=prop.the[1]-q;
 880         prop.the[1]=-prop.the[1]-q;
 881         prop.dl[0]=prop.dist;
 882         prop.dl[1]=prop.dist;
 883
 884         if (np>=2)
 885         {
 886                 sa=0.0;
 887                 sb=prop.dist;
 888                 wq=true;
 889
 890                 for (int i=1; i<np; i++)
 891                 {
 892                         sa+=xi;
 893                         sb-=xi;
 894                         q=pfl[i+2]-(qc*sa+prop.the[0])*sa-za;
 895
 896                         if (q>0.0)
 897                         {
 898                                 prop.the[0]+=q/sa;
 899                                 prop.dl[0]=sa;
 900                                 wq=false;
 901                         }
 902
 903                         if (!wq)
 904                         {
 905                                 q=pfl[i+2]-(qc*sb+prop.the[1])*sb-zb;
 906
 907                                 if (q>0.0)
 908                                 {
 909                                         prop.the[1]+=q/sb;
 910                                         prop.dl[1]=sb;
 911                                 }
 912                         }
 913                 }
 914         }
 915 }
 916
 917 void z1sq1 (double z[], const double &x1, const double &x2, double& z0, double& zn)
 918
 919 {
 920         double xn, xa, xb, x, a, b;
 921         int n, ja, jb;
 922
 923         xn=z[0];
 924         xa=int(FORTRAN_DIM(x1/z[1],0.0));
 925         xb=xn-int(FORTRAN_DIM(xn,x2/z[1]));
 926
 927         if (xb<=xa)
 928         {
 929                 xa=FORTRAN_DIM(xa,1.0);
 930                 xb=xn-FORTRAN_DIM(xn,xb+1.0);
 931         }
 932
 933         ja=(int)xa;
 934         jb=(int)xb;
 935         n=jb-ja;
 936         xa=xb-xa;
 937         x=-0.5*xa;
 938         xb+=x;
 939         a=0.5*(z[ja+2]+z[jb+2]);
 940         b=0.5*(z[ja+2]-z[jb+2])*x;
 941
 942         for (int i=2; i<=n; ++i)
 943         {
 944                 ++ja;
 945                 x+=1.0;
 946                 a+=z[ja+2];
 947                 b+=z[ja+2]*x;
 948         }
 949
 950         a/=xa;
 951         b=b*12.0/((xa*xa+2.0)*xa);
 952         z0=a-b*xb;
 953         zn=a+b*(xn-xb);
 954 }
 955
 956 double qtile (const int &nn, double a[], const int &ir)
 957 {
 958         double q=0.0, r;                /* Initializations by KD2BD */
 959         int m, n, i, j, j1=0, i0=0, k;  /* Initializations by KD2BD */
 960         bool done=false;
 961         bool goto10=true;
 962
 963         m=0;
 964         n=nn;
 965         k=mymin(mymax(0,ir),n);
 966
 967         while (!done)
 968         {
 969                 if (goto10)
 970                 {
 971                         q=a[k];
 972                         i0=m;
 973                         j1=n;
 974                 }
 975
 976                 i=i0;
 977
 978                 while (i<=n && a[i]>=q)
 979                         i++;
 980
 981                 if (i>n)
 982                         i=n;
 983                 j=j1;
 984
 985                 while (j>=m && a[j]<=q)
 986                         j--;
 987
 988                 if (j<m)
 989                         j=m;
 990
 991                 if (i<j)
 992                 {
 993                         r=a[i];
 994                         a[i]=a[j];
 995                         a[j]=r;
 996                         i0=i+1;
 997                         j1=j-1;
 998                         goto10=false;
 999                 }
1000
1001                 else if (i<k)
1002                 {
1003                         a[k]=a[i];
1004                         a[i]=q;
1005                         m=i+1;
1006                         goto10=true;
1007                 }
1008
1009                 else if (j>k)
1010                 {
1011                         a[k]=a[j];
1012                         a[j]=q;
1013                         n=j-1;
1014                         goto10=true;
1015                 }
1016
1017                 else
1018                         done=true;
1019         }
1020
1021         return q;
1022 }
1023
1024 double qerf(const double &z)
1025 {
1026         double b1=0.319381530, b2=-0.356563782, b3=1.781477937;
1027         double b4=-1.821255987, b5=1.330274429;
1028         double rp=4.317008, rrt2pi=0.398942280;
1029         double t, x, qerfv;
1030
1031         x=z;
1032         t=fabs(x);
1033
1034         if (t>=10.0)
1035                 qerfv=0.0;
1036         else
1037         {
1038                 t=rp/(t+rp);
1039                 qerfv=exp(-0.5*x*x)*rrt2pi*((((b5*t+b4)*t+b3)*t+b2)*t+b1)*t;
1040         }
1041
1042         if (x<0.0)
1043                 qerfv=1.0-qerfv;
1044
1045         return qerfv;
1046 }
1047
1048 double d1thx(double pfl[], const double &x1, const double &x2)
1049 {
1050         int np, ka, kb, n, k, j;
1051         double d1thxv, sn, xa, xb;
1052         double *s;
1053
1054         np=(int)pfl[0];
1055         xa=x1/pfl[1];
1056         xb=x2/pfl[1];
1057         d1thxv=0.0;
1058
1059         if (xb-xa<2.0)  // exit out
1060                 return d1thxv;
1061
1062         ka=(int)(0.1*(xb-xa+8.0));
1063         ka=mymin(mymax(4,ka),25);
1064         n=10*ka-5;
1065         kb=n-ka+1;
1066         sn=n-1;
1067         assert((s=new double[n+2])!=0);
1068         s[0]=sn;
1069         s[1]=1.0;
1070         xb=(xb-xa)/sn;
1071         k=(int)(xa+1.0);
1072         xa-=(double)k;
1073
1074         for (j=0; j<n; j++)
1075         {
1076                 while (xa>0.0 && k<np)
1077                 {
1078                         xa-=1.0;
1079                         ++k;
1080                 }
1081
1082                 s[j+2]=pfl[k+2]+(pfl[k+2]-pfl[k+1])*xa;
1083                 xa=xa+xb;
1084         }
1085
1086         z1sq1(s,0.0,sn,xa,xb);
1087         xb=(xb-xa)/sn;
1088
1089         for (j=0; j<n; j++)
1090         {
1091                 s[j+2]-=xa;
1092                 xa=xa+xb;
1093         }
1094
1095         d1thxv=qtile(n-1,s+2,ka-1)-qtile(n-1,s+2,kb-1);
1096         d1thxv/=1.0-0.8*exp(-(x2-x1)/50.0e3);
1097         delete[] s;
1098
1099         return d1thxv;
1100 }
1101
1102 void qlrpfl (double pfl[], int klimx, int mdvarx, prop_type &prop, propa_type &propa, propv_type &propv)
1103 {
1104         int np, j;
1105         double xl[2], q, za, zb;
1106
1107         prop.dist=pfl[0]*pfl[1];
1108         np=(int)pfl[0];
1109         hzns(pfl,prop);
1110
1111         for (j=0; j<2; j++)
1112                 xl[j]=mymin(15.0*prop.hg[j],0.1*prop.dl[j]);
1113
1114         xl[1]=prop.dist-xl[1];
1115         prop.dh=d1thx(pfl,xl[0],xl[1]);
1116
1117         if (prop.dl[0]+prop.dl[1]>1.5*prop.dist)
1118         {
1119                 z1sq1(pfl,xl[0],xl[1],za,zb);
1120                 prop.he[0]=prop.hg[0]+FORTRAN_DIM(pfl[2],za);
1121                 prop.he[1]=prop.hg[1]+FORTRAN_DIM(pfl[np+2],zb);
1122
1123                 for (j=0; j<2; j++)
1124                         prop.dl[j]=sqrt(2.0*prop.he[j]/prop.gme)*exp(-0.07*sqrt(prop.dh/mymax(prop.he[j],5.0)));
1125
1126                 q=prop.dl[0]+prop.dl[1];
1127
1128                 if (q<=prop.dist)
1129                 {
1130                         /* q=pow(prop.dist/q,2.0); */
1131                         q=((prop.dist/q)*(prop.dist/q));
1132
1133                         for (j=0; j<2; j++)
1134                         {
1135                                 prop.he[j]*=q;
1136                                 prop.dl[j]=sqrt(2.0*prop.he[j]/prop.gme)*exp(-0.07*sqrt(prop.dh/mymax(prop.he[j],5.0)));
1137                         }
1138                 }
1139
1140                 for (j=0; j<2; j++)
1141                 {
1142                         q=sqrt(2.0*prop.he[j]/prop.gme);
1143                         prop.the[j]=(0.65*prop.dh*(q/prop.dl[j]-1.0)-2.0*prop.he[j])/q;
1144                 }
1145         }
1146
1147         else
1148         {
1149                 z1sq1(pfl,xl[0],0.9*prop.dl[0],za,q);
1150                 z1sq1(pfl,prop.dist-0.9*prop.dl[1],xl[1],q,zb);
1151                 prop.he[0]=prop.hg[0]+FORTRAN_DIM(pfl[2],za);
1152                 prop.he[1]=prop.hg[1]+FORTRAN_DIM(pfl[np+2],zb);
1153         }
1154
1155         prop.mdp=-1;
1156         propv.lvar=mymax(propv.lvar,3);
1157
1158         if (mdvarx>=0)
1159         {
1160                 propv.mdvar=mdvarx;
1161                 propv.lvar=mymax(propv.lvar,4);
1162         }
1163
1164         if (klimx>0)
1165         {
1166                 propv.klim=klimx;
1167                 propv.lvar=5;
1168         }
1169
1170         lrprop(0.0,prop,propa);
1171 }
1172
1173 double deg2rad(double d)
1174 {
1175         return d*3.1415926535897/180.0;
1176 }
1177
1178 //********************************************************
1179 //* Point-To-Point Mode Calculations                     *
1180 //********************************************************
1181
1182 void point_to_point(double elev[], double tht_m, double rht_m, double eps_dielect, double sgm_conductivity, double eno_ns_surfref, double frq_mhz, int radio_climate, int pol, double conf, double rel, double &dbloss, char *strmode, int &errnum)
1183 {
1184         // pol: 0-Horizontal, 1-Vertical
1185         // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1186         //                4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1187         //                7-Maritime Temperate, Over Sea
1188         // conf, rel: .01 to .99
1189         // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
1190         // errnum: 0- No Error.
1191         //         1- Warning: Some parameters are nearly out of range.
1192         //                     Results should be used with caution.
1193         //         2- Note: Default parameters have been substituted for impossible ones.
1194         //         3- Warning: A combination of parameters is out of range.
1195         //                     Results are probably invalid.
1196         //         Other-  Warning: Some parameters are out of range.
1197         //                          Results are probably invalid.
1198
1199         prop_type   prop;
1200         propv_type  propv;
1201         propa_type  propa;
1202
1203         double zsys=0;
1204         double zc, zr;
1205         double eno, enso, q;
1206         long ja, jb, i, np;
1207         double dkm, xkm;
1208         double fs;
1209
1210         prop.hg[0]=tht_m;
1211         prop.hg[1]=rht_m;
1212         propv.klim=radio_climate;
1213         prop.kwx=0;
1214         propv.lvar=5;
1215         prop.mdp=-1;
1216         zc=qerfi(conf);
1217         zr=qerfi(rel);
1218         np=(long)elev[0];
1219         dkm=(elev[1]*elev[0])/1000.0;
1220         xkm=elev[1]/1000.0;
1221         eno=eno_ns_surfref;
1222         enso=0.0;
1223         q=enso;
1224
1225         if (q<=0.0)
1226         {
1227                 /* ja = 3.0 + 0.1 * elev[0]; */
1228                 ja=(long)(3.0+0.1* elev[0]);
1229
1230                 jb=np-ja+6;
1231
1232                 for (i=ja-1; i<jb; ++i)
1233                         zsys+=elev[i];
1234
1235                 zsys/=(jb-ja+1);
1236                 q=eno;
1237         }
1238
1239         propv.mdvar=12;
1240         qlrps(frq_mhz,zsys,q,pol,eps_dielect,sgm_conductivity,prop);
1241         qlrpfl(elev,propv.klim,propv.mdvar,prop,propa,propv);
1242         fs=32.45+20.0*log10(frq_mhz)+20.0*log10(prop.dist/1000.0);
1243         q=prop.dist-propa.dla;
1244
1245         if (int(q)<0.0)
1246                 strcpy(strmode,"Line-Of-Sight Mode");
1247
1248         else
1249         {
1250                 if (int(q)==0.0)
1251                         strcpy(strmode,"Single Horizon");
1252
1253                 else if (int(q)>0.0)
1254                         strcpy(strmode,"Double Horizon");
1255
1256                 if (prop.dist<=propa.dlsa || prop.dist <= propa.dx)
1257                         strcat(strmode,", Diffraction Dominant");
1258
1259                 else if (prop.dist>propa.dx)
1260                         strcat(strmode, ", Troposcatter Dominant");
1261         }
1262
1263         dbloss=avar(zr,0.0,zc,prop,propv)+fs;
1264         errnum=prop.kwx;
1265 }
1266
1267 void point_to_pointMDH (double elev[], double tht_m, double rht_m, double eps_dielect, double sgm_conductivity, double eno_ns_surfref, double frq_mhz, int radio_climate, int pol, double timepct, double locpct, double confpct, double &dbloss, int &propmode, double &deltaH, int &errnum)
1268 {
1269         // pol: 0-Horizontal, 1-Vertical
1270         // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1271         //                4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1272         //                7-Maritime Temperate, Over Sea
1273         // timepct, locpct, confpct: .01 to .99
1274         // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
1275         // propmode:  Value   Mode
1276         //             -1     mode is undefined
1277         //              0     Line of Sight
1278         //              5     Single Horizon, Diffraction
1279         //              6     Single Horizon, Troposcatter
1280         //              9     Double Horizon, Diffraction
1281         //             10     Double Horizon, Troposcatter
1282         // errnum: 0- No Error.
1283         //         1- Warning: Some parameters are nearly out of range.
1284         //                     Results should be used with caution.
1285         //         2- Note: Default parameters have been substituted for impossible ones.
1286         //         3- Warning: A combination of parameters is out of range.
1287         //                     Results are probably invalid.
1288         //         Other-  Warning: Some parameters are out of range.
1289         //                          Results are probably invalid.
1290
1291         prop_type   prop;
1292         propv_type  propv;
1293         propa_type  propa;
1294         double zsys=0;
1295         double ztime, zloc, zconf;
1296         double eno, enso, q;
1297         long ja, jb, i, np;
1298         double dkm, xkm;
1299         double fs;
1300
1301         propmode=-1;  // mode is undefined
1302         prop.hg[0]=tht_m;
1303         prop.hg[1]=rht_m;
1304         propv.klim=radio_climate;
1305         prop.kwx=0;
1306         propv.lvar=5;
1307         prop.mdp=-1;
1308         ztime=qerfi(timepct);
1309         zloc=qerfi(locpct);
1310         zconf=qerfi(confpct);
1311
1312         np=(long)elev[0];
1313         dkm=(elev[1]*elev[0])/1000.0;
1314         xkm=elev[1]/1000.0;
1315         eno=eno_ns_surfref;
1316         enso=0.0;
1317         q=enso;
1318
1319         if (q<=0.0)
1320         {
1321                 /* ja = 3.0 + 0.1 * elev[0]; */
1322                 ja=(long)(3.0+0.1*elev[0]);
1323                 jb=np-ja+6;
1324
1325                 for (i=ja-1; i<jb; ++i)
1326                         zsys+=elev[i];
1327
1328                 zsys/=(jb-ja+1);
1329                 q=eno;
1330         }
1331
1332         propv.mdvar=12;
1333         qlrps(frq_mhz,zsys,q,pol,eps_dielect,sgm_conductivity,prop);
1334         qlrpfl(elev,propv.klim,propv.mdvar,prop,propa,propv);
1335         fs=32.45+20.0*log10(frq_mhz)+20.0*log10(prop.dist/1000.0);
1336         deltaH=prop.dh;
1337         q=prop.dist-propa.dla;
1338
1339         if (int(q)<0.0)
1340                 propmode=0;  // Line-Of-Sight Mode
1341         else
1342         {
1343                 if (int(q)==0.0)
1344                         propmode=4;  // Single Horizon
1345
1346                 else if (int(q)>0.0)
1347                         propmode=8;  // Double Horizon
1348
1349                 if (prop.dist<=propa.dlsa || prop.dist<=propa.dx)
1350                         propmode+=1; // Diffraction Dominant
1351
1352                 else if (prop.dist>propa.dx)
1353                         propmode+=2; // Troposcatter Dominant
1354         }
1355
1356         dbloss=avar(ztime, zloc, zconf, prop, propv)+fs;      //avar(time,location,confidence)
1357         errnum=prop.kwx;
1358 }
1359
1360 void point_to_pointDH (double elev[], double tht_m, double rht_m, double eps_dielect, double sgm_conductivity, double eno_ns_surfref, double frq_mhz, int radio_climate, int pol, double conf, double rel, double &dbloss, double &deltaH, int &errnum)
1361 {
1362         // pol: 0-Horizontal, 1-Vertical
1363         // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1364         //                4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1365         //                7-Maritime Temperate, Over Sea
1366         // conf, rel: .01 to .99
1367         // elev[]: [num points - 1], [delta dist(meters)], [height(meters) point 1], ..., [height(meters) point n]
1368         // errnum: 0- No Error.
1369         //         1- Warning: Some parameters are nearly out of range.
1370         //                     Results should be used with caution.
1371         //         2- Note: Default parameters have been substituted for impossible ones.
1372         //         3- Warning: A combination of parameters is out of range.
1373         //                     Results are probably invalid.
1374         //         Other-  Warning: Some parameters are out of range.
1375         //                          Results are probably invalid.
1376
1377         char strmode[100];
1378         prop_type   prop;
1379         propv_type  propv;
1380         propa_type  propa;
1381         double zsys=0;
1382         double zc, zr;
1383         double eno, enso, q;
1384         long ja, jb, i, np;
1385         double dkm, xkm;
1386         double fs;
1387
1388         prop.hg[0]=tht_m;
1389         prop.hg[1]=rht_m;
1390         propv.klim=radio_climate;
1391         prop.kwx=0;
1392         propv.lvar=5;
1393         prop.mdp=-1;
1394         zc=qerfi(conf);
1395         zr=qerfi(rel);
1396         np=(long)elev[0];
1397         dkm=(elev[1]*elev[0])/1000.0;
1398         xkm=elev[1]/1000.0;
1399         eno=eno_ns_surfref;
1400         enso=0.0;
1401         q=enso;
1402
1403         if (q<=0.0)
1404         {
1405                 /* ja = 3.0 + 0.1 * elev[0]; */
1406                 ja=(long)(3.0+0.1*elev[0]);
1407
1408                 jb=np-ja+6;
1409
1410                 for (i=ja-1; i<jb; ++i)
1411                         zsys+=elev[i];
1412
1413                 zsys/=(jb-ja+1);
1414                 q=eno;
1415         }
1416
1417         propv.mdvar=12;
1418         qlrps(frq_mhz,zsys,q,pol,eps_dielect,sgm_conductivity,prop);
1419         qlrpfl(elev,propv.klim,propv.mdvar,prop,propa,propv);
1420         fs=32.45+20.0*log10(frq_mhz)+20.0*log10(prop.dist/1000.0);
1421         deltaH=prop.dh;
1422         q=prop.dist-propa.dla;
1423
1424         if (int(q)<0.0)
1425                 strcpy(strmode,"Line-Of-Sight Mode");
1426         else
1427         {
1428                 if (int(q)==0.0)
1429                         strcpy(strmode,"Single Horizon");
1430
1431                 else if (int(q)>0.0)
1432                         strcpy(strmode,"Double Horizon");
1433
1434                 if (prop.dist<=propa.dlsa || prop.dist <= propa.dx)
1435                         strcat(strmode,", Diffraction Dominant");
1436
1437                 else if (prop.dist>propa.dx)
1438                         strcat(strmode, ", Troposcatter Dominant");
1439         }
1440
1441         dbloss=avar(zr,0.0,zc,prop,propv)+fs;      //avar(time,location,confidence)
1442         errnum=prop.kwx;
1443 }
1444
1445
1446 //********************************************************
1447 //* Area Mode Calculations                               *
1448 //********************************************************
1449
1450 void area(long ModVar, double deltaH, double tht_m, double rht_m, double dist_km, int TSiteCriteria, int RSiteCriteria, double eps_dielect, double sgm_conductivity, double eno_ns_surfref, double frq_mhz, int radio_climate, int pol, double pctTime, double pctLoc, double pctConf, double &dbloss, char *strmode, int &errnum)
1451 {
1452         // pol: 0-Horizontal, 1-Vertical
1453         // TSiteCriteria, RSiteCriteria:
1454         //                 0 - random, 1 - careful, 2 - very careful
1455         // radio_climate: 1-Equatorial, 2-Continental Subtropical, 3-Maritime Tropical,
1456         //                4-Desert, 5-Continental Temperate, 6-Maritime Temperate, Over Land,
1457         //                7-Maritime Temperate, Over Sea
1458         // ModVar: 0 - Single: pctConf is "Time/Situation/Location", pctTime, pctLoc not used
1459         //         1 - Individual: pctTime is "Situation/Location", pctConf is "Confidence", pctLoc not used
1460         //         2 - Mobile: pctTime is "Time/Locations (Reliability)", pctConf is "Confidence", pctLoc not used
1461         //         3 - Broadcast: pctTime is "Time", pctLoc is "Location", pctConf is "Confidence"
1462         // pctTime, pctLoc, pctConf: .01 to .99
1463         // errnum: 0- No Error.
1464         //         1- Warning: Some parameters are nearly out of range.
1465         //                     Results should be used with caution.
1466         //         2- Note: Default parameters have been substituted for impossible ones.
1467         //         3- Warning: A combination of parameters is out of range.
1468         //                     Results are probably invalid.
1469         //         Other-  Warning: Some parameters are out of range.
1470         //                          Results are probably invalid.
1471         // NOTE: strmode is not used at this time.
1472
1473         prop_type prop;
1474         propv_type propv;
1475         propa_type propa;
1476         double zt, zl, zc, xlb;
1477         double fs;
1478         long ivar;
1479         double eps, eno, sgm;
1480         long ipol;
1481         int kst[2];
1482
1483         kst[0]=(int)TSiteCriteria;
1484         kst[1]=(int)RSiteCriteria;
1485         zt=qerfi(pctTime);
1486         zl=qerfi(pctLoc);
1487         zc=qerfi(pctConf);
1488         eps=eps_dielect;
1489         sgm=sgm_conductivity;
1490         eno=eno_ns_surfref;
1491         prop.dh=deltaH;
1492         prop.hg[0]=tht_m;
1493         prop.hg[1]=rht_m;
1494         /* propv.klim = (__int32) radio_climate; */
1495         propv.klim=(long)radio_climate;
1496         prop.ens=eno;
1497         prop.kwx=0;
1498         ivar=(long)ModVar;
1499         ipol=(long)pol;
1500         qlrps(frq_mhz, 0.0, eno, ipol, eps, sgm, prop);
1501         qlra(kst, propv.klim, ivar, prop, propv);
1502
1503         if (propv.lvar<1)
1504                 propv.lvar=1;
1505
1506         lrprop(dist_km*1000.0, prop, propa);
1507         fs=32.45+20.0*log10(frq_mhz)+20.0*log10(prop.dist/1000.0);
1508         xlb=fs+avar(zt, zl, zc, prop, propv);
1509         dbloss=xlb;
1510
1511         if (prop.kwx==0)
1512                 errnum=0;
1513         else
1514                 errnum=prop.kwx;
1515 }
1516
1517 double ITMAreadBLoss(long ModVar, double deltaH, double tht_m, double rht_m, double dist_km, int TSiteCriteria, int RSiteCriteria, double eps_dielect, double sgm_conductivity, double eno_ns_surfref, double frq_mhz, int radio_climate, int pol, double pctTime, double pctLoc, double pctConf)
1518 {
1519         char strmode[200];
1520         int errnum;
1521         double dbloss;
1522
1523         area(ModVar,deltaH,tht_m,rht_m,dist_km,TSiteCriteria,RSiteCriteria, eps_dielect,sgm_conductivity,eno_ns_surfref, frq_mhz,radio_climate,pol,pctTime,pctLoc, pctConf,dbloss,strmode,errnum);
1524
1525         return dbloss;
1526 }
1527
1528 double ITMVersion()
1529 {
1530         return 7.0;
1531 }
1532