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

   1 /****************************************************************************
   2 *      SPLAT: An RF Signal Propagation Loss and Terrain Analysis Tool       *
   3 *                         Last update: 15-Mar-2007                          *
   4 *****************************************************************************
   5 *            Project started in 1997 by John A. Magliacane, KD2BD           *
   6 *****************************************************************************
   7 *                                                                           *
   8 *     Extensively modified by J. D. McDonald in Jan. 2004 to include        *
   9 *    the Longley-Rice propagation model using C++ code from NTIA/ITS.       *
  10 *                                                                           *
  11 *              See: http://flattop.its.bldrdoc.gov/itm.html                 *
  12 *                                                                           *
  13 *****************************************************************************
  14 *                                                                           *
  15 * This program is free software; you can redistribute it and/or modify it   *
  16 * under the terms of the GNU General Public License as published by the     *
  17 * Free Software Foundation; either version 2 of the License or any later    *
  18 * version.                                                                  *
  19 *                                                                           *
  20 * This program is distributed in the hope that it will useful, but WITHOUT  *
  21 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or     *
  22 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License     *
  23 * for more details.                                                         *
  24 *                                                                           *
  25 *****************************************************************************
  26  g++ -Wall -O3 -s -lm -lbz2 -fomit-frame-pointer itm.cpp splat.cpp -o splat
  27 *****************************************************************************/
  28
  29 #include <stdio.h>
  30 #include <math.h>
  31 #include <stdlib.h>
  32 #include <string.h>
  33 #include <ctype.h>
  34 #include <bzlib.h>
  35 #include <unistd.h>
  36 #include "fontdata.h"
  37 #include "smallfont.h"
  38
  39 #define GAMMA 2.5
  40 #define MAXSLOTS 9
  41 #define BZBUFFER 65536
  42
  43 #if MAXSLOTS==4
  44 #define ARRAYSIZE 4950
  45 #endif
  46
  47 #if MAXSLOTS==9
  48 #define ARRAYSIZE 10870
  49 #endif
  50
  51 #if MAXSLOTS==16
  52 #define ARRAYSIZE 19240
  53 #endif
  54
  55 #if MAXSLOTS==25
  56 #define ARRAYSIZE 30025
  57 #endif
  58
  59 char    string[255], sdf_path[255], opened=0, *splat_version={"1.2.0b"};
  60
  61 double  TWOPI=6.283185307179586, HALFPI=1.570796326794896,
  62         PI=3.141592653589793, deg2rad=1.74532925199e-02,
  63         EARTHRADIUS=20902230.97, METERS_PER_MILE=1609.344,
  64         METERS_PER_FOOT=0.3048, KM_PER_MILE=1.609344, earthradius,
  65         max_range=0.0;
  66
  67 int     min_north=90, max_north=-90, min_west=360, max_west=-1,
  68         max_elevation=-32768, min_elevation=32768, bzerror, maxdB=230;
  69
  70 unsigned char got_elevation_pattern=0, got_azimuth_pattern=0, metric=0;
  71
  72 struct site {   double lat;
  73                 double lon;
  74                 float alt;
  75                 char name[50];
  76             }   site;
  77
  78 struct path {   double lat[ARRAYSIZE];
  79                 double lon[ARRAYSIZE];
  80                 double elevation[ARRAYSIZE];
  81                 double distance[ARRAYSIZE];
  82                 int length;
  83             }   path;
  84
  85 struct dem {    int min_north;
  86                 int max_north;
  87                 int min_west;
  88                 int max_west;
  89                 int max_el;
  90                 int min_el;
  91                 short data[1200][1200];
  92                 unsigned char mask[1200][1200];
  93            }    dem[MAXSLOTS];
  94
  95 struct LR {     double eps_dielect;
  96                 double sgm_conductivity;
  97                 double eno_ns_surfref;
  98                 double frq_mhz;
  99                 double conf;
 100                 double rel;
 101                 int radio_climate;
 102                 int pol;
 103                 float antenna_pattern[361][1001];
 104           }     LR;
 105
 106 double elev_l[ARRAYSIZE+10];
 107
 108 void point_to_point(double elev[], double tht_m, double rht_m,
 109           double eps_dielect, double sgm_conductivity, double eno_ns_surfref,
 110           double frq_mhz, int radio_climate, int pol, double conf,
 111           double rel, double &dbloss, char *strmode, int &errnum);
 112
 113 double arccos(double x, double y)
 114 {
 115         /* This function implements the arc cosine function,
 116            returning a value between 0 and TWOPI. */
 117
 118         double result=0.0;
 119
 120         if (y>0.0)
 121                 result=acos(x/y);
 122
 123         if (y<0.0)
 124                 result=PI+acos(x/y);
 125
 126         return result;
 127 }
 128
 129 int ReduceAngle(double angle)
 130 {
 131         /* This function normalizes the argument to
 132            an integer angle between 0 and 180 degrees */
 133
 134         double temp;
 135
 136         temp=acos(cos(angle*deg2rad));
 137
 138         return (int)rint(temp/deg2rad);
 139 }
 140
 141 char *dec2dms(double decimal)
 142 {
 143         /* Converts decimal degrees to degrees, minutes, seconds,
 144            (DMS) and returns the result as a character string. */
 145
 146         char    sign;
 147         int     degrees, minutes, seconds;
 148         double  a, b, c, d;
 149
 150         if (decimal<0.0)
 151         {
 152                 decimal=-decimal;
 153                 sign=-1;
 154         }
 155
 156         else
 157                 sign=1;
 158
 159         a=floor(decimal);
 160         b=60.0*(decimal-a);
 161         c=floor(b);
 162         d=60.0*(b-c);
 163
 164         degrees=(int)a;
 165         minutes=(int)c;
 166         seconds=(int)d;
 167
 168         if (seconds<0)
 169                 seconds=0;
 170
 171         if (seconds>59)
 172                 seconds=59;
 173
 174         string[0]=0;
 175         sprintf(string,"%d%c %d\' %d\"", degrees*sign, 176, minutes, seconds);
 176         return (string);
 177 }
 178
 179 int OrMask(double lat, double lon, int value)
 180 {
 181         /* Lines, text, markings, and coverage areas are stored in a
 182            mask that is combined with topology data when topographic
 183            maps are generated by SPLAT!.  This function sets bits in
 184            the mask based on the latitude and longitude of the area
 185            pointed to. */
 186
 187         int     x, y, indx, minlat, minlon;
 188         char    found;
 189
 190         minlat=(int)floor(lat);
 191         minlon=(int)floor(lon);
 192
 193         for (indx=0, found=0; indx<MAXSLOTS && found==0;)
 194                 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
 195                         found=1;
 196                 else
 197                         indx++;
 198
 199         if (found)
 200         {
 201                 x=(int)(1199.0*(lat-floor(lat)));
 202                 y=(int)(1199.0*(lon-floor(lon)));
 203
 204                 dem[indx].mask[x][y]|=value;
 205
 206                 return (dem[indx].mask[x][y]);
 207         }
 208
 209         else
 210                 return -1;
 211 }
 212
 213 int GetMask(double lat, double lon)
 214 {
 215         /* This function returns the mask bits based on the latitude
 216            and longitude given. */
 217
 218         return (OrMask(lat,lon,0));
 219 }
 220
 221 double GetElevation(struct site location)
 222 {
 223         /* This function returns the elevation (in feet) of any location
 224            represented by the digital elevation model data in memory.
 225            Function returns -5000.0 for locations not found in memory. */
 226
 227         char    found;
 228         int     x, y, indx, minlat, minlon;
 229         double  elevation;
 230
 231         elevation=-5000.0;
 232
 233         minlat=(int)floor(location.lat);
 234         minlon=(int)floor(location.lon);
 235
 236         x=(int)(1199.0*(location.lat-floor(location.lat)));
 237         y=(int)(1199.0*(location.lon-floor(location.lon)));
 238
 239         for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
 240         {
 241                 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
 242                 {
 243
 244                         elevation=3.28084*dem[indx].data[x][y];
 245                         found=1;
 246                 }
 247         }
 248
 249         return elevation;
 250 }
 251
 252 int AddElevation(double lat, double lon, double height)
 253 {
 254         /* This function adds a user-defined terrain feature
 255            (in meters AGL) to the digital elevation model data
 256            in memory.  Does nothing and returns 0 for locations
 257            not found in memory. */
 258
 259         char    found;
 260         int     x, y, indx, minlat, minlon;
 261
 262         minlat=(int)floor(lat);
 263         minlon=(int)floor(lon);
 264
 265         x=(int)(1199.0*(lat-floor(lat)));
 266         y=(int)(1199.0*(lon-floor(lon)));
 267
 268         for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
 269         {
 270                 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
 271                 {
 272
 273                         dem[indx].data[x][y]+=(short)rint(height);
 274                         found=1;
 275                 }
 276         }
 277
 278         return found;
 279 }
 280
 281 double Distance(struct site site1, struct site site2)
 282 {
 283         /* This function returns the great circle distance
 284            in miles between any two site locations. */
 285
 286         double  lat1, lon1, lat2, lon2, distance;
 287
 288         lat1=site1.lat*deg2rad;
 289         lon1=site1.lon*deg2rad;
 290         lat2=site2.lat*deg2rad;
 291         lon2=site2.lon*deg2rad;
 292
 293         distance=3959.0*acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos((lon1)-(lon2)));
 294
 295         return distance;
 296 }
 297
 298 double Azimuth(struct site source, struct site destination)
 299 {
 300         /* This function returns the azimuth (in degrees) to the
 301            destination as seen from the location of the source. */
 302
 303         double  dest_lat, dest_lon, src_lat, src_lon,
 304                 beta, azimuth, diff, num, den, fraction;
 305
 306         dest_lat=destination.lat*deg2rad;
 307         dest_lon=destination.lon*deg2rad;
 308
 309         src_lat=source.lat*deg2rad;
 310         src_lon=source.lon*deg2rad;
 311
 312         /* Calculate Surface Distance */
 313
 314         beta=acos(sin(src_lat)*sin(dest_lat)+cos(src_lat)*cos(dest_lat)*cos(src_lon-dest_lon));
 315
 316         /* Calculate Azimuth */
 317
 318         num=sin(dest_lat)-(sin(src_lat)*cos(beta));
 319         den=cos(src_lat)*sin(beta);
 320         fraction=num/den;
 321
 322         /* Trap potential problems in acos() due to rounding */
 323
 324         if (fraction>=1.0)
 325                 fraction=1.0;
 326
 327         if (fraction<=-1.0)
 328                 fraction=-1.0;
 329
 330         /* Calculate azimuth */
 331
 332         azimuth=acos(fraction);
 333
 334         /* Reference it to True North */
 335
 336         diff=dest_lon-src_lon;
 337
 338         if (diff<=-PI)
 339                 diff+=TWOPI;
 340
 341         if (diff>=PI)
 342                 diff-=TWOPI;
 343
 344         if (diff>0.0)
 345                 azimuth=TWOPI-azimuth;
 346
 347         return (azimuth/deg2rad);
 348 }
 349
 350 double ElevationAngle(struct site source, struct site destination)
 351 {
 352         /* This function returns the angle of elevation (in degrees)
 353            of the destination as seen from the source location.
 354            A positive result represents an angle of elevation (uptilt),
 355            while a negative result represents an angle of depression
 356            (downtilt), as referenced to a normal to the center of
 357            the earth. */
 358
 359         register double a, b, dx;
 360
 361         a=GetElevation(destination)+destination.alt+earthradius;
 362         b=GetElevation(source)+source.alt+earthradius;
 363
 364         dx=5280.0*Distance(source,destination);
 365
 366         /* Apply the Law of Cosines */
 367
 368         return ((180.0*(acos(((b*b)+(dx*dx)-(a*a))/(2.0*b*dx)))/PI)-90.0);
 369 }
 370
 371 void ReadPath(struct site source, struct site destination)
 372 {
 373         /* This function generates a sequence of latitude and
 374            longitude positions between source and destination
 375            locations along a great circle path, and stores
 376            elevation and distance information for points
 377            along that path in the "path" structure. */
 378
 379         int     c;
 380         double  azimuth, distance, lat1, lon1, beta, den, num,
 381                 lat2, lon2, total_distance, x, y, path_length,
 382                 increment;
 383         struct  site tempsite;
 384
 385         lat1=source.lat*deg2rad;
 386         lon1=source.lon*deg2rad;
 387
 388         lat2=destination.lat*deg2rad;
 389         lon2=destination.lon*deg2rad;
 390
 391         azimuth=Azimuth(source,destination)*deg2rad;
 392
 393         total_distance=Distance(source,destination);
 394
 395         x=68755.0*acos(cos(lon1-lon2));         /* 1200 samples per degree */
 396         y=68755.0*acos(cos(lat1-lat2));         /* 68755 samples per radian */
 397
 398         path_length=sqrt((x*x)+(y*y));          /* Total number of samples */
 399
 400         increment=total_distance/path_length;   /* Miles per sample */
 401
 402         for (distance=0, c=0; distance<=total_distance; distance+=increment)
 403         {
 404                 beta=distance/3959.0;
 405                 lat2=asin(sin(lat1)*cos(beta)+cos(azimuth)*sin(beta)*cos(lat1));
 406                 num=cos(beta)-(sin(lat1)*sin(lat2));
 407                 den=cos(lat1)*cos(lat2);
 408
 409                 if (azimuth==0.0 && (beta>HALFPI-lat1))
 410                         lon2=lon1+PI;
 411
 412                 else if (azimuth==HALFPI && (beta>HALFPI+lat1))
 413                         lon2=lon1+PI;
 414
 415                 else if (fabs(num/den)>1.0)
 416                         lon2=lon1;
 417
 418                 else
 419                 {
 420                         if ((PI-azimuth)>=0.0)
 421                                 lon2=lon1-arccos(num,den);
 422                         else
 423                                 lon2=lon1+arccos(num,den);
 424                 }
 425
 426                 while (lon2<0.0)
 427                         lon2+=TWOPI;
 428
 429                 while (lon2>TWOPI)
 430                         lon2-=TWOPI;
 431
 432                 lat2=lat2/deg2rad;
 433                 lon2=lon2/deg2rad;
 434
 435                 if (c<ARRAYSIZE)
 436                 {
 437                         path.lat[c]=lat2;
 438                         path.lon[c]=lon2;
 439                         tempsite.lat=lat2;
 440                         tempsite.lon=lon2;
 441                         path.elevation[c]=GetElevation(tempsite);
 442                         path.distance[c]=distance;
 443                         c++;
 444                 }
 445         }
 446
 447         /* Make sure exact destination point is recorded at path.length-1 */
 448
 449         if (c<ARRAYSIZE)
 450         {
 451                 path.lat[c]=destination.lat;
 452                 path.lon[c]=destination.lon;
 453                 path.elevation[c]=GetElevation(destination);
 454                 path.distance[c]=total_distance;
 455                 c++;
 456         }
 457
 458         if (c<ARRAYSIZE)
 459                 path.length=c;
 460         else
 461                 path.length=ARRAYSIZE-1;
 462 }
 463
 464 double ElevationAngle2(struct site source, struct site destination, double er)
 465 {
 466         /* This function returns the angle of elevation (in degrees)
 467            of the destination as seen from the source location, UNLESS
 468            the path between the sites is obstructed, in which case, the
 469            elevation angle to the first obstruction is returned instead.
 470            "er" represents the earth radius. */
 471
 472         int     x;
 473         char    block=0;
 474         double  source_alt, destination_alt, cos_xmtr_angle,
 475                 cos_test_angle, test_alt, elevation, distance,
 476                 source_alt2, first_obstruction_angle=0.0;
 477         struct  path temp;
 478
 479         temp=path;
 480
 481         ReadPath(source,destination);
 482
 483         distance=5280.0*Distance(source,destination);
 484         source_alt=er+source.alt+GetElevation(source);
 485         destination_alt=er+destination.alt+GetElevation(destination);
 486         source_alt2=source_alt*source_alt;
 487
 488         /* Calculate the cosine of the elevation angle of the
 489            destination (receiver) as seen by the source (transmitter). */
 490
 491         cos_xmtr_angle=((source_alt2)+(distance*distance)-(destination_alt*destination_alt))/(2.0*source_alt*distance);
 492
 493         /* Test all points in between source and destination locations to
 494            see if the angle to a topographic feature generates a higher
 495            elevation angle than that produced by the destination.  Begin
 496            at the source since we're interested in identifying the FIRST
 497            obstruction along the path between source and destination. */
 498
 499         for (x=2, block=0; x<path.length && block==0; x++)
 500         {
 501                 distance=5280.0*path.distance[x];
 502
 503                 test_alt=earthradius+path.elevation[x];
 504
 505                 cos_test_angle=((source_alt2)+(distance*distance)-(test_alt*test_alt))/(2.0*source_alt*distance);
 506
 507                 /* Compare these two angles to determine if
 508                    an obstruction exists.  Since we're comparing
 509                    the cosines of these angles rather than
 510                    the angles themselves, the sense of the
 511                    following "if" statement is reversed from
 512                    what it would be if the angles themselves
 513                    were compared. */
 514
 515                 if (cos_xmtr_angle>cos_test_angle)
 516                 {
 517                         block=1;
 518                         first_obstruction_angle=((acos(cos_test_angle))/deg2rad)-90.0;
 519                 }
 520         }
 521
 522         if (block)
 523                 elevation=first_obstruction_angle;
 524
 525         else
 526                 elevation=((acos(cos_xmtr_angle))/deg2rad)-90.0;
 527
 528         path=temp;
 529
 530         return elevation;
 531 }
 532
 533 double AverageTerrain(struct site source, double azimuthx, double start_distance, double end_distance)
 534 {
 535         /* This function returns the average terrain calculated in
 536            the direction of "azimuth" (degrees) between "start_distance"
 537            and "end_distance" (miles) from the source location.  If
 538            the terrain is all water (non-critical error), -5000.0 is
 539            returned.  If not enough SDF data has been loaded into
 540            memory to complete the survey (critical error), then
 541            -9999.0 is returned. */
 542
 543         int     c, samples, endpoint;
 544         double  beta, lat1, lon1, lat2, lon2, num, den, azimuth, terrain=0.0;
 545         struct  site destination;
 546
 547         lat1=source.lat*deg2rad;
 548         lon1=source.lon*deg2rad;
 549
 550         /* Generate a path of elevations between the source
 551            location and the remote location provided. */
 552
 553         beta=end_distance/3959.0;
 554
 555         azimuth=deg2rad*azimuthx;
 556
 557         lat2=asin(sin(lat1)*cos(beta)+cos(azimuth)*sin(beta)*cos(lat1));
 558         num=cos(beta)-(sin(lat1)*sin(lat2));
 559         den=cos(lat1)*cos(lat2);
 560
 561         if (azimuth==0.0 && (beta>HALFPI-lat1))
 562                 lon2=lon1+PI;
 563
 564         else if (azimuth==HALFPI && (beta>HALFPI+lat1))
 565                 lon2=lon1+PI;
 566
 567         else if (fabs(num/den)>1.0)
 568                 lon2=lon1;
 569
 570         else
 571         {
 572                 if ((PI-azimuth)>=0.0)
 573                         lon2=lon1-arccos(num,den);
 574                 else
 575                         lon2=lon1+arccos(num,den);
 576         }
 577
 578         while (lon2<0.0)
 579                 lon2+=TWOPI;
 580
 581         while (lon2>TWOPI)
 582                 lon2-=TWOPI;
 583
 584         lat2=lat2/deg2rad;
 585         lon2=lon2/deg2rad;
 586
 587         destination.lat=lat2;
 588         destination.lon=lon2;
 589
 590         /* If SDF data is missing for the endpoint of
 591            the radial, then the average terrain cannot
 592            be accurately calculated.  Return -9999.0 */
 593
 594         if (GetElevation(destination)<-4999.0)
 595                 return (-9999.0);
 596         else
 597         {
 598                 ReadPath(source,destination);
 599
 600                 endpoint=path.length;
 601
 602                 /* Shrink the length of the radial if the
 603                    outermost portion is not over U.S. land. */
 604
 605                 for (c=endpoint-1; c>=0 && path.elevation[c]==0.0; c--);
 606
 607                 endpoint=c+1;
 608
 609                 for (c=0, samples=0; c<endpoint; c++)
 610                 {
 611                         if (path.distance[c]>=start_distance)
 612                         {
 613                                 terrain+=path.elevation[c];
 614                                 samples++;
 615                         }
 616                 }
 617
 618                 if (samples==0)
 619                         terrain=-5000.0;  /* No land */
 620                 else
 621                         terrain=(terrain/(double)samples);
 622
 623                 return terrain;
 624         }
 625 }
 626
 627 double haat(struct site antenna)
 628 {
 629         /* This function returns the antenna's Height Above Average
 630            Terrain (HAAT) based on FCC Part 73.313(d).  If a critical
 631            error occurs, such as a lack of SDF data to complete the
 632            survey, -5000.0 is returned. */
 633
 634         int     azi, c;
 635         char    error=0;
 636         double  terrain, avg_terrain, haat, sum=0.0;
 637
 638         /* Calculate the average terrain between 2 and 10 miles
 639            from the antenna site at azimuths of 0, 45, 90, 135,
 640            180, 225, 270, and 315 degrees. */
 641
 642         for (c=0, azi=0; azi<=315 && error==0; azi+=45)
 643         {
 644                 terrain=AverageTerrain(antenna, (double)azi, 2.0, 10.0);
 645
 646                 if (terrain<-9998.0)  /* SDF data is missing */
 647                         error=1;
 648
 649                 if (terrain>-4999.0)  /* It's land, not water */
 650                 {
 651                         sum+=terrain;  /* Sum of averages */
 652                         c++;
 653                 }
 654         }
 655
 656         if (error)
 657                 return -5000.0;
 658         else
 659         {
 660                 avg_terrain=(sum/(double)c);
 661                 haat=(antenna.alt+GetElevation(antenna))-avg_terrain;
 662                 return haat;
 663         }
 664 }
 665
 666 float LonDiff(float lon1, float lon2)
 667 {
 668         /* This function returns the short path longitudinal
 669            difference between longitude1 and longitude2
 670            as an angle between -180.0 and +180.0 degrees.
 671            If lon1 is west of lon2, the result is positive.
 672            If lon1 is east of lon2, the result is negative. */
 673
 674         float diff;
 675
 676         diff=lon1-lon2;
 677
 678         if (diff<=-180.0)
 679                 diff+=360.0;
 680
 681         if (diff>=180.0)
 682                 diff-=360.0;
 683
 684         return diff;
 685 }
 686
 687 void PlaceMarker(struct site location)
 688 {
 689         /* This function places text and marker data in the mask array
 690            for illustration on topographic maps generated by SPLAT!.
 691            By default, SPLAT! centers text information BELOW the marker,
 692            but may move it above, to the left, or to the right of the
 693            marker depending on how much room is available on the map,
 694            or depending on whether the area is already occupied by
 695            another marker or label.  If no room or clear space is
 696            available on the map to place the marker and its associated
 697            text, then the marker and text are not written to the map. */
 698
 699         int     a, b, c, byte;
 700         char    ok2print, occupied;
 701         double  x, y, lat, lon, textx=0.0, texty=0.0, xmin, xmax,
 702                 ymin, ymax, p1, p3, p6, p8, p12, p16, p24, label_length;
 703
 704         xmin=min_north;
 705         xmax=max_north;
 706         ymin=min_west;
 707         ymax=max_west;
 708         lat=location.lat;
 709         lon=location.lon;
 710
 711         if (lat<xmax && lat>xmin && (LonDiff(lon,ymax)<0.0) && (LonDiff(lon,ymin)>0.0))
 712         {
 713                 p1=1.0/1200.0;
 714                 p3=3.0/1200.0;
 715                 p6=6.0/1200.0;
 716                 p8=8.0/1200.0;
 717                 p12=12.0/1200.0;
 718                 p16=16.0/1200.0;
 719                 p24=24.0/1200.0;
 720                 ok2print=0;
 721                 occupied=0;
 722
 723                 /* Is Marker Position Clear Of Text Or Other Markers? */
 724
 725                 for (x=lat-p3; (x<=xmax && x>=xmin && x<=lat+p3); x+=p1)
 726                         for (y=lon-p3; (LonDiff(y,ymax)<=0.0) && (LonDiff(y,ymin)>=0.0) && (LonDiff(y,lon+p3)<=0.0); y+=p1)
 727                                 occupied|=(GetMask(x,y)&2);
 728
 729                 if (occupied==0)
 730                 {
 731                         /* Determine Where Text Can Be Positioned */
 732
 733                         /* label_length=length in pixels.
 734                            Each character is 8 pixels wide. */
 735
 736                         label_length=p1*(double)(strlen(location.name)<<3);
 737
 738                         if ((LonDiff(lon+label_length,ymax)<=0.0) && (LonDiff(lon-label_length,ymin)>=0.0))
 739                         {
 740                                 /* Default: Centered Text */
 741
 742                                 texty=lon+label_length/2.0;
 743
 744                                 if ((lat-p8)>=p16)
 745                                 {
 746                                         /* Position Text Below The Marker */
 747
 748                                         textx=lat-p8;
 749
 750                                         x=textx;
 751                                         y=texty;
 752
 753                                         /* Is This Position Clear Of
 754                                            Text Or Other Markers? */
 755
 756                                         for (a=0, occupied=0; a<16; a++)
 757                                         {
 758                                                 for (b=0; b<(int)strlen(location.name); b++)
 759                                                         for (c=0; c<8; c++, y-=p1)
 760                                                                 occupied|=(GetMask(x,y)&2);
 761                                                 x-=p1;
 762                                                 y=texty;
 763                                         }
 764
 765                                         x=textx;
 766                                         y=texty;
 767
 768                                         if (occupied==0)
 769                                                 ok2print=1;
 770                                 }
 771
 772                                 else
 773                                 {
 774                                         /* Position Text Above The Marker */
 775
 776                                         textx=lat+p24;
 777
 778                                         x=textx;
 779                                         y=texty;
 780
 781                                         /* Is This Position Clear Of
 782                                            Text Or Other Markers? */
 783
 784                                         for (a=0, occupied=0; a<16; a++)
 785                                         {
 786                                                 for (b=0; b<(int)strlen(location.name); b++)
 787                                                         for (c=0; c<8; c++, y-=p1)
 788                                                                 occupied|=(GetMask(x,y)&2);
 789                                                 x-=p1;
 790                                                 y=texty;
 791                                         }
 792
 793                                         x=textx;
 794                                         y=texty;
 795
 796                                         if (occupied==0)
 797                                                 ok2print=1;
 798                                 }
 799                         }
 800
 801                         if (ok2print==0)
 802                         {
 803                                 if (LonDiff(lon-label_length,ymin)>=0.0)
 804                                 {
 805                                         /* Position Text To The
 806                                            Right Of The Marker */
 807
 808                                         textx=lat+p6;
 809                                         texty=lon-p12;
 810
 811                                         x=textx;
 812                                         y=texty;
 813
 814                                         /* Is This Position Clear Of
 815                                            Text Or Other Markers? */
 816
 817                                         for (a=0, occupied=0; a<16; a++)
 818                                         {
 819                                                 for (b=0; b<(int)strlen(location.name); b++)
 820                                                         for (c=0; c<8; c++, y-=p1)
 821                                                                 occupied|=(GetMask(x,y)&2);
 822                                                 x-=p1;
 823                                                 y=texty;
 824                                         }
 825
 826                                         x=textx;
 827                                         y=texty;
 828
 829                                         if (occupied==0)
 830                                                 ok2print=1;
 831                                 }
 832
 833                                 else
 834                                 {
 835                                         /* Position Text To The
 836                                            Left Of The Marker */
 837
 838                                         textx=lat+p6;
 839                                         texty=lon+p8+(label_length);
 840
 841                                         x=textx;
 842                                         y=texty;
 843
 844                                         /* Is This Position Clear Of
 845                                            Text Or Other Markers? */
 846
 847                                         for (a=0, occupied=0; a<16; a++)
 848                                         {
 849                                                 for (b=0; b<(int)strlen(location.name); b++)
 850                                                         for (c=0; c<8; c++, y-=p1)
 851                                                                 occupied|=(GetMask(x,y)&2);
 852                                                 x-=p1;
 853                                                 y=texty;
 854                                         }
 855
 856                                         x=textx;
 857                                         y=texty;
 858
 859                                         if (occupied==0)
 860                                                 ok2print=1;
 861                                 }
 862                         }
 863
 864                         /* textx and texty contain the latitude and longitude
 865                            coordinates that describe the placement of the text
 866                            on the map. */
 867
 868                         if (ok2print)
 869                         {
 870                                 /* Draw Text */
 871
 872                                 x=textx;
 873                                 y=texty;
 874
 875                                 for (a=0; a<16 && ok2print; a++)
 876                                 {
 877                                         for (b=0; b<(int)strlen(location.name); b++)
 878                                         {
 879                                                 byte=fontdata[16*(location.name[b])+a];
 880
 881                                                 for (c=128; c>0; c=c>>1, y-=p1)
 882                                                         if (byte&c)
 883                                                                 OrMask(x,y,2);
 884                                         }
 885
 886                                         x-=p1;
 887                                         y=texty;
 888                                 }
 889
 890                                 /* Draw Square Marker Centered
 891                                    On Location Specified */
 892
 893                                 for (x=lat-p3; (x<=xmax && x>=xmin && x<=lat+p3); x+=p1)
 894                                         for (y=lon-p3; (LonDiff(y,ymax)<=0.0) && (LonDiff(y,ymin)>=0.0) && (LonDiff(y,lon+p3)<=0.0); y+=p1)
 895                                                 OrMask(x,y,2);
 896                         }
 897                 }
 898         }
 899 }
 900
 901 double ReadBearing(char *input)
 902 {
 903         /* This function takes numeric input in the form of a character
 904            string, and returns an equivalent bearing in degrees as a
 905            decimal number (double).  The input may either be expressed
 906            in decimal format (40.139722) or degree, minute, second
 907            format (40 08 23).  This function also safely handles
 908            extra spaces found either leading, trailing, or
 909            embedded within the numbers expressed in the
 910            input string.  Decimal seconds are permitted. */
 911
 912         double  seconds, bearing=0.0;
 913         char    string[20];
 914         int     a, b, length, degrees, minutes;
 915
 916         /* Copy "input" to "string", and ignore any extra
 917            spaces that might be present in the process. */
 918
 919         string[0]=0;
 920         length=strlen(input);
 921
 922         for (a=0, b=0; a<length && a<18; a++)
 923         {
 924                 if ((input[a]!=32 && input[a]!='\n') || (input[a]==32 && input[a+1]!=32 && input[a+1]!='\n' && b!=0))
 925                 {
 926                         string[b]=input[a];
 927                         b++;
 928                 }
 929         }
 930
 931         string[b]=0;
 932
 933         /* Count number of spaces in the clean string. */
 934
 935         length=strlen(string);
 936
 937         for (a=0, b=0; a<length; a++)
 938                 if (string[a]==32)
 939                         b++;
 940
 941         if (b==0)  /* Decimal Format (40.139722) */
 942                 sscanf(string,"%lf",&bearing);
 943
 944         if (b==2)  /* Degree, Minute, Second Format (40 08 23) */
 945         {
 946                 sscanf(string,"%d %d %lf",&degrees, &minutes, &seconds);
 947
 948                 bearing=(double)abs(degrees)+((double)abs(minutes)/60)+(fabs(seconds)/3600);
 949
 950                 if ((degrees<0) || (minutes<0) || (seconds<0.0))
 951                         bearing=-bearing;
 952         }
 953
 954         /* Anything else returns a 0.0 */
 955
 956         if (bearing>360.0 || bearing<-90.0)
 957                 bearing=0.0;
 958
 959         return bearing;
 960 }
 961
 962 struct site LoadQTH(char *filename)
 963 {
 964         /* This function reads SPLAT! .qth (site location) files.
 965            The latitude and longitude may be expressed either in
 966            decimal degrees, or in degree, minute, second format.
 967            Antenna height is assumed to be expressed in feet above
 968            ground level (AGL), unless followed by the letter 'M',
 969            or 'm', or by the word "meters" or "Meters", in which
 970            case meters is assumed, and is handled accordingly. */
 971
 972         int     x;
 973         char    string[50], qthfile[255];
 974         struct  site tempsite;
 975         FILE    *fd=NULL;
 976
 977         for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
 978                 qthfile[x]=filename[x];
 979
 980         qthfile[x]='.';
 981         qthfile[x+1]='q';
 982         qthfile[x+2]='t';
 983         qthfile[x+3]='h';
 984         qthfile[x+4]=0;
 985
 986         tempsite.lat=91.0;
 987         tempsite.lon=361.0;
 988         tempsite.alt=0.0;
 989         tempsite.name[0]=0;
 990
 991         fd=fopen(qthfile,"r");
 992
 993         if (fd!=NULL)
 994         {
 995                 /* Site Name */
 996                 fgets(string,49,fd);
 997
 998                 /* Strip <CR> and/or <LF> from end of site name */
 999
1000                 for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0; tempsite.name[x]=string[x], x++);
1001
1002                 tempsite.name[x]=0;
1003
1004                 /* Site Latitude */
1005                 fgets(string,49,fd);
1006                 tempsite.lat=ReadBearing(string);
1007
1008                 /* Site Longitude */
1009                 fgets(string,49,fd);
1010                 tempsite.lon=ReadBearing(string);
1011
1012                 /* Antenna Height */
1013                 fgets(string,49,fd);
1014                 fclose(fd);
1015
1016                 /* Remove <CR> and/or <LF> from antenna height string */
1017                 for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0; x++);
1018
1019                 string[x]=0;
1020
1021                 /* Antenna height may either be in feet or meters.
1022                    If the letter 'M' or 'm' is discovered in
1023                    the string, then this is an indication that
1024                    the value given is expressed in meters, and
1025                    must be converted to feet before exiting. */
1026
1027                 for (x=0; string[x]!='M' && string[x]!='m' && string[x]!=0 && x<48; x++);
1028                 if (string[x]=='M' || string[x]=='m')
1029                 {
1030                         string[x]=0;
1031                         sscanf(string,"%f",&tempsite.alt);
1032                         tempsite.alt*=3.28084;
1033                 }
1034
1035                 else
1036                 {
1037                         string[x]=0;
1038                         sscanf(string,"%f",&tempsite.alt);
1039                 }
1040         }
1041
1042         return tempsite;
1043 }
1044
1045 void LoadPAT(char *filename)
1046 {
1047         /* This function reads and processes antenna pattern (.az
1048            and .el) files that correspond in name to previously
1049            loaded SPLAT! .lrp files.  */
1050
1051         int     a, b, w, x, y, z, last_index, next_index, span;
1052         char    string[255], azfile[255], elfile[255], *pointer=NULL;
1053         float   az, xx, elevation, amplitude, rotation, valid1, valid2,
1054                 delta, azimuth[361], azimuth_pattern[361], el_pattern[10001],
1055                 elevation_pattern[361][1001], slant_angle[361], tilt,
1056                 mechanical_tilt, tilt_azimuth, tilt_increment, sum;
1057         FILE    *fd=NULL;
1058         unsigned char read_count[10001];
1059
1060         for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
1061         {
1062                 azfile[x]=filename[x];
1063                 elfile[x]=filename[x];
1064         }
1065
1066         azfile[x]='.';
1067         azfile[x+1]='a';
1068         azfile[x+2]='z';
1069         azfile[x+3]=0;
1070
1071         elfile[x]='.';
1072         elfile[x+1]='e';
1073         elfile[x+2]='l';
1074         elfile[x+3]=0;
1075
1076         rotation=0.0;
1077
1078         /* Load .az antenna pattern file */
1079
1080         fd=fopen(azfile,"r");
1081
1082         if (fd!=NULL)
1083         {
1084                 /* Clear azimuth pattern array */
1085
1086                 for (x=0; x<=360; x++)
1087                 {
1088                         azimuth[x]=0.0;
1089                         read_count[x]=0;
1090                 }
1091
1092
1093                 /* Read azimuth pattern rotation
1094                    in degrees measured clockwise
1095                    from true North. */
1096
1097                 fgets(string,254,fd);
1098                 pointer=strchr(string,';');
1099
1100                 if (pointer!=NULL)
1101                         *pointer=0;
1102
1103                 sscanf(string,"%f",&rotation);
1104
1105
1106                 /* Read azimuth (degrees) and corresponding
1107                    normalized field radiation pattern amplitude
1108                    (0.0 to 1.0) until EOF is reached. */
1109
1110                 fgets(string,254,fd);
1111                 pointer=strchr(string,';');
1112
1113                 if (pointer!=NULL)
1114                         *pointer=0;
1115
1116                 sscanf(string,"%f %f",&az, &amplitude);
1117
1118                 do
1119                 {
1120                         x=(int)rintf(az);
1121
1122                         if (x>=0 && x<=360 && fd!=NULL)
1123                         {
1124                                 azimuth[x]+=amplitude;
1125                                 read_count[x]++;
1126                         }
1127
1128                         fgets(string,254,fd);
1129                         pointer=strchr(string,';');
1130
1131                         if (pointer!=NULL)
1132                                 *pointer=0;
1133
1134                         sscanf(string,"%f %f",&az, &amplitude);
1135
1136                 } while (feof(fd)==0);
1137
1138                 fclose(fd);
1139
1140
1141                 /* Handle 0=360 degree ambiguity */
1142
1143                 if ((read_count[0]==0) && (read_count[360]!=0))
1144                 {
1145                         read_count[0]=read_count[360];
1146                         azimuth[0]=azimuth[360];
1147                 }
1148
1149                 if ((read_count[0]!=0) && (read_count[360]==0))
1150                 {
1151                         read_count[360]=read_count[0];
1152                         azimuth[360]=azimuth[0];
1153                 }
1154
1155                 /* Average pattern values in case more than
1156                     one was read for each degree of azimuth. */
1157
1158                 for (x=0; x<=360; x++)
1159                 {
1160                         if (read_count[x]>1)
1161                                 azimuth[x]/=(float)read_count[x];
1162                 }
1163
1164                 /* Interpolate missing azimuths
1165                    to completely fill the array */
1166
1167                 last_index=-1;
1168                 next_index=-1;
1169
1170                 for (x=0; x<=360; x++)
1171                 {
1172                         if (read_count[x]!=0)
1173                         {
1174                                 if (last_index==-1)
1175                                         last_index=x;
1176                                 else
1177                                         next_index=x;
1178                         }
1179
1180                         if (last_index!=-1 && next_index!=-1)
1181                         {
1182                                 valid1=azimuth[last_index];
1183                                 valid2=azimuth[next_index];
1184
1185                                 span=next_index-last_index;
1186                                 delta=(valid2-valid1)/(float)span;
1187
1188                                 for (y=last_index+1; y<next_index; y++)
1189                                         azimuth[y]=azimuth[y-1]+delta;
1190
1191                                 last_index=y;
1192                                 next_index=-1;
1193                         }
1194                 }
1195
1196                 /* Perform azimuth pattern rotation
1197                    and load azimuth_pattern[361] with
1198                    azimuth pattern data in its final form. */
1199
1200                 for (x=0; x<360; x++)
1201                 {
1202                         y=x+(int)rintf(rotation);
1203
1204                         if (y>=360)
1205                                 y-=360;
1206
1207                         azimuth_pattern[y]=azimuth[x];
1208                 }
1209
1210                 azimuth_pattern[360]=azimuth_pattern[0];
1211
1212                 got_azimuth_pattern=255;
1213         }
1214
1215         /* Read and process .el file */
1216
1217         fd=fopen(elfile,"r");
1218
1219         if (fd!=NULL)
1220         {
1221                 for (x=0; x<=10000; x++)
1222                 {
1223                         el_pattern[x]=0.0;
1224                         read_count[x]=0;
1225                 }
1226
1227                 /* Read mechanical tilt (degrees) and
1228                    tilt azimuth in degrees measured
1229                    clockwise from true North. */
1230
1231                 fgets(string,254,fd);
1232                 pointer=strchr(string,';');
1233
1234                 if (pointer!=NULL)
1235                         *pointer=0;
1236
1237                 sscanf(string,"%f %f",&mechanical_tilt, &tilt_azimuth);
1238
1239                 /* Read elevation (degrees) and corresponding
1240                    normalized field radiation pattern amplitude
1241                    (0.0 to 1.0) until EOF is reached. */
1242
1243                 fgets(string,254,fd);
1244                 pointer=strchr(string,';');
1245
1246                 if (pointer!=NULL)
1247                         *pointer=0;
1248
1249                 sscanf(string,"%f %f", &elevation, &amplitude);
1250
1251                 while (feof(fd)==0)
1252                 {
1253                         /* Read in normalized radiated field values
1254                            for every 0.01 degrees of elevation between
1255                            -10.0 and +90.0 degrees */
1256
1257                         x=(int)rintf(100.0*(elevation+10.0));
1258
1259                         if (x>=0 && x<=10000)
1260                         {
1261                                 el_pattern[x]+=amplitude;
1262                                 read_count[x]++;
1263                         }
1264
1265                         fgets(string,254,fd);
1266                         pointer=strchr(string,';');
1267
1268                         if (pointer!=NULL)
1269                                 *pointer=0;
1270
1271                         sscanf(string,"%f %f", &elevation, &amplitude);
1272                 }
1273
1274                 fclose(fd);
1275
1276                 /* Average the field values in case more than
1277                    one was read for each 0.01 degrees of elevation. */
1278
1279                 for (x=0; x<=10000; x++)
1280                 {
1281                         if (read_count[x]>1)
1282                                 el_pattern[x]/=(float)read_count[x];
1283                 }
1284
1285                 /* Interpolate between missing elevations (if
1286                    any) to completely fill the array and provide
1287                    radiated field values for every 0.01 degrees of
1288                    elevation. */
1289
1290                 last_index=-1;
1291                 next_index=-1;
1292
1293                 for (x=0; x<=10000; x++)
1294                 {
1295                         if (read_count[x]!=0)
1296                         {
1297                                 if (last_index==-1)
1298                                         last_index=x;
1299                                 else
1300                                         next_index=x;
1301                         }
1302
1303                         if (last_index!=-1 && next_index!=-1)
1304                         {
1305                                 valid1=el_pattern[last_index];
1306                                 valid2=el_pattern[next_index];
1307
1308                                 span=next_index-last_index;
1309                                 delta=(valid2-valid1)/(float)span;
1310
1311                                 for (y=last_index+1; y<next_index; y++)
1312                                         el_pattern[y]=el_pattern[y-1]+delta;
1313
1314                                 last_index=y;
1315                                 next_index=-1;
1316                         }
1317                 }
1318
1319                 /* Fill slant_angle[] array with offset angles based
1320                    on the antenna's mechanical beam tilt (if any)
1321                    and tilt direction (azimuth). */
1322
1323                 if (mechanical_tilt==0.0)
1324                 {
1325                         for (x=0; x<=360; x++)
1326                                 slant_angle[x]=0.0;
1327                 }
1328
1329                 else
1330                 {
1331                         tilt_increment=mechanical_tilt/90.0;
1332
1333                         for (x=0; x<=360; x++)
1334                         {
1335                                 xx=(float)x;
1336                                 y=(int)rintf(tilt_azimuth+xx);
1337
1338                                 while (y>=360)
1339                                         y-=360;
1340
1341                                 while (y<0)
1342                                         y+=360;
1343
1344                                 if (x<=180)
1345                                         slant_angle[y]=-(tilt_increment*(90.0-xx));
1346
1347                                 if (x>180)
1348                                         slant_angle[y]=-(tilt_increment*(xx-270.0));
1349                         }
1350                 }
1351
1352                 slant_angle[360]=slant_angle[0];   /* 360 degree wrap-around */
1353
1354                 for (w=0; w<=360; w++)
1355                 {
1356                         tilt=slant_angle[w];
1357
1358                         /** Convert tilt angle to
1359                             an array index offset **/
1360
1361                         y=(int)rintf(100.0*tilt);
1362
1363                         /* Copy shifted el_pattern[10001] field
1364                            values into elevation_pattern[361][1001]
1365                            at the corresponding azimuth, downsampling
1366                            (averaging) along the way in chunks of 10. */
1367
1368                         for (x=y, z=0; z<=1000; x+=10, z++)
1369                         {
1370                                 for (sum=0.0, a=0; a<10; a++)
1371                                 {
1372                                         b=a+x;
1373
1374                                         if (b>=0 && b<=10000)
1375                                                 sum+=el_pattern[b];
1376                                         if (b<0)
1377                                                 sum+=el_pattern[0];
1378                                         if (b>10000)
1379                                                 sum+=el_pattern[10000];
1380                                 }
1381
1382                                 elevation_pattern[w][z]=sum/10.0;
1383                         }
1384                 }
1385
1386                 got_elevation_pattern=255;
1387         }
1388
1389         for (x=0; x<=360; x++)
1390         {
1391                 for (y=0; y<=1000; y++)
1392                 {
1393                         if (got_elevation_pattern)
1394                                 elevation=elevation_pattern[x][y];
1395                         else
1396                                 elevation=1.0;
1397
1398                         if (got_azimuth_pattern)
1399                                 az=azimuth_pattern[x];
1400                         else
1401                                 az=1.0;
1402
1403                         LR.antenna_pattern[x][y]=az*elevation;
1404                 }
1405         }
1406 }
1407
1408 int LoadSDF_SDF(char *name)
1409 {
1410         /* This function reads uncompressed SPLAT Data Files (.sdf)
1411            containing digital elevation model data into memory.
1412            Elevation data, maximum and minimum elevations, and
1413            quadrangle limits are stored in the first available
1414            dem[] structure. */
1415
1416         int     x, y, data, indx, minlat, minlon, maxlat, maxlon;
1417         char    found, free_slot=0, line[20], sdf_file[255],
1418                 path_plus_name[255];
1419         FILE    *fd;
1420
1421         for (x=0; name[x]!='.' && name[x]!=0 && x<250; x++)
1422                 sdf_file[x]=name[x];
1423
1424         sdf_file[x]=0;
1425
1426         /* Parse filename for minimum latitude and longitude values */
1427
1428         sscanf(sdf_file,"%d:%d:%d:%d",&minlat,&maxlat,&minlon,&maxlon);
1429
1430         sdf_file[x]='.';
1431         sdf_file[x+1]='s';
1432         sdf_file[x+2]='d';
1433         sdf_file[x+3]='f';
1434         sdf_file[x+4]=0;
1435
1436         /* Is it already in memory? */
1437
1438         for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
1439         {
1440                 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west && maxlat==dem[indx].max_north && maxlon==dem[indx].max_west)
1441                         found=1;
1442         }
1443
1444         /* Is room available to load it? */
1445
1446         if (found==0)
1447         {
1448                 for (indx=0, free_slot=0; indx<MAXSLOTS && free_slot==0; indx++)
1449                         if (dem[indx].max_north==-90)
1450                                 free_slot=1;
1451         }
1452
1453         indx--;
1454
1455         if (free_slot && found==0 && indx>=0 && indx<MAXSLOTS)
1456         {
1457                 /* Search for SDF file in current working directory first */
1458
1459                 strncpy(path_plus_name,sdf_file,255);
1460
1461                 fd=fopen(path_plus_name,"rb");
1462
1463                 if (fd==NULL)
1464                 {
1465                         /* Next, try loading SDF file from path specified
1466                            in $HOME/.splat_path file or by -d argument */
1467
1468                         strncpy(path_plus_name,sdf_path,255);
1469                         strncat(path_plus_name,sdf_file,255);
1470
1471                         fd=fopen(path_plus_name,"rb");
1472                 }
1473
1474                 if (fd!=NULL)
1475                 {
1476                         fprintf(stdout,"Loading \"%s\" into slot %d...",path_plus_name,indx+1);
1477                         fflush(stdout);
1478
1479                         fgets(line,19,fd);
1480                         sscanf(line,"%d",&dem[indx].max_west);
1481
1482                         fgets(line,19,fd);
1483                         sscanf(line,"%d",&dem[indx].min_north);
1484
1485                         fgets(line,19,fd);
1486                         sscanf(line,"%d",&dem[indx].min_west);
1487
1488                         fgets(line,19,fd);
1489                         sscanf(line,"%d",&dem[indx].max_north);
1490
1491                         for (x=0; x<1200; x++)
1492                                 for (y=0; y<1200; y++)
1493                                 {
1494                                         fgets(line,19,fd);
1495                                         sscanf(line,"%d",&data);
1496
1497                                         dem[indx].data[x][y]=data;
1498
1499                                         if (data>dem[indx].max_el)
1500                                                 dem[indx].max_el=data;
1501
1502                                         if (data<dem[indx].min_el)
1503                                                 dem[indx].min_el=data;
1504                                 }
1505
1506                         fclose(fd);
1507
1508                         if (dem[indx].min_el<min_elevation)
1509                                 min_elevation=dem[indx].min_el;
1510
1511                         if (dem[indx].max_el>max_elevation)
1512                                 max_elevation=dem[indx].max_el;
1513
1514                         if (max_north==-90)
1515                                 max_north=dem[indx].max_north;
1516
1517                         else if (dem[indx].max_north>max_north)
1518                                 max_north=dem[indx].max_north;
1519
1520                         if (min_north==90)
1521                                 min_north=dem[indx].min_north;
1522
1523                         else if (dem[indx].min_north<min_north)
1524                                 min_north=dem[indx].min_north;
1525
1526                         if (max_west==-1)
1527                                 max_west=dem[indx].max_west;
1528
1529                         else
1530                         {
1531                                 if (abs(dem[indx].max_west-max_west)<180)
1532                                 {
1533                                         if (dem[indx].max_west>max_west)
1534                                                 max_west=dem[indx].max_west;
1535                                 }
1536
1537                                 else
1538                                 {
1539                                         if (dem[indx].max_west<max_west)
1540                                                 max_west=dem[indx].max_west;
1541                                 }
1542                         }
1543
1544                         if (min_west==360)
1545                                 min_west=dem[indx].min_west;
1546
1547                         else
1548                         {
1549                                 if (abs(dem[indx].min_west-min_west)<180)
1550                                 {
1551                                         if (dem[indx].min_west<min_west)
1552                                                 min_west=dem[indx].min_west;
1553                                 }
1554
1555                                 else
1556                                 {
1557                                         if (dem[indx].min_west>min_west)
1558                                                 min_west=dem[indx].min_west;
1559                                 }
1560                         }
1561
1562                         fprintf(stdout," Done!\n");
1563                         fflush(stdout);
1564                         return 1;
1565                 }
1566
1567                 else
1568                         return -1;
1569         }
1570
1571         else
1572                 return 0;
1573 }
1574
1575 char *BZfgets(BZFILE *bzfd, unsigned length)
1576 {
1577         /* This function returns at most one less than 'length' number
1578            of characters from a bz2 compressed file whose file descriptor
1579            is pointed to by *bzfd.  In operation, a buffer is filled with
1580            uncompressed data (size = BZBUFFER), which is then parsed
1581            and doled out as NULL terminated character strings every time
1582            this function is invoked.  A NULL string indicates an EOF
1583            or error condition. */
1584
1585         static int x, y, nBuf;
1586         static char buffer[BZBUFFER+1], output[BZBUFFER+1];
1587         char done=0;
1588
1589         if (opened!=1 && bzerror==BZ_OK)
1590         {
1591                 /* First time through.  Initialize everything! */
1592
1593                 x=0;
1594                 y=0;
1595                 nBuf=0;
1596                 opened=1;
1597                 output[0]=0;
1598         }
1599
1600         do
1601         {
1602                 if (x==nBuf && bzerror!=BZ_STREAM_END && bzerror==BZ_OK && opened)
1603                 {
1604                         /* Uncompress data into a static buffer */
1605
1606                         nBuf=BZ2_bzRead(&bzerror, bzfd, buffer, BZBUFFER);
1607                         buffer[nBuf]=0;
1608                         x=0;
1609                 }
1610
1611                 /* Build a string from buffer contents */
1612
1613                 output[y]=buffer[x];
1614
1615                 if (output[y]=='\n' || output[y]==0 || y==(int)length-1)
1616                 {
1617                         output[y+1]=0;
1618                         done=1;
1619                         y=0;
1620                 }
1621
1622                 else
1623                         y++;
1624                 x++;
1625
1626         } while (done==0);
1627
1628         if (output[0]==0)
1629                 opened=0;
1630
1631         return (output);
1632 }
1633
1634 int LoadSDF_BZ(char *name)
1635 {
1636         /* This function reads .bz2 compressed SPLAT Data Files containing
1637            digital elevation model data into memory.  Elevation data,
1638            maximum and minimum elevations, and quadrangle limits are
1639            stored in the first available dem[] structure. */
1640
1641         int     x, y, data, indx, minlat, minlon, maxlat, maxlon;
1642         char    found, free_slot=0, sdf_file[255], path_plus_name[255];
1643         FILE    *fd;
1644         BZFILE  *bzfd;
1645
1646         for (x=0; name[x]!='.' && name[x]!=0 && x<247; x++)
1647                 sdf_file[x]=name[x];
1648
1649         sdf_file[x]=0;
1650
1651         /* Parse sdf_file name for minimum latitude and longitude values */
1652
1653         sscanf(sdf_file,"%d:%d:%d:%d",&minlat,&maxlat,&minlon,&maxlon);
1654
1655         sdf_file[x]='.';
1656         sdf_file[x+1]='s';
1657         sdf_file[x+2]='d';
1658         sdf_file[x+3]='f';
1659         sdf_file[x+4]='.';
1660         sdf_file[x+5]='b';
1661         sdf_file[x+6]='z';
1662         sdf_file[x+7]='2';
1663         sdf_file[x+8]=0;
1664
1665         /* Is it already in memory? */
1666
1667         for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
1668         {
1669                 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west && maxlat==dem[indx].max_north && maxlon==dem[indx].max_west)
1670                         found=1;
1671         }
1672
1673         /* Is room available to load it? */
1674
1675         if (found==0)
1676         {
1677                 for (indx=0, free_slot=0; indx<MAXSLOTS && free_slot==0; indx++)
1678                         if (dem[indx].max_north==-90)
1679                                 free_slot=1;
1680         }
1681
1682         indx--;
1683
1684         if (free_slot && found==0 && indx>=0 && indx<MAXSLOTS)
1685         {
1686                 /* Search for SDF file in current working directory first */
1687
1688                 strncpy(path_plus_name,sdf_file,255);
1689
1690                 fd=fopen(path_plus_name,"rb");
1691                 bzfd=BZ2_bzReadOpen(&bzerror,fd,0,0,NULL,0);
1692
1693                 if (fd==NULL || bzerror!=BZ_OK)
1694                 {
1695                         /* Next, try loading SDF file from path specified
1696                            in $HOME/.splat_path file or by -d argument */
1697
1698                         strncpy(path_plus_name,sdf_path,255);
1699                         strncat(path_plus_name,sdf_file,255);
1700
1701                         fd=fopen(path_plus_name,"rb");
1702                         bzfd=BZ2_bzReadOpen(&bzerror,fd,0,0,NULL,0);
1703                 }
1704
1705                 if (fd!=NULL && bzerror==BZ_OK)
1706                 {
1707                         fprintf(stdout,"Loading \"%s\" into slot %d...",path_plus_name,indx+1);
1708                         fflush(stdout);
1709
1710                         sscanf(BZfgets(bzfd,255),"%d",&dem[indx].max_west);
1711                         sscanf(BZfgets(bzfd,255),"%d",&dem[indx].min_north);
1712                         sscanf(BZfgets(bzfd,255),"%d",&dem[indx].min_west);
1713                         sscanf(BZfgets(bzfd,255),"%d",&dem[indx].max_north);
1714
1715                         for (x=0; x<1200; x++)
1716                                 for (y=0; y<1200; y++)
1717                                 {
1718                                         sscanf(BZfgets(bzfd,20),"%d",&data);
1719
1720                                         dem[indx].data[x][y]=data;
1721
1722                                         if (data>dem[indx].max_el)
1723                                                 dem[indx].max_el=data;
1724
1725                                         if (data<dem[indx].min_el)
1726                                                 dem[indx].min_el=data;
1727                                 }
1728
1729                         fclose(fd);
1730
1731                         BZ2_bzReadClose(&bzerror,bzfd);
1732
1733                         if (dem[indx].min_el<min_elevation)
1734                                 min_elevation=dem[indx].min_el;
1735
1736                         if (dem[indx].max_el>max_elevation)
1737                                 max_elevation=dem[indx].max_el;
1738
1739                         if (max_north==-90)
1740                                 max_north=dem[indx].max_north;
1741
1742                         else if (dem[indx].max_north>max_north)
1743                                 max_north=dem[indx].max_north;
1744
1745                         if (min_north==90)
1746                                 min_north=dem[indx].min_north;
1747
1748                         else if (dem[indx].min_north<min_north)
1749                                 min_north=dem[indx].min_north;
1750
1751                         if (max_west==-1)
1752                                 max_west=dem[indx].max_west;
1753
1754                         else
1755                         {
1756                                 if (abs(dem[indx].max_west-max_west)<180)
1757                                 {
1758                                         if (dem[indx].max_west>max_west)
1759                                                 max_west=dem[indx].max_west;
1760                                 }
1761
1762                                 else
1763                                 {
1764                                         if (dem[indx].max_west<max_west)
1765                                                 max_west=dem[indx].max_west;
1766                                 }
1767                         }
1768
1769                         if (min_west==360)
1770                                 min_west=dem[indx].min_west;
1771
1772                         else
1773                         {
1774                                 if (abs(dem[indx].min_west-min_west)<180)
1775                                 {
1776                                         if (dem[indx].min_west<min_west)
1777                                                 min_west=dem[indx].min_west;
1778                                 }
1779
1780                                 else
1781                                 {
1782                                         if (dem[indx].min_west>min_west)
1783                                                 min_west=dem[indx].min_west;
1784                                 }
1785                         }
1786
1787                         fprintf(stdout," Done!\n");
1788                         fflush(stdout);
1789                         return 1;
1790                 }
1791
1792                 else
1793                         return -1;
1794         }
1795
1796         else
1797                 return 0;
1798 }
1799
1800 char LoadSDF(char *name)
1801 {
1802         /* This function loads the requested SDF file from the filesystem.
1803            It first tries to invoke the LoadSDF_SDF() function to load an
1804            uncompressed SDF file (since uncompressed files load slightly
1805            faster).  If that attempt fails, then it tries to load a
1806            compressed SDF file by invoking the LoadSDF_BZ() function.
1807            If that fails, then we can assume that no elevation data
1808            exists for the region requested, and that the region
1809            requested must be entirely over water. */
1810
1811         int     x, y, indx, minlat, minlon, maxlat, maxlon;
1812         char    found, free_slot=0;
1813         int     return_value=-1;
1814
1815         /* Try to load an uncompressed SDF first. */
1816
1817         return_value=LoadSDF_SDF(name);
1818
1819         /* If that fails, try loading a compressed SDF. */
1820
1821         if (return_value==0 || return_value==-1)
1822                 return_value=LoadSDF_BZ(name);
1823
1824         /* If neither format can be found, then assume the area is water. */
1825
1826         if (return_value==0 || return_value==-1)
1827         {
1828                 /* Parse SDF name for minimum latitude and longitude values */
1829
1830                 sscanf(name,"%d:%d:%d:%d",&minlat,&maxlat,&minlon,&maxlon);
1831
1832                 /* Is it already in memory? */
1833
1834                 for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
1835                 {
1836                         if (minlat==dem[indx].min_north && minlon==dem[indx].min_west && maxlat==dem[indx].max_north && maxlon==dem[indx].max_west)
1837                                 found=1;
1838                 }
1839
1840                 /* Is room available to load it? */
1841
1842                 if (found==0)
1843                 {
1844                         for (indx=0, free_slot=0; indx<MAXSLOTS && free_slot==0; indx++)
1845                                 if (dem[indx].max_north==-90)
1846                                         free_slot=1;
1847                 }
1848
1849                 indx--;
1850
1851                 if (free_slot && found==0 && indx>=0 && indx<MAXSLOTS)
1852                 {
1853                         fprintf(stdout,"Region  \"%s\" assumed as sea-level into slot %d...",name,indx+1);
1854                         fflush(stdout);
1855
1856                         dem[indx].max_west=maxlon;
1857                         dem[indx].min_north=minlat;
1858                         dem[indx].min_west=minlon;
1859                         dem[indx].max_north=maxlat;
1860
1861                         /* Fill DEM with sea-level topography */
1862
1863                         for (x=0; x<1200; x++)
1864                                 for (y=0; y<1200; y++)
1865                                 {
1866                                         dem[indx].data[x][y]=0;
1867
1868                                         if (dem[indx].min_el>0)
1869                                                 dem[indx].min_el=0;
1870                                 }
1871
1872                         if (dem[indx].min_el<min_elevation)
1873                                 min_elevation=dem[indx].min_el;
1874
1875                         if (dem[indx].max_el>max_elevation)
1876                                 max_elevation=dem[indx].max_el;
1877
1878                         if (max_north==-90)
1879                                 max_north=dem[indx].max_north;
1880
1881                         else if (dem[indx].max_north>max_north)
1882                                 max_north=dem[indx].max_north;
1883
1884                         if (min_north==90)
1885                                 min_north=dem[indx].min_north;
1886
1887                         else if (dem[indx].min_north<min_north)
1888                                 min_north=dem[indx].min_north;
1889
1890                         if (max_west==-1)
1891                                 max_west=dem[indx].max_west;
1892
1893                         else
1894                         {
1895                                 if (abs(dem[indx].max_west-max_west)<180)
1896                                 {
1897                                         if (dem[indx].max_west>max_west)
1898                                                 max_west=dem[indx].max_west;
1899                                 }
1900
1901                                 else
1902                                 {
1903                                         if (dem[indx].max_west<max_west)
1904                                                 max_west=dem[indx].max_west;
1905                                 }
1906                         }
1907
1908                         if (min_west==360)
1909                                 min_west=dem[indx].min_west;
1910
1911                         else
1912                         {
1913                                 if (abs(dem[indx].min_west-min_west)<180)
1914                                 {
1915                                         if (dem[indx].min_west<min_west)
1916                                                 min_west=dem[indx].min_west;
1917                                 }
1918
1919                                 else
1920                                 {
1921                                         if (dem[indx].min_west>min_west)
1922                                                 min_west=dem[indx].min_west;
1923                                 }
1924                         }
1925
1926                         fprintf(stdout," Done!\n");
1927                         fflush(stdout);
1928
1929                         return_value=1;
1930                 }
1931         }
1932
1933         return return_value;
1934 }
1935
1936 void LoadCities(char *filename)
1937 {
1938         /* This function reads SPLAT! city/site files, and plots
1939            the locations and names of the cities and site locations
1940            read on topographic maps generated by SPLAT! */
1941
1942         int     x, y, z;
1943         char    input[80], str[3][80];
1944         struct  site city_site;
1945         FILE    *fd=NULL;
1946
1947         fd=fopen(filename,"r");
1948
1949         if (fd!=NULL)
1950         {
1951                 fgets(input,78,fd);
1952
1953                 fprintf(stdout,"Reading \"%s\"... ",filename);
1954                 fflush(stdout);
1955
1956                 while (fd!=NULL && feof(fd)==0)
1957                 {
1958                         /* Parse line for name, latitude, and longitude */
1959
1960                         for (x=0, y=0, z=0; x<78 && input[x]!=0 && z<3; x++)
1961                         {
1962                                 if (input[x]!=',' && y<78)
1963                                 {
1964                                         str[z][y]=input[x];
1965                                         y++;
1966                                 }
1967
1968                                 else
1969                                 {
1970                                         str[z][y]=0;
1971                                         z++;
1972                                         y=0;
1973                                 }
1974                         }
1975
1976                         strncpy(city_site.name,str[0],49);
1977                         city_site.lat=ReadBearing(str[1]);
1978                         city_site.lon=ReadBearing(str[2]);
1979                         city_site.alt=0.0;
1980
1981                         PlaceMarker(city_site);
1982
1983                         fgets(input,78,fd);
1984                 }
1985
1986                 fclose(fd);
1987                 fprintf(stdout,"Done!\n");
1988                 fflush(stdout);
1989         }
1990
1991         else
1992                 fprintf(stderr,"*** ERROR: \"%s\": not found!\n",filename);
1993 }
1994
1995 void LoadUDT(char *filename)
1996 {
1997         /* This function reads a file containing User-Defined Terrain
1998            features for their addition to the digital elevation model
1999            data used by SPLAT!.  Elevations in the UDT file are evaluated
2000            and then copied into a temporary file under /tmp.  Then the
2001            contents of the temp file are scanned, and if found to be unique,
2002            are added to the ground elevations described by the digital
2003            elevation data already loaded into memory. */
2004
2005         int     i, x, y, z, fd=0;
2006         char    input[80], str[3][80], tempname[15], *pointer=NULL;
2007         double  latitude, longitude, height, templat, templon,
2008                 tempheight, one_pixel;
2009         FILE    *fd1=NULL, *fd2=NULL;
2010
2011         strcpy(tempname,"/tmp/XXXXXX\0");
2012         one_pixel=1.0/1200.0;
2013
2014         fd1=fopen(filename,"r");
2015
2016         if (fd1!=NULL)
2017         {
2018                 fd=mkstemp(tempname);
2019                 fd2=fopen(tempname,"w");
2020
2021                 fgets(input,78,fd1);
2022
2023                 pointer=strchr(input,';');
2024
2025                 if (pointer!=NULL)
2026                         *pointer=0;
2027
2028                 fprintf(stdout,"Reading \"%s\"... ",filename);
2029                 fflush(stdout);
2030
2031                 while (feof(fd1)==0)
2032                 {
2033                         /* Parse line for latitude, longitude, height */
2034
2035                         for (x=0, y=0, z=0; x<78 && input[x]!=0 && z<3; x++)
2036                         {
2037                                 if (input[x]!=',' && y<78)
2038                                 {
2039                                         str[z][y]=input[x];
2040                                         y++;
2041                                 }
2042
2043                                 else
2044                                 {
2045                                         str[z][y]=0;
2046                                         z++;
2047                                         y=0;
2048                                 }
2049                         }
2050
2051                         latitude=ReadBearing(str[0]);
2052                         longitude=ReadBearing(str[1]);
2053
2054                         /* Remove <CR> and/or <LF> from antenna height string */
2055
2056                         for (i=0; str[2][i]!=13 && str[2][i]!=10 && str[2][i]!=0; i++);
2057
2058                         str[2][i]=0;
2059
2060                         /* The terrain feature may be expressed in either
2061                            feet or meters.  If the letter 'M' or 'm' is
2062                            discovered in the string, then this is an
2063                            indication that the value given is expressed
2064                            in meters.  Otherwise the height is interpreted
2065                            as being expressed in feet.  */
2066
2067                         for (i=0; str[2][i]!='M' && str[2][i]!='m' && str[2][i]!=0 && i<48; i++);
2068
2069                         if (str[2][i]=='M' || str[2][i]=='m')
2070                         {
2071                                 str[2][i]=0;
2072                                 height=rint(atof(str[2]));
2073                         }
2074
2075                         else
2076                         {
2077                                 str[2][i]=0;
2078                                 height=rint(3.28084*atof(str[2]));
2079                         }
2080
2081                         if (height>0.0)
2082                                 fprintf(fd2,"%f, %f, %f\n",latitude, longitude, height);
2083
2084                         fgets(input,78,fd1);
2085
2086                         pointer=strchr(input,';');
2087
2088                         if (pointer!=NULL)
2089                                 *pointer=0;
2090                 }
2091
2092                 fclose(fd1);
2093                 fclose(fd2);
2094                 close(fd);
2095
2096                 fprintf(stdout,"Done!\n");
2097                 fflush(stdout);
2098
2099                 fd1=fopen(tempname,"r");
2100                 fd2=fopen(tempname,"r");
2101
2102                 fscanf(fd1,"%lf, %lf, %lf", &latitude, &longitude, &height);
2103
2104                 for (y=0; feof(fd1)==0; y++)
2105                 {
2106                         rewind(fd2);
2107
2108                         fscanf(fd2,"%lf, %lf, %lf", &templat, &templon, &tempheight);
2109
2110                         for (x=0, z=0; feof(fd2)==0; x++)
2111                         {
2112                                 if (x>y)
2113                                         if (fabs(latitude-templat)<=one_pixel && fabs(longitude-templon)<=one_pixel)
2114                                                 z=1;
2115
2116                                 fscanf(fd2,"%lf, %lf, %lf", &templat, &templon, &tempheight);
2117                         }
2118
2119                         if (z==0)
2120                                 AddElevation(latitude, longitude, height);
2121
2122                         fscanf(fd1,"%lf, %lf, %lf", &latitude, &longitude, &height);
2123                 }
2124
2125                 fclose(fd1);
2126                 fclose(fd2);
2127                 unlink(tempname);
2128         }
2129
2130         else
2131                 fprintf(stderr,"*** ERROR: \"%s\": not found!\n",filename);
2132 }
2133
2134 void LoadBoundaries(char *filename)
2135 {
2136         /* This function reads Cartographic Boundary Files available from
2137            the U.S. Census Bureau, and plots the data contained in those
2138            files on the PPM Map generated by SPLAT!.  Such files contain
2139            the coordinates that describe the boundaries of cities,
2140            counties, and states. */
2141
2142         int     x;
2143         double  lat0, lon0, lat1, lon1;
2144         char    string[80];
2145         struct  site source, destination;
2146         FILE    *fd=NULL;
2147
2148         fd=fopen(filename,"r");
2149
2150         if (fd!=NULL)
2151         {
2152                 fgets(string,78,fd);
2153
2154                 fprintf(stdout,"Reading \"%s\"... ",filename);
2155                 fflush(stdout);
2156
2157                 do
2158                 {
2159                         fgets(string,78,fd);
2160                         sscanf(string,"%lf %lf", &lon0, &lat0);
2161                         fgets(string,78,fd);
2162
2163                         do
2164                         {
2165                                 sscanf(string,"%lf %lf", &lon1, &lat1);
2166
2167                                 lon0=fabs(lon0);
2168                                 lon1=fabs(lon1);
2169
2170                                 source.lat=lat0;
2171                                 source.lon=lon0;
2172                                 destination.lat=lat1;
2173                                 destination.lon=lon1;
2174
2175                                 ReadPath(source,destination);
2176
2177                                 for (x=0; x<path.length; x++)
2178                                         OrMask(path.lat[x],path.lon[x],4);
2179
2180                                 lat0=lat1;
2181                                 lon0=lon1;
2182
2183                                 fgets(string,78,fd);
2184
2185                         } while (strncmp(string,"END",3)!=0 && feof(fd)==0);
2186
2187                         fgets(string,78,fd);
2188
2189                 } while (strncmp(string,"END",3)!=0 && feof(fd)==0);
2190
2191                 fclose(fd);
2192
2193                 fprintf(stdout,"Done!\n");
2194                 fflush(stdout);
2195         }
2196
2197         else
2198                 fprintf(stderr,"*** ERROR: \"%s\": not found!\n",filename);
2199 }
2200
2201 void ReadLRParm(char *txsite_filename)
2202 {
2203         /* This function reads Longley-Rice parameter data for the
2204            transmitter site.  The file name is the same as the txsite,
2205            except the filename extension is .lrp.  If the needed file
2206            is not found, then the file "splat.lrp" is read from the
2207            current working directory.  Failure to load this file will
2208            result in the default parameters hard coded into this
2209            function to be used and written to "splat.lrp". */
2210
2211         double  din;
2212         char    filename[255], string[80], *pointer=NULL;
2213         int     iin, ok=0, x;
2214         FILE    *fd=NULL, *outfile=NULL;
2215
2216         /* Default parameters in case things go bad */
2217
2218         LR.eps_dielect=15.0;
2219         LR.sgm_conductivity=0.005;
2220         LR.eno_ns_surfref=301.0;
2221         LR.frq_mhz=300.0;
2222         LR.radio_climate=5;
2223         LR.pol=0;
2224         LR.conf=0.50;
2225         LR.rel=0.50;
2226
2227         /* Modify txsite filename to one with a .lrp extension. */
2228
2229         strncpy(filename,txsite_filename,255);
2230
2231         for (x=0; filename[x]!='.' && filename[x]!=0 && filename[x]!='\n' && x<249; x++);
2232
2233         filename[x]='.';
2234         filename[x+1]='l';
2235         filename[x+2]='r';
2236         filename[x+3]='p';
2237         filename[x+4]=0;
2238
2239         fd=fopen(filename,"r");
2240
2241         if (fd==NULL)
2242         {
2243                 /* Load default "splat.lrp" file */
2244
2245                 strncpy(filename,"splat.lrp\0",10);
2246                 fd=fopen(filename,"r");
2247         }
2248
2249         if (fd!=NULL)
2250         {
2251                 fgets(string,80,fd);
2252
2253                 pointer=strchr(string,';');
2254
2255                 if (pointer!=NULL)
2256                         *pointer=0;
2257
2258                 ok=sscanf(string,"%lf", &din);
2259
2260                 if (ok)
2261                 {
2262                         LR.eps_dielect=din;
2263
2264                         fgets(string,80,fd);
2265
2266                         pointer=strchr(string,';');
2267
2268                         if (pointer!=NULL)
2269                                 *pointer=0;
2270
2271                         ok=sscanf(string,"%lf", &din);
2272                 }
2273
2274                 if (ok)
2275                 {
2276                         LR.sgm_conductivity=din;
2277
2278                         fgets(string,80,fd);
2279
2280                         pointer=strchr(string,';');
2281
2282                         if (pointer!=NULL)
2283                                 *pointer=0;
2284
2285                         ok=sscanf(string,"%lf", &din);
2286                 }
2287
2288                 if (ok)
2289                 {
2290                         LR.eno_ns_surfref=din;
2291
2292                         fgets(string,80,fd);
2293
2294                         pointer=strchr(string,';');
2295
2296                         if (pointer!=NULL)
2297                                 *pointer=0;
2298
2299                         ok=sscanf(string,"%lf", &din);
2300                 }
2301
2302                 if (ok)
2303                 {
2304                         LR.frq_mhz=din;
2305
2306                         fgets(string,80,fd);
2307
2308                         pointer=strchr(string,';');
2309
2310                         if (pointer!=NULL)
2311                                 *pointer=0;
2312
2313                         ok=sscanf(string,"%d", &iin);
2314                 }
2315
2316                 if (ok)
2317                 {
2318                         LR.radio_climate=iin;
2319
2320                         fgets(string,80,fd);
2321
2322                         pointer=strchr(string,';');
2323
2324                         if (pointer!=NULL)
2325                                 *pointer=0;
2326
2327                         ok=sscanf(string,"%d", &iin);
2328                 }
2329
2330                 if (ok)
2331                 {
2332                         LR.pol=iin;
2333
2334                         fgets(string,80,fd);
2335
2336                         pointer=strchr(string,';');
2337
2338                         if (pointer!=NULL)
2339                                 *pointer=0;
2340
2341                         ok=sscanf(string,"%lf", &din);
2342                 }
2343
2344                 if (ok)
2345                 {
2346                         LR.conf=din;
2347
2348                         fgets(string,80,fd);
2349
2350                         pointer=strchr(string,';');
2351
2352                         if (pointer!=NULL)
2353                                 *pointer=0;
2354
2355                         ok=sscanf(string,"%lf", &din);
2356                 }
2357
2358                 fclose(fd);
2359
2360                 if (ok)
2361                 {
2362                         LR.rel=din;
2363                         LoadPAT(filename);
2364                 }
2365         }
2366
2367         if (fd==NULL)
2368         {
2369                 /* Create a "splat.lrp" file since one
2370                    could not be successfully loaded. */
2371
2372                 outfile=fopen("splat.lrp","w");
2373
2374                 fprintf(outfile,"%.3f\t; Earth Dielectric Constant (Relative permittivity)\n",LR.eps_dielect);
2375                 fprintf(outfile,"%.3f\t; Earth Conductivity (Siemens per meter)\n", LR.sgm_conductivity);
2376                 fprintf(outfile,"%.3f\t; Atmospheric Bending Constant (N-Units)\n",LR.eno_ns_surfref);
2377                 fprintf(outfile,"%.3f\t; Frequency in MHz (20 MHz to 20 GHz)\n", LR.frq_mhz);
2378                 fprintf(outfile,"%d\t; Radio Climate\n",LR.radio_climate);
2379                 fprintf(outfile,"%d\t; Polarization (0 = Horizontal, 1 = Vertical)\n", LR.pol);
2380                 fprintf(outfile,"%.2f\t; Fraction of situations\n",LR.conf);
2381                 fprintf(outfile, "%.2f\t; Fraction of time\n",LR.rel);
2382                 fprintf(outfile,"\nPlease consult SPLAT! documentation for the meaning and use of this data.\n");
2383
2384                 fclose(outfile);
2385
2386                 fprintf(stderr,"\n%c*** There were problems reading your \"%s\" file! ***\nA \"splat.lrp\" file was written to your directory with default data.\n",7,filename);
2387         }
2388
2389         if (fd==NULL || ok==0)
2390                 fprintf(stderr,"Longley-Rice default parameters have been assumed for this analysis.\n");
2391 }
2392
2393 void PlotPath(struct site source, struct site destination, char mask_value)
2394 {
2395         /* This function analyzes the path between the source and
2396            destination locations.  It determines which points along
2397            the path have line-of-sight visibility to the source.
2398            Points along with path having line-of-sight visibility
2399            to the source at an AGL altitude equal to that of the
2400            destination location are stored by setting bit 1 in the
2401            mask[][] array, which are displayed in green when PPM
2402            maps are later generated by SPLAT!. */
2403
2404         char block;
2405         int x, y;
2406         register double cos_xmtr_angle, cos_test_angle, test_alt;
2407         double distance, rx_alt, tx_alt;
2408
2409         ReadPath(source,destination);
2410
2411         for (y=0; y<path.length; y++)
2412         {
2413                 /* Test this point only if it hasn't been already
2414                    tested and found to be free of obstructions. */
2415
2416                 if ((GetMask(path.lat[y],path.lon[y])&mask_value)==0)
2417                 {
2418                         distance=5280.0*path.distance[y];
2419                         tx_alt=earthradius+source.alt+path.elevation[0];
2420                         rx_alt=earthradius+destination.alt+path.elevation[y];
2421
2422                         /* Calculate the cosine of the elevation of the
2423                            transmitter as seen at the temp rx point. */
2424
2425                         cos_xmtr_angle=((rx_alt*rx_alt)+(distance*distance)-(tx_alt*tx_alt))/(2.0*rx_alt*distance);
2426
2427                         for (x=y, block=0; x>=0 && block==0; x--)
2428                         {
2429                                 distance=5280.0*(path.distance[y]-path.distance[x]);
2430                                 test_alt=earthradius+path.elevation[x];
2431
2432                                 cos_test_angle=((rx_alt*rx_alt)+(distance*distance)-(test_alt*test_alt))/(2.0*rx_alt*distance);
2433
2434                                 /* Compare these two angles to determine if
2435                                    an obstruction exists.  Since we're comparing
2436                                    the cosines of these angles rather than
2437                                    the angles themselves, the following "if"
2438                                    statement is reversed from what it would
2439                                    be if the actual angles were compared. */
2440
2441                                 if (cos_xmtr_angle>cos_test_angle)
2442                                         block=1;
2443                         }
2444
2445                         if (block==0)
2446                                 OrMask(path.lat[y],path.lon[y],mask_value);
2447                 }
2448         }
2449 }
2450
2451 void PlotLRPath(struct site source, struct site destination, FILE *fd)
2452 {
2453         /* This function plots the RF path loss between source and
2454            destination points based on the Longley-Rice propagation
2455            model, taking into account antenna pattern data, if available. */
2456
2457         char    block=0, strmode[100];
2458         int     x, y, errnum;
2459         double  loss, azimuth, pattern=0.0,
2460                 source_alt, dest_alt, source_alt2, dest_alt2,
2461                 cos_xmtr_angle, cos_test_angle=0.0, test_alt,
2462                 elevation, distance=0.0, four_thirds_earth;
2463         struct  site temp;
2464
2465         ReadPath(source,destination);
2466
2467         four_thirds_earth=EARTHRADIUS*(4.0/3.0);
2468
2469         /* Copy elevations along path into the elev_l[] array. */
2470
2471         for (x=0; x<path.length; x++)
2472                 elev_l[x+2]=path.elevation[x]*METERS_PER_FOOT;
2473
2474         /* Since the only energy the Longley-Rice model considers
2475            reaching the destination is based on what is scattered
2476            or deflected from the first obstruction along the path,
2477            we first need to find the location and elevation angle
2478            of that first obstruction (if it exists).  This is done
2479            using a 4/3rds Earth radius to match the model used by
2480            Longley-Rice.  This information is required for properly
2481            integrating the antenna's elevation pattern into the
2482            calculation for overall path loss.  (Using path.length-1
2483            below avoids a Longley-Rice model error from occuring at
2484            the destination point.) */
2485
2486         for (y=2; (y<(path.length-1) && path.distance[y]<=max_range); y++)
2487         {
2488                 /* Process this point only if it
2489                    has not already been processed. */
2490
2491                 if (GetMask(path.lat[y],path.lon[y])==0)
2492                 {
2493                         distance=5280.0*path.distance[y];
2494                         source_alt=four_thirds_earth+source.alt+path.elevation[0];
2495                         dest_alt=four_thirds_earth+destination.alt+path.elevation[y];
2496                         dest_alt2=dest_alt*dest_alt;
2497                         source_alt2=source_alt*source_alt;
2498
2499                         /* Calculate the cosine of the elevation of
2500                            the receiver as seen by the transmitter. */
2501
2502                         cos_xmtr_angle=((source_alt2)+(distance*distance)-(dest_alt2))/(2.0*source_alt*distance);
2503
2504                         if (got_elevation_pattern || fd!=NULL)
2505                         {
2506                                 /* If no antenna elevation pattern is available, and
2507                                    no output file is designated, the following code
2508                                    that determines the elevation angle to the first
2509                                    obstruction along the path is bypassed. */
2510
2511                                 for (x=2, block=0; (x<y && block==0); x++)
2512                                 {
2513                                         distance=5280.0*path.distance[x];
2514                                         test_alt=four_thirds_earth+path.elevation[x];
2515
2516                                         /* Calculate the cosine of the elevation
2517                                            angle of the terrain (test point)
2518                                            as seen by the transmitter. */
2519
2520                                         cos_test_angle=((source_alt2)+(distance*distance)-(test_alt*test_alt))/(2.0*source_alt*distance);
2521
2522                                         /* Compare these two angles to determine if
2523                                            an obstruction exists.  Since we're comparing
2524                                            the cosines of these angles rather than
2525                                            the angles themselves, the sense of the
2526                                            following "if" statement is reversed from
2527                                            what it would be if the angles themselves
2528                                            were compared. */
2529
2530                                         if (cos_xmtr_angle>cos_test_angle)
2531                                                 block=1;
2532                                 }
2533
2534                                 /* At this point, we have the elevation angle
2535                                    to the first obstruction (if it exists). */
2536                         }
2537
2538                         /* Determine attenuation for each point along the
2539                            path using Longley-Rice's point_to_point mode
2540                            starting at y=2 (number_of_points = 1), the
2541                            shortest distance terrain can play a role in
2542                            path loss. */
2543
2544                         elev_l[0]=y-1;  /* (number of points - 1) */
2545
2546                         /* Distance between elevation samples */
2547                         elev_l[1]=METERS_PER_MILE*(path.distance[y]-path.distance[y-1]);
2548
2549                         point_to_point(elev_l,source.alt*METERS_PER_FOOT,
2550                         destination.alt*METERS_PER_FOOT, LR.eps_dielect,
2551                         LR.sgm_conductivity, LR.eno_ns_surfref, LR.frq_mhz,
2552                         LR.radio_climate, LR.pol, LR.conf, LR.rel, loss,
2553                         strmode, errnum);
2554
2555                         if (block)
2556                                 elevation=((acos(cos_test_angle))/deg2rad)-90.0;
2557
2558                         else
2559                                 elevation=((acos(cos_xmtr_angle))/deg2rad)-90.0;
2560
2561                         temp.lat=path.lat[y];
2562                         temp.lon=path.lon[y];
2563
2564                         azimuth=(Azimuth(source,temp));
2565
2566                         if (fd!=NULL)
2567                         {
2568                                 /* Write path loss data to output file */
2569
2570                                 fprintf(fd,"%.7f, %.7f, %.3f, %.3f, %.2f\n",path.lat[y], path.lon[y], azimuth, elevation, loss);
2571                         }
2572
2573                         /* Integrate the antenna's radiation
2574                            pattern into the overall path loss. */
2575
2576                         x=(int)rint(10.0*(10.0-elevation));
2577
2578                         if (x>=0 && x<=1000)
2579                         {
2580                                 azimuth=rint(azimuth);
2581
2582                                 pattern=(double)LR.antenna_pattern[(int)azimuth][x];
2583
2584                                 if (pattern!=0.0)
2585                                 {
2586                                         pattern=20.0*log10(pattern);
2587                                         loss-=pattern;
2588                                 }
2589                         }
2590
2591                         if (loss>225.0)
2592                                 loss=225.0;
2593
2594                         if (loss<75.0)
2595                                 loss=75.0;
2596
2597                         loss-=75.0;
2598                         loss/=10.0;
2599                         loss+=1.0;
2600
2601                         OrMask(path.lat[y],path.lon[y],((unsigned char)(loss))<<3);
2602                 }
2603
2604                 else if (GetMask(path.lat[y],path.lon[y])==0 && path.distance[y]>max_range)
2605                         OrMask(path.lat[y],path.lon[y],1);
2606         }
2607 }
2608
2609 void PlotCoverage(struct site source, double altitude)
2610 {
2611         /* This function performs a 360 degree sweep around the
2612            transmitter site (source location), and plots the
2613            line-of-sight coverage of the transmitter on the SPLAT!
2614            generated topographic map based on a receiver located
2615            at the specified altitude (in feet AGL).  Results
2616            are stored in memory, and written out in the form
2617            of a topographic map when the WritePPM() function
2618            is later invoked. */
2619
2620         float lat, lon, one_pixel;
2621         static unsigned char mask_value;
2622         int z, count;
2623         struct site edge;
2624         unsigned char symbol[4], x;
2625
2626         /* Initialize mask_value */
2627
2628         if (mask_value!=8 && mask_value!=16 && mask_value!=32)
2629                 mask_value=1;
2630
2631         one_pixel=1.0/1200.0;
2632
2633         symbol[0]='.';
2634         symbol[1]='o';
2635         symbol[2]='O';
2636         symbol[3]='o';
2637
2638         count=0;
2639
2640         fprintf(stdout,"\nComputing line-of-sight coverage of %s with an RX antenna\nat %.2f %s AGL:\n\n 0%c to  25%c ",source.name,metric?altitude*METERS_PER_FOOT:altitude,metric?"meters":"feet",37,37);
2641         fflush(stdout);
2642
2643         /* 18.75=1200 pixels/degree divided by 64 loops
2644            per progress indicator symbol (.oOo) printed. */
2645
2646         z=(int)(18.75*ReduceAngle(max_west-min_west));
2647
2648         for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2649         {
2650                 if (lon>=360.0)
2651                         lon-=360.0;
2652
2653                 edge.lat=max_north;
2654                 edge.lon=lon;
2655                 edge.alt=altitude;
2656
2657                 PlotPath(source,edge,mask_value);
2658                 count++;
2659
2660                 if (count==z)
2661                 {
2662                         fprintf(stdout,"%c",symbol[x]);
2663                         fflush(stdout);
2664                         count=0;
2665
2666                         if (x==3)
2667                                 x=0;
2668                         else
2669                                 x++;
2670                 }
2671         }
2672
2673         count=0;
2674         fprintf(stdout,"\n25%c to  50%c ",37,37);
2675         fflush(stdout);
2676
2677         z=(int)(18.75*(max_north-min_north));
2678
2679         for (lat=max_north, x=0; lat>=min_north; lat-=one_pixel)
2680         {
2681                 edge.lat=lat;
2682                 edge.lon=min_west;
2683                 edge.alt=altitude;
2684
2685                 PlotPath(source,edge,mask_value);
2686                 count++;
2687
2688                 if (count==z)
2689                 {
2690                         fprintf(stdout,"%c",symbol[x]);
2691                         fflush(stdout);
2692                         count=0;
2693
2694                         if (x==3)
2695                                 x=0;
2696                         else
2697                                 x++;
2698                 }
2699         }
2700
2701         count=0;
2702         fprintf(stdout,"\n50%c to  75%c ",37,37);
2703         fflush(stdout);
2704
2705         z=(int)(18.75*ReduceAngle(max_west-min_west));
2706
2707         for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2708         {
2709                 if (lon>=360.0)
2710                         lon-=360.0;
2711
2712                 edge.lat=min_north;
2713                 edge.lon=lon;
2714                 edge.alt=altitude;
2715
2716                 PlotPath(source,edge,mask_value);
2717                 count++;
2718
2719                 if (count==z)
2720                 {
2721                         fprintf(stdout,"%c",symbol[x]);
2722                         fflush(stdout);
2723                         count=0;
2724
2725                         if (x==3)
2726                                 x=0;
2727                         else
2728                                 x++;
2729                 }
2730         }
2731
2732         count=0;
2733         fprintf(stdout,"\n75%c to 100%c ",37,37);
2734         fflush(stdout);
2735
2736         z=(int)(18.75*(max_north-min_north));
2737
2738         for (lat=min_north, x=0; lat<=max_north; lat+=one_pixel)
2739         {
2740                 edge.lat=lat;
2741                 edge.lon=max_west;
2742                 edge.alt=altitude;
2743
2744                 PlotPath(source,edge,mask_value);
2745                 count++;
2746
2747                 if (count==z)
2748                 {
2749                         fprintf(stdout,"%c",symbol[x]);
2750                         fflush(stdout);
2751                         count=0;
2752
2753                         if (x==3)
2754                                 x=0;
2755                         else
2756                                 x++;
2757                 }
2758         }
2759
2760         fprintf(stdout,"\nDone!\n");
2761         fflush(stdout);
2762
2763         /* Assign next mask value */
2764
2765         switch (mask_value)
2766         {
2767                 case 1:
2768                         mask_value=8;
2769                         break;
2770
2771                 case 8:
2772                         mask_value=16;
2773                         break;
2774
2775                 case 16:
2776                         mask_value=32;
2777         }
2778 }
2779
2780 void PlotLRMap(struct site source, double altitude, char *plo_filename)
2781 {
2782         /* This function performs a 360 degree sweep around the
2783            transmitter site (source location), and plots the
2784            Longley-Rice attenuation on the SPLAT! generated
2785            topographic map based on a receiver located at
2786            the specified altitude (in feet AGL).  Results
2787            are stored in memory, and written out in the form
2788            of a topographic map when the WritePPMLR() function
2789            is later invoked. */
2790
2791         int z, count;
2792         struct site edge;
2793         float lat, lon, one_pixel;
2794         unsigned char symbol[4], x;
2795         FILE *fd=NULL;
2796
2797         one_pixel=1.0/1200.0;
2798
2799         symbol[0]='.';
2800         symbol[1]='o';
2801         symbol[2]='O';
2802         symbol[3]='o';
2803
2804         count=0;
2805
2806         fprintf(stdout,"\nComputing Longley-Rice coverage of %s ", source.name);
2807
2808         fprintf(stdout,"out to a radius\nof %.2f %s with an RX antenna at %.2f %s AGL:\n\n 0%c to  25%c ",metric?max_range*KM_PER_MILE:max_range,metric?"kilometers":"miles",metric?altitude*METERS_PER_FOOT:altitude,metric?"meters":"feet",37,37);
2809         fflush(stdout);
2810
2811         if (plo_filename[0]!=0)
2812                 fd=fopen(plo_filename,"wb");
2813
2814         if (fd!=NULL)
2815         {
2816                 /* Write header information to output file */
2817
2818                 fprintf(fd,"%d, %d\t; max_west, min_west\n%d, %d\t; max_north, min_north\n",max_west, min_west, max_north, min_north);
2819         }
2820
2821         /* 18.75=1200 pixels/degree divided by 64 loops
2822            per progress indicator symbol (.oOo) printed. */
2823
2824         z=(int)(18.75*ReduceAngle(max_west-min_west));
2825
2826         for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2827         {
2828                 if (lon>=360.0)
2829                         lon-=360.0;
2830
2831                 edge.lat=max_north;
2832                 edge.lon=lon;
2833                 edge.alt=altitude;
2834
2835                 PlotLRPath(source,edge,fd);
2836                 count++;
2837
2838                 if (count==z)
2839                 {
2840                         fprintf(stdout,"%c",symbol[x]);
2841                         fflush(stdout);
2842                         count=0;
2843
2844                         if (x==3)
2845                                 x=0;
2846                         else
2847                                 x++;
2848                 }
2849         }
2850
2851         count=0;
2852         fprintf(stdout,"\n25%c to  50%c ",37,37);
2853         fflush(stdout);
2854
2855         z=(int)(18.75*(max_north-min_north));
2856
2857         for (lat=max_north, x=0; lat>=min_north; lat-=one_pixel)
2858         {
2859                 edge.lat=lat;
2860                 edge.lon=min_west;
2861                 edge.alt=altitude;
2862
2863                 PlotLRPath(source,edge,fd);
2864                 count++;
2865
2866                 if (count==z)
2867                 {
2868                         fprintf(stdout,"%c",symbol[x]);
2869                         fflush(stdout);
2870                         count=0;
2871
2872                         if (x==3)
2873                                 x=0;
2874                         else
2875                                 x++;
2876                 }
2877         }
2878
2879         count=0;
2880         fprintf(stdout,"\n50%c to  75%c ",37,37);
2881         fflush(stdout);
2882
2883         z=(int)(18.75*ReduceAngle(max_west-min_west));
2884
2885         for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2886         {
2887                 if (lon>=360.0)
2888                         lon-=360.0;
2889
2890                 edge.lat=min_north;
2891                 edge.lon=lon;
2892                 edge.alt=altitude;
2893
2894                 PlotLRPath(source,edge,fd);
2895                 count++;
2896
2897                 if (count==z)
2898                 {
2899                         fprintf(stdout,"%c",symbol[x]);
2900                         fflush(stdout);
2901                         count=0;
2902
2903                         if (x==3)
2904                                 x=0;
2905                         else
2906                                 x++;
2907                 }
2908         }
2909
2910         count=0;
2911         fprintf(stdout,"\n75%c to 100%c ",37,37);
2912         fflush(stdout);
2913
2914         z=(int)(18.75*(max_north-min_north));
2915
2916         for (lat=min_north, x=0; lat<=max_north; lat+=one_pixel)
2917         {
2918                 edge.lat=lat;
2919                 edge.lon=max_west;
2920                 edge.alt=altitude;
2921
2922                 PlotLRPath(source,edge,fd);
2923                 count++;
2924
2925                 if (count==z)
2926                 {
2927                         fprintf(stdout,"%c",symbol[x]);
2928                         fflush(stdout);
2929                         count=0;
2930
2931                         if (x==3)
2932                                 x=0;
2933                         else
2934                                 x++;
2935                 }
2936         }
2937
2938         if (fd!=NULL)
2939                 fclose(fd);
2940
2941         fprintf(stdout,"\nDone!\n");
2942         fflush(stdout);
2943 }
2944
2945 void WritePPM(char *filename, unsigned char geo)
2946 {
2947         /* This function generates a topographic map in Portable Pix Map
2948            (PPM) format based on logarithmically scaled topology data,
2949            as well as the content of flags held in the mask[][] array.
2950            The image created is rotated counter-clockwise 90 degrees
2951            from its representation in dem[][] so that north points
2952            up and east points right in the image generated. */
2953
2954         char mapfile[255], geofile[255];
2955         unsigned char found, mask;
2956         unsigned width, height, output;
2957         int indx, x, y, x0=0, y0=0, minlat, minlon;
2958         float lat, lon, one_pixel, conversion, one_over_gamma;
2959         FILE *fd;
2960
2961         one_pixel=1.0/1200.0;
2962         one_over_gamma=1.0/GAMMA;
2963         conversion=255.0/pow((double)(max_elevation-min_elevation),one_over_gamma);
2964
2965         width=(unsigned)(1200*ReduceAngle(max_west-min_west));
2966         height=(unsigned)(1200*ReduceAngle(max_north-min_north));
2967
2968         if (filename[0]==0)
2969                 strncpy(mapfile, "map.ppm\0",8);
2970         else
2971         {
2972                 for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
2973                 {
2974                         mapfile[x]=filename[x];
2975                         geofile[x]=filename[x];
2976                 }
2977
2978                 mapfile[x]='.';
2979                 geofile[x]='.';
2980                 mapfile[x+1]='p';
2981                 geofile[x+1]='g';
2982                 mapfile[x+2]='p';
2983                 geofile[x+2]='e';
2984                 mapfile[x+3]='m';
2985                 geofile[x+3]='o';
2986                 mapfile[x+4]=0;
2987                 geofile[x+4]=0;
2988         }
2989
2990         if (geo)
2991         {
2992                 fd=fopen(geofile,"wb");
2993
2994                 fprintf(fd,"FILENAME\t%s\n",mapfile);
2995                 fprintf(fd,"#\t\tX\tY\tLong\t\tLat\n");
2996                 fprintf(fd,"TIEPOINT\t0\t0\t%d.000\t\t%d.000\n",(max_west<180?-max_west:360-max_west),max_north);
2997                 fprintf(fd,"TIEPOINT\t%u\t%u\t%d.000\t\t%d.000\n",width-1,height-1,(min_west<180?-min_west:360-min_west),min_north);
2998                 fprintf(fd,"IMAGESIZE\t%u\t%u\n",width,height);
2999                 fprintf(fd,"#\n# Auto Generated by SPLAT! v%s\n#\n",splat_version);
3000
3001                 fclose(fd);
3002         }
3003
3004         fd=fopen(mapfile,"wb");
3005
3006         fprintf(fd,"P6\n%u %u\n255\n",width,height);
3007
3008         fprintf(stdout,"\nWriting \"%s\" (%ux%u pixmap image)... ",mapfile,width,height);
3009         fflush(stdout);
3010
3011         for (y=0, lat=((double)max_north)-one_pixel; y<(int)height; y++, lat-=one_pixel)
3012         {
3013                 minlat=(int)floor(lat);
3014
3015                 for (x=0, lon=((double)max_west)-one_pixel; x<(int)width; x++, lon-=one_pixel)
3016                 {
3017                         if (lon<0.0)
3018                                 lon+=360.0;
3019
3020                         minlon=(int)floor(lon);
3021
3022                         for (indx=0, found=0; indx<MAXSLOTS && found==0;)
3023                                 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
3024                                         found=1;
3025                                 else
3026                                         indx++;
3027
3028                         if (found)
3029                         {
3030                                 x0=(int)(1199.0*(lat-floor(lat)));
3031                                 y0=(int)(1199.0*(lon-floor(lon)));
3032
3033                                 mask=dem[indx].mask[x0][y0];
3034
3035                                 if (mask&2)
3036                                         /* Text Labels: Red */
3037                                         fprintf(fd,"%c%c%c",255,0,0);
3038
3039                                 else if (mask&4)
3040                                         /* County Boundaries: Light Cyan */
3041                                         fprintf(fd,"%c%c%c",128,128,255);
3042
3043                                 else switch (mask&57)
3044                                 {
3045                                         case 1:
3046                                         /* TX1: Green */
3047                                         fprintf(fd,"%c%c%c",0,255,0);
3048                                         break;
3049
3050                                         case 8:
3051                                         /* TX2: Cyan */
3052                                         fprintf(fd,"%c%c%c",0,255,255);
3053                                         break;
3054
3055                                         case 9:
3056                                         /* TX1 + TX2: Yellow */
3057                                         fprintf(fd,"%c%c%c",255,255,0);
3058                                         break;
3059
3060                                         case 16:
3061                                         /* TX3: Medium Violet */
3062                                         fprintf(fd,"%c%c%c",147,112,219);
3063                                         break;
3064
3065                                         case 17:
3066                                         /* TX1 + TX3: Pink */
3067                                         fprintf(fd,"%c%c%c",255,192,203);
3068                                         break;
3069
3070                                         case 24:
3071                                         /* TX2 + TX3: Orange */
3072                                         fprintf(fd,"%c%c%c",255,165,0);
3073                                         break;
3074
3075                                         case 25:
3076                                         /* TX1 + TX2 + TX3: Dark Green */
3077                                         fprintf(fd,"%c%c%c",0,100,0);
3078                                         break;
3079
3080                                         case 32:
3081                                         /* TX4: Sienna 1 */
3082                                         fprintf(fd,"%c%c%c",255,130,71);
3083                                         break;
3084
3085                                         case 33:
3086                                         /* TX1 + TX4: Green Yellow */
3087                                         fprintf(fd,"%c%c%c",173,255,47);
3088                                         break;
3089
3090                                         case 40:
3091                                         /* TX2 + TX4: Dark Sea Green 1 */
3092                                         fprintf(fd,"%c%c%c",193,255,193);
3093                                         break;
3094
3095                                         case 41:
3096                                         /* TX1 + TX2 + TX4: Blanched Almond */
3097                                         fprintf(fd,"%c%c%c",255,235,205);
3098                                         break;
3099
3100                                         case 48:
3101                                         /* TX3 + TX4: Dark Turquoise */
3102                                         fprintf(fd,"%c%c%c",0,206,209);
3103                                         break;
3104
3105                                         case 49:
3106                                         /* TX1 + TX3 + TX4: Medium Spring Green */
3107                                         fprintf(fd,"%c%c%c",0,250,154);
3108                                         break;
3109
3110                                         case 56:
3111                                         /* TX2 + TX3 + TX4: Tan */
3112                                         fprintf(fd,"%c%c%c",210,180,140);
3113                                         break;
3114
3115                                         case 57:
3116                                         /* TX1 + TX2 + TX3 + TX4: Gold2 */
3117                                         fprintf(fd,"%c%c%c",238,201,0);
3118                                         break;
3119
3120                                         default:
3121                                         /* Water: Medium Blue */
3122                                         if (dem[indx].data[x0][y0]==0)
3123                                                 fprintf(fd,"%c%c%c",0,0,170);
3124                                         else
3125                                         {
3126                                                 /* Elevation: Greyscale */
3127                                                 output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
3128                                                 fprintf(fd,"%c%c%c",output,output,output);
3129                                         }
3130                                 }
3131                         }
3132
3133                         else
3134                         {
3135                                 /* We should never get here, but if */
3136                                 /* we do, display the region as black */
3137
3138                                 fprintf(fd,"%c%c%c",0,0,0);
3139                         }
3140                 }
3141         }
3142
3143         fclose(fd);
3144         fprintf(stdout,"Done!\n");
3145         fflush(stdout);
3146 }
3147
3148 void WritePPMLR(char *filename, unsigned char geo)
3149 {
3150         /* This function generates a topographic map in Portable Pix Map
3151            (PPM) format based on the content of flags held in the mask[][]
3152            array (only).  The image created is rotated counter-clockwise
3153            90 degrees from its representation in dem[][] so that north
3154            points up and east points right in the image generated. */
3155
3156         char mapfile[255], geofile[255];
3157         unsigned width, height, output;
3158         unsigned char found, mask, cityorcounty;
3159         int indx, x, y, t, t2, x0, y0, minlat, minlon, loss;
3160         float lat, lon, one_pixel, conversion, one_over_gamma;
3161         FILE *fd;
3162
3163         one_pixel=1.0/1200.0;
3164         one_over_gamma=1.0/GAMMA;
3165         conversion=255.0/pow((double)(max_elevation-min_elevation),one_over_gamma);
3166
3167         width=(unsigned)(1200*ReduceAngle(max_west-min_west));
3168         height=(unsigned)(1200*ReduceAngle(max_north-min_north));
3169
3170         if (filename[0]==0)
3171                 strncpy(mapfile, "map.ppm\0",8);
3172         else
3173         {
3174                 for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
3175                 {
3176                         mapfile[x]=filename[x];
3177                         geofile[x]=filename[x];
3178                 }
3179
3180                 mapfile[x]='.';
3181                 geofile[x]='.';
3182                 mapfile[x+1]='p';
3183                 geofile[x+1]='g';
3184                 mapfile[x+2]='p';
3185                 geofile[x+2]='e';
3186                 mapfile[x+3]='m';
3187                 geofile[x+3]='o';
3188                 mapfile[x+4]=0;
3189                 geofile[x+4]=0;
3190         }
3191
3192         if (geo)
3193         {
3194                 fd=fopen(geofile,"wb");
3195
3196                 fprintf(fd,"FILENAME\t%s\n",mapfile);
3197                 fprintf(fd,"#\t\tX\tY\tLong\t\tLat\n");
3198                 fprintf(fd,"TIEPOINT\t0\t0\t%d.000\t\t%d.000\n",(max_west<180?-max_west:360-max_west),max_north);
3199                 fprintf(fd,"TIEPOINT\t%u\t%u\t%d.000\t\t%.3f\n",width-1,height+29,(min_west<180?-min_west:360-min_west),(double)(min_north-0.025));
3200                 fprintf(fd,"IMAGESIZE\t%u\t%u\n",width,height+30);
3201                 fprintf(fd,"#\n# Auto Generated by SPLAT! v%s\n#\n",splat_version);
3202
3203                 fclose(fd);
3204         }
3205
3206         fd=fopen(mapfile,"wb");
3207
3208         fprintf(fd,"P6\n%u %u\n255\n",width,height+30);
3209
3210         fprintf(stdout,"\nWriting \"%s\" (%ux%u pixmap image)... ",mapfile,width,height+30);
3211         fflush(stdout);
3212
3213         for (y=0, lat=((double)max_north)-one_pixel; y<(int)height; y++, lat-=one_pixel)
3214         {
3215                 minlat=(int)floor(lat);
3216
3217                 for (x=0, lon=((double)max_west)-one_pixel; x<(int)width; x++, lon-=one_pixel)
3218                 {
3219                         if (lon<0.0)
3220                                 lon+=360.0;
3221
3222                         minlon=(int)floor(lon);
3223
3224                         for (indx=0, found=0; indx<MAXSLOTS && found==0;)
3225                                 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
3226                                         found=1;
3227                                 else
3228                                         indx++;
3229                         if (found)
3230                         {
3231                                 x0=(int)(1199.0*(lat-floor(lat)));
3232                                 y0=(int)(1199.0*(lon-floor(lon)));
3233
3234                                 mask=dem[indx].mask[x0][y0];
3235                                 loss=70+(10*(int)((mask&248)>>3));
3236                                 cityorcounty=0;
3237
3238                                 if (mask&2)
3239                                 {
3240                                         /* Text Labels - Black or Red */
3241
3242                                         if ((mask&120) && (loss<=90))
3243                                                 fprintf(fd,"%c%c%c",0,0,0);
3244                                         else
3245                                                 fprintf(fd,"%c%c%c",255,0,0);
3246
3247                                         cityorcounty=1;
3248                                 }
3249
3250                                 else if (mask&4)
3251                                 {
3252                                         /* County Boundaries: Black */
3253
3254                                         fprintf(fd,"%c%c%c",0,0,0);
3255
3256                                         cityorcounty=1;
3257                                 }
3258
3259                                 if (cityorcounty==0)
3260                                 {
3261                                         if (loss>maxdB)
3262
3263                                         { /* Display land or sea elevation */
3264
3265                                                 if (dem[indx].data[x0][y0]==0)
3266                                                         fprintf(fd,"%c%c%c",0,0,170);
3267                                                 else
3268                                                 {
3269                                                         /* Elevation: Greyscale */
3270                                                         output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
3271                                                         fprintf(fd,"%c%c%c",output,output,output);
3272                                                 }
3273                                         }
3274
3275                                         else switch (loss)
3276                                         {
3277                                                 /* Plot signal loss in color */
3278
3279                                                 case 80:
3280                                                 fprintf(fd,"%c%c%c",255,0,0);
3281                                                 break;
3282
3283                                                 case 90:
3284                                                 fprintf(fd,"%c%c%c",255,128,0);
3285                                                 break;
3286
3287                                                 case 100:
3288                                                 fprintf(fd,"%c%c%c",255,165,0);
3289                                                 break;
3290
3291                                                 case 110:
3292                                                 fprintf(fd,"%c%c%c",255,206,0);
3293                                                 break;
3294
3295                                                 case 120:
3296                                                 fprintf(fd,"%c%c%c",255,255,0);
3297                                                 break;
3298
3299                                                 case 130:
3300                                                 fprintf(fd,"%c%c%c",184,255,0);
3301                                                 break;
3302
3303                                                 case 140:
3304                                                 fprintf(fd,"%c%c%c",0,255,0);
3305                                                 break;
3306
3307                                                 case 150:
3308                                                 fprintf(fd,"%c%c%c",0,208,0);
3309                                                 break;
3310
3311                                                 case 160:
3312                                                 fprintf(fd,"%c%c%c",0,196,196);
3313                                                 break;
3314
3315                                                 case 170:
3316                                                 fprintf(fd,"%c%c%c",0,148,255);
3317                                                 break;
3318
3319                                                 case 180:
3320                                                 fprintf(fd,"%c%c%c",80,80,255);
3321                                                 break;
3322
3323                                                 case 190:
3324                                                 fprintf(fd,"%c%c%c",0,38,255);
3325                                                 break;
3326
3327                                                 case 200:
3328                                                 fprintf(fd,"%c%c%c",142,63,255);
3329                                                 break;
3330
3331                                                 case 210:
3332                                                 fprintf(fd,"%c%c%c",196,54,255);
3333                                                 break;
3334
3335                                                 case 220:
3336                                                 fprintf(fd,"%c%c%c",255,0,255);
3337                                                 break;
3338
3339                                                 case 230:
3340                                                 fprintf(fd,"%c%c%c",255,194,204);
3341                                                 break;
3342
3343                                                 default:
3344
3345                                                 if (dem[indx].data[x0][y0]==0)
3346                                                         fprintf(fd,"%c%c%c",0,0,170);
3347                                                 else
3348                                                 {
3349                                                         /* Elevation: Greyscale */
3350                                                         output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
3351                                                         fprintf(fd,"%c%c%c",output,output,output);
3352                                                 }
3353                                         }
3354                                 }
3355                         }
3356
3357                         else
3358                         {
3359                                 /* We should never get here, but if */
3360                                 /* we do, display the region as black */
3361
3362                                 fprintf(fd,"%c%c%c",0,0,0);
3363                         }
3364                 }
3365         }
3366
3367         /* Display legend along bottom of image */
3368
3369         x0=width/16;
3370
3371         for (y0=0; y0<30; y0++)
3372         {
3373                 for (indx=0; indx<16; indx++)
3374                 {
3375                         for (x=0; x<x0; x++)
3376                         {
3377                                 t=indx;
3378                                 t2=indx+8;
3379
3380                                 if (y0>=10 && y0<=18)
3381                                 {
3382                                         if (t2>9)
3383                                         {
3384                                                 if (x>=11 && x<=17)
3385                                                         if (smallfont[t2/10][y0-10][x-11])
3386                                                                 t=255;
3387                                         }
3388
3389                                         if (x>=19 && x<=25)
3390                                                 if (smallfont[t2%10][y0-10][x-19])
3391                                                         t=255;
3392
3393                                         if (x>=27 && x<=33)
3394                                                 if (smallfont[0][y0-10][x-27])
3395                                                         t=255;
3396                                 }
3397
3398                                 switch (t)
3399                                 {
3400                                         case 0:
3401                                         fprintf(fd,"%c%c%c",255,0,0);
3402                                         break;
3403
3404                                         case 1:
3405                                         fprintf(fd,"%c%c%c",255,128,0);
3406                                         break;
3407
3408                                         case 2:
3409                                         fprintf(fd,"%c%c%c",255,165,0);
3410                                         break;
3411
3412                                         case 3:
3413                                         fprintf(fd,"%c%c%c",255,206,0);
3414                                         break;
3415
3416                                         case 4:
3417                                         fprintf(fd,"%c%c%c",255,255,0);
3418                                         break;
3419
3420                                         case 5:
3421                                         fprintf(fd,"%c%c%c",184,255,0);
3422                                         break;
3423
3424                                         case 6:
3425                                         fprintf(fd,"%c%c%c",0,255,0);
3426                                         break;
3427
3428                                         case 7:
3429                                         fprintf(fd,"%c%c%c",0,208,0);
3430                                         break;
3431
3432                                         case 8:
3433                                         fprintf(fd,"%c%c%c",0,196,196);
3434                                         break;
3435
3436                                         case 9:
3437                                         fprintf(fd,"%c%c%c",0,148,255);
3438                                         break;
3439
3440                                         case 10:
3441                                         fprintf(fd,"%c%c%c",80,80,255);
3442                                         break;
3443
3444                                         case 11:
3445                                         fprintf(fd,"%c%c%c",0,38,255);
3446                                         break;
3447
3448                                         case 12:
3449                                         fprintf(fd,"%c%c%c",142,63,255);
3450                                         break;
3451
3452                                         case 13:
3453                                         fprintf(fd,"%c%c%c",196,54,255);
3454                                         break;
3455
3456                                         case 14:
3457                                         fprintf(fd,"%c%c%c",255,0,255);
3458                                         break;
3459
3460                                         case 255:
3461                                         /* Black */
3462                                         fprintf(fd,"%c%c%c",0,0,0);
3463                                         break;
3464
3465                                         default:
3466                                         fprintf(fd,"%c%c%c",255,194,204);
3467                                 }
3468                         }
3469                 }
3470         }
3471
3472         fclose(fd);
3473         fprintf(stdout,"Done!\n");
3474         fflush(stdout);
3475 }
3476
3477 void GraphTerrain(struct site source, struct site destination, char *name)
3478 {
3479         /* This function invokes gnuplot to generate an appropriate
3480            output file indicating the terrain profile between the source
3481            and destination locations.  "filename" is the name assigned
3482            to the output file generated by gnuplot.  The filename extension
3483            is used to set gnuplot's terminal setting and output file type.
3484            If no extension is found, .png is assumed.  */
3485
3486         int     x, y, z;
3487         char    filename[255], term[30], ext[15];
3488         FILE    *fd=NULL;
3489
3490         ReadPath(destination,source);
3491
3492         fd=fopen("profile.gp","wb");
3493
3494         for (x=0; x<path.length; x++)
3495         {
3496                 if (metric)
3497                         fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*path.elevation[x]);
3498
3499                 else
3500                         fprintf(fd,"%f\t%f\n",path.distance[x],path.elevation[x]);
3501         }
3502
3503         fclose(fd);
3504
3505         if (name[0]==0)
3506         {
3507                 /* Default filename and output file type */
3508
3509                 strncpy(filename,"profile\0",8);
3510                 strncpy(term,"png\0",4);
3511                 strncpy(ext,"png\0",4);
3512         }
3513
3514         else
3515         {
3516                 /* Grab extension and terminal type from "name" */
3517
3518                 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
3519                         filename[x]=name[x];
3520
3521                 if (name[x]=='.')
3522                 {
3523                         for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
3524                         {
3525                                 term[y]=tolower(name[x]);
3526                                 ext[y]=term[y];
3527                         }
3528
3529                         ext[y]=0;
3530                         term[y]=0;
3531                         filename[z]=0;
3532                 }
3533
3534                 else
3535                 {       /* No extension -- Default is png */
3536
3537                         filename[x]=0;
3538                         strncpy(term,"png\0",4);
3539                         strncpy(ext,"png\0",4);
3540                 }
3541         }
3542
3543         /* Either .ps or .postscript may be used
3544            as an extension for postscript output. */
3545
3546         if (strncmp(term,"postscript",10)==0)
3547                 strncpy(ext,"ps\0",3);
3548
3549         else if (strncmp(ext,"ps",2)==0)
3550                 strncpy(term,"postscript enhanced color\0",26);
3551
3552         fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
3553         fflush(stdout);
3554
3555         fd=fopen("splat.gp","w");
3556         fprintf(fd,"set grid\n");
3557         fprintf(fd,"set autoscale\n");
3558         fprintf(fd,"set encoding iso_8859_1\n");
3559         fprintf(fd,"set term %s\n",term);
3560         fprintf(fd,"set title \"SPLAT! Terrain Profile Between %s and %s (%.2f%c Azimuth)\"\n",destination.name, source.name, Azimuth(destination,source),176);
3561
3562         if (metric)
3563         {
3564                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(source,destination));
3565                 fprintf(fd,"set ylabel \"Ground Elevation Above Sea Level (meters)\"\n");
3566
3567
3568         }
3569
3570         else
3571         {
3572                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(source,destination));
3573                 fprintf(fd,"set ylabel \"Ground Elevation Above Sea Level (feet)\"\n");
3574         }
3575
3576         fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
3577         fprintf(fd,"plot \"profile.gp\" title \"\" with lines\n");
3578         fclose(fd);
3579
3580         x=system("gnuplot splat.gp");
3581
3582         if (x!=-1)
3583         {
3584                 unlink("splat.gp");
3585                 unlink("profile.gp");
3586                 fprintf(stdout," Done!\n");
3587                 fflush(stdout);
3588         }
3589
3590         else
3591                 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
3592 }
3593
3594 void GraphElevation(struct site source, struct site destination, char *name)
3595 {
3596         /* This function invokes gnuplot to generate an appropriate
3597            output file indicating the terrain profile between the source
3598            and destination locations.  "filename" is the name assigned
3599            to the output file generated by gnuplot.  The filename extension
3600            is used to set gnuplot's terminal setting and output file type.
3601            If no extension is found, .png is assumed. */
3602
3603         int     x, y, z;
3604         char    filename[255], term[30], ext[15];
3605         double  angle, refangle, maxangle=-90.0;
3606         struct  site remote;
3607         FILE    *fd=NULL, *fd2=NULL;
3608
3609         ReadPath(destination,source);  /* destination=RX, source=TX */
3610         refangle=ElevationAngle(destination,source);
3611
3612         fd=fopen("profile.gp","wb");
3613         fd2=fopen("reference.gp","wb");
3614
3615         for (x=1; x<path.length-1; x++)
3616         {
3617                 remote.lat=path.lat[x];
3618                 remote.lon=path.lon[x];
3619                 remote.alt=0.0;
3620                 angle=ElevationAngle(destination,remote);
3621
3622                 if (metric)
3623                 {
3624                         fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[x],angle);
3625                         fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[x],refangle);
3626                 }
3627
3628                 else
3629                 {
3630                         fprintf(fd,"%f\t%f\n",path.distance[x],angle);
3631                         fprintf(fd2,"%f\t%f\n",path.distance[x],refangle);
3632                 }
3633
3634                 if (angle>maxangle)
3635                         maxangle=angle;
3636         }
3637
3638         if (metric)
3639         {
3640                 fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],refangle);
3641                 fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],refangle);
3642         }
3643
3644         else
3645         {
3646                 fprintf(fd,"%f\t%f\n",path.distance[path.length-1],refangle);
3647                 fprintf(fd2,"%f\t%f\n",path.distance[path.length-1],refangle);
3648         }
3649
3650         fclose(fd);
3651         fclose(fd2);
3652
3653         if (name[0]==0)
3654         {
3655                 /* Default filename and output file type */
3656
3657                 strncpy(filename,"profile\0",8);
3658                 strncpy(term,"png\0",4);
3659                 strncpy(ext,"png\0",4);
3660         }
3661
3662         else
3663         {
3664                 /* Grab extension and terminal type from "name" */
3665
3666                 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
3667                         filename[x]=name[x];
3668
3669                 if (name[x]=='.')
3670                 {
3671                         for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
3672                         {
3673                                 term[y]=tolower(name[x]);
3674                                 ext[y]=term[y];
3675                         }
3676
3677                         ext[y]=0;
3678                         term[y]=0;
3679                         filename[z]=0;
3680                 }
3681
3682                 else
3683                 {       /* No extension -- Default is png */
3684
3685                         filename[x]=0;
3686                         strncpy(term,"png\0",4);
3687                         strncpy(ext,"png\0",4);
3688                 }
3689         }
3690
3691         /* Either .ps or .postscript may be used
3692            as an extension for postscript output. */
3693
3694         if (strncmp(term,"postscript",10)==0)
3695                 strncpy(ext,"ps\0",3);
3696
3697         else if (strncmp(ext,"ps",2)==0)
3698                 strncpy(term,"postscript enhanced color\0",26);
3699
3700         fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
3701         fflush(stdout);
3702
3703         fd=fopen("splat.gp","w");
3704
3705         fprintf(fd,"set grid\n");
3706         fprintf(fd,"set yrange [%2.3f to %2.3f]\n", (-fabs(refangle)-0.25), maxangle+0.25);
3707         fprintf(fd,"set encoding iso_8859_1\n");
3708         fprintf(fd,"set term %s\n",term);
3709         fprintf(fd,"set title \"SPLAT! Elevation Profile Between %s and %s (%.2f%c azimuth)\"\n",destination.name,source.name,Azimuth(destination,source),176);
3710
3711         if (metric)
3712                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(source,destination));
3713         else
3714                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(source,destination));
3715
3716
3717         fprintf(fd,"set ylabel \"Elevation Angle Along LOS Path Between %s and %s (degrees)\"\n",destination.name,source.name);
3718         fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
3719         fprintf(fd,"plot \"profile.gp\" title \"Real Earth Profile\" with lines, \"reference.gp\" title \"Line of Sight Path (%.2f%c elevation)\" with lines\n",refangle,176);
3720
3721         fclose(fd);
3722
3723         x=system("gnuplot splat.gp");
3724
3725         if (x!=-1)
3726         {
3727                 unlink("splat.gp");
3728                 unlink("profile.gp");
3729                 unlink("reference.gp");
3730
3731                 fprintf(stdout," Done!\n");
3732                 fflush(stdout);
3733         }
3734
3735         else
3736                 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
3737 }
3738
3739 void GraphHeight(struct site source, struct site destination, char *name, double f, unsigned char n)
3740 {
3741         /* This function invokes gnuplot to generate an appropriate
3742            output file indicating the terrain profile between the source
3743            and destination locations referenced to the line-of-sight path
3744            between the receive and transmit sites.  "filename" is the name
3745            assigned to the output file generated by gnuplot.  The filename
3746            extension is used to set gnuplot's terminal setting and output
3747            file type.  If no extension is found, .png is assumed. */
3748
3749         int     x, y, z;
3750         char    filename[255], term[30], ext[15];
3751         double  a, b, c, height=0.0, refangle, cangle, maxheight=-100000.0,
3752                 minheight=100000.0, lambda=0.0, f_zone=0.0, fpt6_zone=0.0,
3753                 nm=0.0, nb=0.0, ed=0.0, es=0.0, r=0.0, d=0.0, d1=0.0,
3754                 terrain, azimuth, distance, dheight=0.0, minterrain=100000.0,
3755                 minearth=100000.0, miny, maxy, min2y, max2y;
3756         struct  site remote;
3757         FILE    *fd=NULL, *fd2=NULL, *fd3=NULL, *fd4=NULL, *fd5=NULL;
3758
3759         ReadPath(destination,source);  /* destination=RX, source=TX */
3760         azimuth=Azimuth(destination,source);
3761         distance=Distance(destination,source);
3762         refangle=ElevationAngle(destination,source);
3763         b=GetElevation(destination)+destination.alt+earthradius;
3764
3765         /* Wavelength and path distance (great circle) in feet. */
3766
3767         if (f)
3768         {
3769                 lambda=9.8425e8/(f*1e6);
3770                 d=5280.0*path.distance[path.length-1];
3771         }
3772
3773         if (n)
3774         {
3775                 ed=GetElevation(destination);
3776                 es=GetElevation(source);
3777                 nb=-destination.alt-ed;
3778                 nm=(-source.alt-es-nb)/(path.distance[path.length-1]);
3779         }
3780
3781         fd=fopen("profile.gp","wb");
3782         fd2=fopen("reference.gp","wb");
3783         fd5=fopen("curvature.gp", "wb");
3784
3785         if (f)
3786         {
3787                 fd3=fopen("fresnel.gp", "wb");
3788                 fd4=fopen("fresnel_pt_6.gp", "wb");
3789         }
3790
3791         for (x=0; x<path.length; x++)
3792         {
3793                 remote.lat=path.lat[x];
3794                 remote.lon=path.lon[x];
3795                 remote.alt=0.0;
3796
3797                 terrain=GetElevation(remote);
3798
3799                 if (x==0)
3800                         terrain+=destination.alt;  /* RX antenna spike */
3801
3802                 a=terrain+earthradius;
3803                 cangle=5280.0*Distance(destination,remote)/earthradius;
3804                 c=b*sin(refangle*deg2rad+HALFPI)/sin(HALFPI-refangle*deg2rad-cangle);
3805
3806                 height=a-c;
3807
3808                 /* Per Fink and Christiansen, Electronics
3809                  * Engineers' Handbook, 1989:
3810                  *
3811                  *   H = sqrt(lamba * d1 * (d - d1)/d)
3812                  *
3813                  * where H is the distance from the LOS
3814                  * path to the first Fresnel zone boundary.
3815                  */
3816
3817                 if (f)
3818                 {
3819                         d1=5280.0*path.distance[x];
3820                         f_zone=-1*sqrt(lambda*d1*(d-d1)/d);
3821                         fpt6_zone=0.6*f_zone;
3822                 }
3823
3824                 if (n)
3825                 {
3826                         r=-(nm*path.distance[x])-nb;
3827                         height+=r;
3828
3829                         if (f>0)
3830                         {
3831                                 f_zone+=r;
3832                                 fpt6_zone+=r;
3833                         }
3834                 }
3835
3836                 else
3837                         r=0.0;
3838
3839                 if (metric)
3840                 {
3841                         fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*height);
3842                         fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*r);
3843                         fprintf(fd5,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*(height-terrain));
3844                 }
3845
3846                 else
3847                 {
3848                         fprintf(fd,"%f\t%f\n",path.distance[x],height);
3849                         fprintf(fd2,"%f\t%f\n",path.distance[x],r);
3850                         fprintf(fd5,"%f\t%f\n",path.distance[x],height-terrain);
3851                 }
3852
3853                 if (f)
3854                 {
3855                         if (metric)
3856                         {
3857                                 fprintf(fd3,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*f_zone);
3858                                 fprintf(fd4,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*fpt6_zone);
3859                         }
3860
3861                         else
3862                         {
3863                                 fprintf(fd3,"%f\t%f\n",path.distance[x],f_zone);
3864                                 fprintf(fd4,"%f\t%f\n",path.distance[x],fpt6_zone);
3865                         }
3866
3867                         if (f_zone<minheight)
3868                                 minheight=f_zone;
3869                 }
3870
3871                 if (height>maxheight)
3872                         maxheight=height;
3873
3874                 if (height<minheight)
3875                         minheight=height;
3876
3877                 if (r>maxheight)
3878                         maxheight=r;
3879
3880                 if (terrain<minterrain)
3881                         minterrain=terrain;
3882
3883                 if ((height-terrain)<minearth)
3884                         minearth=height-terrain;
3885         }
3886
3887         if (n)
3888                 r=-(nm*path.distance[path.length-1])-nb;
3889         else
3890                 r=0.0;
3891
3892         if (metric)
3893         {
3894                 fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3895                 fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3896         }
3897
3898         else
3899         {
3900                 fprintf(fd,"%f\t%f\n",path.distance[path.length-1],r);
3901                 fprintf(fd2,"%f\t%f\n",path.distance[path.length-1],r);
3902         }
3903
3904         if (f)
3905         {
3906                 if (metric)
3907                 {
3908                         fprintf(fd3,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3909                         fprintf(fd4,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3910                 }
3911
3912                 else
3913                 {
3914                         fprintf(fd3,"%f\t%f\n",path.distance[path.length-1],r);
3915                         fprintf(fd4,"%f\t%f\n",path.distance[path.length-1],r);
3916                 }
3917         }
3918
3919         if (r>maxheight)
3920                 maxheight=r;
3921
3922         if (r<minheight)
3923                 minheight=r;
3924
3925         fclose(fd);
3926         fclose(fd2);
3927         fclose(fd5);
3928
3929         if (f)
3930         {
3931                 fclose(fd3);
3932                 fclose(fd4);
3933         }
3934
3935         if (name[0]==0)
3936         {
3937                 /* Default filename and output file type */
3938
3939                 strncpy(filename,"height\0",8);
3940                 strncpy(term,"png\0",4);
3941                 strncpy(ext,"png\0",4);
3942         }
3943
3944         else
3945         {
3946                 /* Grab extension and terminal type from "name" */
3947
3948                 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
3949                         filename[x]=name[x];
3950
3951                 if (name[x]=='.')
3952                 {
3953                         for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
3954                         {
3955                                 term[y]=tolower(name[x]);
3956                                 ext[y]=term[y];
3957                         }
3958
3959                         ext[y]=0;
3960                         term[y]=0;
3961                         filename[z]=0;
3962                 }
3963
3964                 else
3965                 {       /* No extension -- Default is png */
3966
3967                         filename[x]=0;
3968                         strncpy(term,"png\0",4);
3969                         strncpy(ext,"png\0",4);
3970                 }
3971         }
3972
3973         /* Either .ps or .postscript may be used
3974            as an extension for postscript output. */
3975
3976         if (strncmp(term,"postscript",10)==0)
3977                 strncpy(ext,"ps\0",3);
3978
3979         else if (strncmp(ext,"ps",2)==0)
3980                 strncpy(term,"postscript enhanced color\0",26);
3981
3982         fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
3983         fflush(stdout);
3984
3985         fd=fopen("splat.gp","w");
3986
3987         dheight=maxheight-minheight;
3988         miny=minheight-0.15*dheight;
3989         maxy=maxheight+0.05*dheight;
3990
3991         if (maxy<20.0)
3992                 maxy=20.0;
3993
3994         dheight=maxheight-minheight;
3995         min2y=miny-minterrain+0.05*dheight;
3996
3997         if (minearth<min2y)
3998         {
3999                 miny-=min2y-minearth+0.05*dheight;
4000                 min2y=minearth-0.05*dheight;
4001         }
4002
4003         max2y=min2y+maxy-miny;
4004
4005         fprintf(fd,"set grid\n");
4006         fprintf(fd,"set yrange [%2.3f to %2.3f]\n", metric?miny*METERS_PER_FOOT:miny, metric?maxy*METERS_PER_FOOT:maxy);
4007         fprintf(fd,"set y2range [%2.3f to %2.3f]\n", metric?min2y*METERS_PER_FOOT:min2y, metric?max2y*METERS_PER_FOOT:max2y);
4008         fprintf(fd,"set xrange [-0.5 to %2.3f]\n",metric?KM_PER_MILE*rint(distance+0.5):rint(distance+0.5));
4009         fprintf(fd,"set encoding iso_8859_1\n");
4010         fprintf(fd,"set term %s\n",term);
4011
4012         if (f)
4013                 fprintf(fd,"set title \"SPLAT! Path Profile Between %s and %s (%.2f%c azimuth)\\nWith First Fresnel Zone\"\n",destination.name, source.name, azimuth,176);
4014
4015         else
4016                 fprintf(fd,"set title \"SPLAT! Height Profile Between %s and %s (%.2f%c azimuth)\"\n",destination.name, source.name, azimuth,176);
4017
4018         if (metric)
4019                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(source,destination));
4020         else
4021                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(source,destination));
4022
4023         if (n)
4024         {
4025                 if (metric)
4026                         fprintf(fd,"set ylabel \"Normalized Height Referenced To LOS Path Between\\n%s and %s (meters)\"\n",destination.name,source.name);
4027
4028                 else
4029                         fprintf(fd,"set ylabel \"Normalized Height Referenced To LOS Path Between\\n%s and %s (feet)\"\n",destination.name,source.name);
4030
4031         }
4032
4033         else
4034         {
4035                 if (metric)
4036                         fprintf(fd,"set ylabel \"Height Referenced To LOS Path Between %s and %s (meters)\"\n",destination.name,source.name);
4037
4038                 else
4039                         fprintf(fd,"set ylabel \"Height Referenced To LOS Path Between %s and %s (feet)\"\n",destination.name,source.name);
4040         }
4041
4042         fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
4043
4044         if (f)
4045                 fprintf(fd,"plot \"profile.gp\" title \"Point-to-Point Profile\" with lines, \"reference.gp\" title \"Line of Sight Path\" with lines, \"curvature.gp\" axes x1y2 title \"Earth's Curvature Contour\" with lines, \"fresnel.gp\" axes x1y1 title \"First Fresnel Zone (%.3f MHz)\" with lines, \"fresnel_pt_6.gp\" title \"60%% of First Fresnel Zone\" with lines\n",f);
4046
4047         else
4048                 fprintf(fd,"plot \"profile.gp\" title \"Point-to-Point Profile\" with lines, \"reference.gp\" title \"Line Of Sight Path\" with lines, \"curvature.gp\" axes x1y2 title \"Earth's Curvature Contour\" with lines\n");
4049
4050         fclose(fd);
4051
4052         x=system("gnuplot splat.gp");
4053
4054         if (x!=-1)
4055         {
4056                 unlink("splat.gp");
4057                 unlink("profile.gp");
4058                 unlink("reference.gp");
4059                 unlink("curvature.gp");
4060
4061                 if (f)
4062                 {
4063                         unlink("fresnel.gp");
4064                         unlink("fresnel_pt_6.gp");
4065                 }
4066
4067                 fprintf(stdout," Done!\n");
4068                 fflush(stdout);
4069         }
4070
4071         else
4072                 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
4073 }
4074
4075 void GraphLongley(struct site source, struct site destination, char *name)
4076 {
4077         /* This function invokes gnuplot to generate an appropriate
4078            output file indicating the Longley-Rice model loss between
4079            the source and destination locations.  "filename" is the
4080            name assigned to the output file generated by gnuplot.
4081            The filename extension is used to set gnuplot's terminal
4082            setting and output file type.  If no extension is found,
4083            .png is assumed. */
4084
4085         int     x, y, z, errnum, errflag=0;
4086         char    filename[255], term[30], ext[15], strmode[100],
4087                 report_name[80], block=0;
4088         double  maxloss=-100000.0, minloss=100000.0, loss, haavt,
4089                 angle1, angle2, azimuth, pattern=1.0, patterndB=0.0,
4090                 total_loss=0.0, cos_xmtr_angle, cos_test_angle=0.0,
4091                 source_alt, test_alt, dest_alt, source_alt2, dest_alt2,
4092                 distance, elevation, four_thirds_earth;
4093         FILE    *fd=NULL, *fd2=NULL;
4094
4095         sprintf(report_name,"%s-to-%s.lro",source.name,destination.name);
4096
4097         four_thirds_earth=EARTHRADIUS*(4.0/3.0);
4098
4099         for (x=0; report_name[x]!=0; x++)
4100                 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
4101                         report_name[x]='_';
4102
4103         fd2=fopen(report_name,"w");
4104
4105         fprintf(fd2,"\n\t--==[ SPLAT! v%s Longley-Rice Model Path Loss Report ]==--\n\n",splat_version);
4106         fprintf(fd2,"Analysis of RF path conditions between %s and %s:\n",source.name, destination.name);
4107         fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4108         fprintf(fd2,"Transmitter site: %s\n",source.name);
4109
4110         if (source.lat>=0.0)
4111         {
4112                 fprintf(fd2,"Site location: %.4f North / %.4f West",source.lat, source.lon);
4113                 fprintf(fd2, " (%s N / ", dec2dms(source.lat));
4114         }
4115
4116         else
4117         {
4118
4119                 fprintf(fd2,"Site location: %.4f South / %.4f West",-source.lat, source.lon);
4120                 fprintf(fd2, " (%s S / ", dec2dms(source.lat));
4121         }
4122
4123         fprintf(fd2, "%s W)\n", dec2dms(source.lon));
4124
4125         if (metric)
4126         {
4127                 fprintf(fd2,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(source));
4128                 fprintf(fd2,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*source.alt,METERS_PER_FOOT*(source.alt+GetElevation(source)));
4129         }
4130
4131         else
4132         {
4133                 fprintf(fd2,"Ground elevation: %.2f feet AMSL\n",GetElevation(source));
4134                 fprintf(fd2,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",source.alt, source.alt+GetElevation(source));
4135         }
4136
4137         haavt=haat(source);
4138
4139         if (haavt>-4999.0)
4140         {
4141                 if (metric)
4142                         fprintf(fd2,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4143                 else
4144                         fprintf(fd2,"Antenna height above average terrain: %.2f feet\n",haavt);
4145         }
4146
4147         azimuth=Azimuth(source,destination);
4148         angle1=ElevationAngle(source,destination);
4149         angle2=ElevationAngle2(source,destination,earthradius);
4150
4151         if (got_azimuth_pattern || got_elevation_pattern)
4152         {
4153                 x=(int)rint(10.0*(10.0-angle2));
4154
4155                 if (x>=0 && x<=1000)
4156                         pattern=(double)LR.antenna_pattern[(int)rint(azimuth)][x];
4157
4158                 patterndB=20.0*log10(pattern);
4159
4160                 fprintf(fd2,"Antenna pattern between %s and %s: %.3f (%.2f dB)\n",source.name, destination.name, pattern, patterndB);
4161
4162         }
4163
4164         if (metric)
4165                 fprintf(fd2,"Distance to %s: %.2f kilometers\n",destination.name,METERS_PER_FOOT*Distance(source,destination));
4166
4167         else
4168                 fprintf(fd2,"Distance to %s: %.2f miles.\n",destination.name,Distance(source,destination));
4169
4170         fprintf(fd2,"Azimuth to %s: %.2f degrees\n",destination.name,azimuth);
4171
4172         if (angle1>=0.0)
4173                 fprintf(fd2,"Elevation angle to %s: %+.4f degrees\n",destination.name,angle1);
4174
4175         else
4176                 fprintf(fd2,"Depression angle to %s: %+.4f degrees\n",destination.name,angle1);
4177
4178         if (angle1!=angle2)
4179         {
4180                 if (angle2<0.0)
4181                         fprintf(fd2,"Depression");
4182                 else
4183                         fprintf(fd2,"Elevation");
4184
4185                 fprintf(fd2," angle to the first obstruction: %+.4f degrees\n",angle2);
4186         }
4187
4188         fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4189
4190         /* Receiver */
4191
4192         fprintf(fd2,"Receiver site: %s\n",destination.name);
4193
4194         if (destination.lat>=0.0)
4195         {
4196                 fprintf(fd2,"Site location: %.4f North / %.4f West",destination.lat, destination.lon);
4197                 fprintf(fd2, " (%s N / ", dec2dms(destination.lat));
4198         }
4199
4200         else
4201         {
4202                 fprintf(fd2,"Site location: %.4f South / %.4f West",-destination.lat, destination.lon);
4203                 fprintf(fd2, " (%s S / ", dec2dms(destination.lat));
4204         }
4205
4206         fprintf(fd2, "%s W)\n", dec2dms(destination.lon));
4207
4208         if (metric)
4209         {
4210                 fprintf(fd2,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(destination));
4211                 fprintf(fd2,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*destination.alt, METERS_PER_FOOT*(destination.alt+GetElevation(destination)));
4212         }
4213
4214         else
4215         {
4216                 fprintf(fd2,"Ground elevation: %.2f feet AMSL\n",GetElevation(destination));
4217                 fprintf(fd2,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",destination.alt, destination.alt+GetElevation(destination));
4218         }
4219
4220         haavt=haat(destination);
4221
4222         if (haavt>-4999.0)
4223         {
4224                 if (metric)
4225                         fprintf(fd2,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4226                 else
4227                         fprintf(fd2,"Antenna height above average terrain: %.2f feet\n",haavt);
4228         }
4229
4230         if (metric)
4231                 fprintf(fd2,"Distance to %s: %.2f kilometers\n",source.name,KM_PER_MILE*Distance(source,destination));
4232
4233         else
4234                 fprintf(fd2,"Distance to %s: %.2f miles\n",source.name,Distance(source,destination));
4235
4236         azimuth=Azimuth(destination,source);
4237
4238         angle1=ElevationAngle(destination,source);
4239         angle2=ElevationAngle2(destination,source,earthradius);
4240
4241         fprintf(fd2,"Azimuth to %s: %.2f degrees.\n",source.name,azimuth);
4242
4243         if (angle1>=0.0)
4244                 fprintf(fd2,"Elevation angle to %s: %+.4f degrees\n",source.name,angle1);
4245
4246         else
4247                 fprintf(fd2,"Depression angle to %s: %+.4f degrees\n",source.name,angle1);
4248
4249         if (angle1!=angle2)
4250         {
4251                 if (angle2<0.0)
4252                         fprintf(fd2,"Depression");
4253                 else
4254                         fprintf(fd2,"Elevation");
4255
4256                 fprintf(fd2," angle to the first obstruction: %+.4f degrees\n",angle2);
4257         }
4258
4259         fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4260
4261         fprintf(fd2,"Longley-Rice path calculation parameters used in this analysis:\n\n");
4262         fprintf(fd2,"Earth's Dielectric Constant: %.3lf\n",LR.eps_dielect);
4263         fprintf(fd2,"Earth's Conductivity: %.3lf\n",LR.sgm_conductivity);
4264         fprintf(fd2,"Atmospheric Bending Constant (N): %.3lf\n",LR.eno_ns_surfref);
4265         fprintf(fd2,"Frequency: %.3lf (MHz)\n",LR.frq_mhz);
4266         fprintf(fd2,"Radio Climate: %d (",LR.radio_climate);
4267
4268         switch (LR.radio_climate)
4269         {
4270                 case 1:
4271                 fprintf(fd2,"Equatorial");
4272                 break;
4273
4274                 case 2:
4275                 fprintf(fd2,"Continental Subtropical");
4276                 break;
4277
4278                 case 3:
4279                 fprintf(fd2,"Maritime Subtropical");
4280                 break;
4281
4282                 case 4:
4283                 fprintf(fd2,"Desert");
4284                 break;
4285
4286                 case 5:
4287                 fprintf(fd2,"Continental Temperate");
4288                 break;
4289
4290                 case 6:
4291                 fprintf(fd2,"Martitime Temperate, Over Land");
4292                 break;
4293
4294                 case 7:
4295                 fprintf(fd2,"Maritime Temperate, Over Sea");
4296                 break;
4297
4298                 default:
4299                 fprintf(fd2,"Unknown");
4300         }
4301
4302         fprintf(fd2,")\nPolarization: %d (",LR.pol);
4303
4304         if (LR.pol==0)
4305                 fprintf(fd2,"Horizontal");
4306
4307         if (LR.pol==1)
4308                 fprintf(fd2,"Vertical");
4309
4310         fprintf(fd2,")\nFraction of Situations: %.1lf%c\n",LR.conf*100.0,37);
4311         fprintf(fd2,"Fraction of Time: %.1lf%c\n",LR.rel*100.0,37);
4312         fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4313
4314         fprintf(fd2,"Analysis Results:\n\n");
4315
4316         ReadPath(source, destination);  /* source=TX, destination=RX */
4317
4318         /* Copy elevations along path into the elev_l[] array. */
4319
4320         for (x=0; x<path.length; x++)
4321                 elev_l[x+2]=path.elevation[x]*METERS_PER_FOOT;
4322
4323         if (metric)
4324                 fprintf(fd2,"Distance (km)");
4325         else
4326                 fprintf(fd2,"Distance (mi)");
4327
4328         fprintf(fd2,"\tLoss (dB)\tErrnum\tComment\n\n");
4329
4330         fd=fopen("profile.gp","w");
4331
4332         azimuth=rint(Azimuth(source,destination));
4333
4334         for (y=2; y<(path.length-1); y++)  /* path.length-1 avoids LR error */
4335         {
4336                 distance=5280.0*path.distance[y];
4337                 source_alt=four_thirds_earth+source.alt+path.elevation[0];
4338                 dest_alt=four_thirds_earth+destination.alt+path.elevation[y];
4339                 dest_alt2=dest_alt*dest_alt;
4340                 source_alt2=source_alt*source_alt;
4341
4342                 /* Calculate the cosine of the elevation of
4343                    the receiver as seen by the transmitter. */
4344
4345                 cos_xmtr_angle=((source_alt2)+(distance*distance)-(dest_alt2))/(2.0*source_alt*distance);
4346
4347                 if (got_elevation_pattern)
4348                 {
4349                         /* If an antenna elevation pattern is available, the
4350                            following code determines the elevation angle to
4351                            the first obstruction along the path. */
4352
4353                         for (x=2, block=0; x<y && block==0; x++)
4354                         {
4355                                 distance=5280.0*(path.distance[y]-path.distance[x]);
4356                                 test_alt=four_thirds_earth+path.elevation[x];
4357
4358                                 /* Calculate the cosine of the elevation
4359                                    angle of the terrain (test point)
4360                                    as seen by the transmitter. */
4361
4362                                 cos_test_angle=((source_alt2)+(distance*distance)-(test_alt*test_alt))/(2.0*source_alt*distance);
4363
4364                                 /* Compare these two angles to determine if
4365                                    an obstruction exists.  Since we're comparing
4366                                    the cosines of these angles rather than
4367                                    the angles themselves, the sense of the
4368                                    following "if" statement is reversed from
4369                                    what it would be if the angles themselves
4370                                    were compared. */
4371
4372                                 if (cos_xmtr_angle>cos_test_angle)
4373                                         block=1;
4374                         }
4375
4376                         /* At this point, we have the elevation angle
4377                            to the first obstruction (if it exists). */
4378                 }
4379
4380                 /* Determine path loss for each point along the
4381                    path using Longley-Rice's point_to_point mode
4382                    starting at x=2 (number_of_points = 1), the
4383                    shortest distance terrain can play a role in
4384                    path loss. */
4385
4386                 elev_l[0]=y-1;  /* (number of points - 1) */
4387
4388                 /* Distance between elevation samples */
4389                 elev_l[1]=METERS_PER_MILE*(path.distance[y]-path.distance[y-1]);
4390
4391                 point_to_point(elev_l, source.alt*METERS_PER_FOOT,
4392                 destination.alt*METERS_PER_FOOT, LR.eps_dielect,
4393                 LR.sgm_conductivity, LR.eno_ns_surfref, LR.frq_mhz,
4394                 LR.radio_climate, LR.pol, LR.conf, LR.rel, loss,
4395                 strmode, errnum);
4396
4397                 if (block)
4398                         elevation=((acos(cos_test_angle))/deg2rad)-90.0;
4399                 else
4400                         elevation=((acos(cos_xmtr_angle))/deg2rad)-90.0;
4401
4402                 /* Integrate the antenna's radiation
4403                    pattern into the overall path loss. */
4404
4405                 x=(int)rint(10.0*(10.0-elevation));
4406
4407                 if (x>=0 && x<=1000)
4408                 {
4409                         pattern=(double)LR.antenna_pattern[(int)azimuth][x];
4410
4411                         if (pattern!=0.0)
4412                                 patterndB=20.0*log10(pattern);
4413                 }
4414
4415                 else
4416                         patterndB=0.0;
4417
4418                 total_loss=loss-patterndB;
4419
4420                 if (metric)
4421                 {
4422                         fprintf(fd,"%f\t%f\n",KM_PER_MILE*(path.distance[path.length-1]-path.distance[y]),total_loss);
4423                         fprintf(fd2,"%7.2f\t\t%7.2f\t\t  %d\t%s\n",KM_PER_MILE*path.distance[y],total_loss, errnum, strmode);
4424                 }
4425
4426                 else
4427                 {
4428                         fprintf(fd,"%f\t%f\n",path.distance[path.length-1]-path.distance[y],total_loss);
4429                         fprintf(fd2,"%7.2f\t\t%7.2f\t\t  %d\t%s\n",path.distance[y],total_loss, errnum, strmode);
4430                 }
4431
4432                 errflag|=errnum;
4433
4434                 if (total_loss>maxloss)
4435                         maxloss=total_loss;
4436
4437                 if (total_loss<minloss)
4438                         minloss=total_loss;
4439         }
4440
4441         fclose(fd);
4442
4443         fprintf(fd2,"\nLongley-Rice path loss between %s and %s is %.2f dB.\n",source.name, destination.name, loss);
4444
4445         if (patterndB!=0.0)
4446                 fprintf(fd2,"Total path loss including TX antenna pattern is %.2f dB.\n",total_loss);
4447
4448         if (errflag)
4449         {
4450                 fprintf(fd2,"\nNotes on \"errnum\":\n\n");
4451                 fprintf(fd2,"  0: No error.  :-)\n");
4452                 fprintf(fd2,"  1: Warning!  Some parameters are nearly out of range.\n");
4453                 fprintf(fd2,"     Results should be used with caution.\n");
4454                 fprintf(fd2,"  2: Note: Default parameters have been substituted for impossible ones.\n");
4455                 fprintf(fd2,"  3: Warning!  A combination of parameters is out of range.\n");
4456                 fprintf(fd2,"     Results are probably invalid.\n");
4457                 fprintf(fd2,"  Other: Warning!  Some parameters are out of range.\n");
4458                 fprintf(fd2,"    Results are probably invalid.\n\nEnd of Report\n");
4459         }
4460
4461         fprintf(stdout,"Longley-Rice path loss between %s and %s is %.2f dB.\n",source.name, destination.name, loss);
4462
4463         if (patterndB!=0.0)
4464                 fprintf(stdout,"Total path loss including TX antenna pattern is %.2f dB.\n",total_loss);
4465
4466         fprintf(stdout,"Path Loss Report written to: \"%s\"\n",report_name);
4467         fflush(stdout);
4468
4469         fclose(fd2);
4470
4471         if (name[0]==0)
4472         {
4473                 /* Default filename and output file type */
4474
4475                 strncpy(filename,"loss\0",5);
4476                 strncpy(term,"png\0",4);
4477                 strncpy(ext,"png\0",4);
4478         }
4479
4480         else
4481         {
4482                 /* Grab extension and terminal type from "name" */
4483
4484                 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
4485                         filename[x]=name[x];
4486
4487                 if (name[x]=='.')
4488                 {
4489                         for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
4490                         {
4491                                 term[y]=tolower(name[x]);
4492                                 ext[y]=term[y];
4493                         }
4494
4495                         ext[y]=0;
4496                         term[y]=0;
4497                         filename[z]=0;
4498                 }
4499
4500                 else
4501                 {       /* No extension -- Default is png */
4502
4503                         filename[x]=0;
4504                         strncpy(term,"png\0",4);
4505                         strncpy(ext,"png\0",4);
4506                 }
4507         }
4508
4509         /* Either .ps or .postscript may be used
4510            as an extension for postscript output. */
4511
4512         if (strncmp(term,"postscript",10)==0)
4513                 strncpy(ext,"ps\0",3);
4514
4515         else if (strncmp(ext,"ps",2)==0)
4516                 strncpy(term,"postscript enhanced color\0",26);
4517
4518         fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
4519         fflush(stdout);
4520
4521         fd=fopen("splat.gp","w");
4522
4523         fprintf(fd,"set grid\n");
4524         fprintf(fd,"set yrange [%2.3f to %2.3f]\n", minloss, maxloss);
4525         fprintf(fd,"set encoding iso_8859_1\n");
4526         fprintf(fd,"set term %s\n",term);
4527         fprintf(fd,"set title \"SPLAT! Loss Profile Along Path Between %s and %s (%.2f%c azimuth)\"\n",destination.name, source.name, Azimuth(destination,source),176);
4528
4529         if (metric)
4530                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(destination,source));
4531         else
4532                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(destination,source));
4533
4534         if (got_azimuth_pattern || got_elevation_pattern)
4535                 fprintf(fd,"set ylabel \"Total Path Loss (including TX antenna pattern) (dB)");
4536         else
4537                 fprintf(fd,"set ylabel \"Longley-Rice Path Loss (dB)");
4538
4539         fprintf(fd,"\"\nset output \"%s.%s\"\n",filename,ext);
4540         fprintf(fd,"plot \"profile.gp\" title \"Path Loss\" with lines\n");
4541
4542         fclose(fd);
4543
4544         x=system("gnuplot splat.gp");
4545
4546         if (x!=-1)
4547         {
4548                 unlink("splat.gp");
4549                 unlink("profile.gp");
4550                 unlink("reference.gp");
4551
4552                 fprintf(stdout," Done!\n");
4553                 fflush(stdout);
4554         }
4555
4556         else
4557                 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
4558 }
4559
4560 void ObstructionReport(struct site xmtr, struct site rcvr, char report, double f)
4561 {
4562         int     x;
4563         struct  site site_x;
4564         double  h_r, h_t, h_x, h_r_orig, cos_tx_angle, cos_test_angle,
4565                 cos_tx_angle_f1, cos_tx_angle_fpt6, haavt, d_tx, d_x,
4566                 h_r_f1, h_r_fpt6, h_f, h_los, lambda=0.0, azimuth,
4567                 pattern, patterndB, distance, angle1, angle2;
4568         char    report_name[80], string[255], string_fpt6[255],
4569                 string_f1[255];
4570         FILE    *fd;
4571
4572         sprintf(report_name,"%s-to-%s.txt",xmtr.name,rcvr.name);
4573
4574         for (x=0; report_name[x]!=0; x++)
4575                 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
4576                         report_name[x]='_';
4577
4578         fd=fopen(report_name,"w");
4579
4580         fprintf(fd,"\n\t\t--==[ SPLAT! v%s Obstruction Report ]==--\n\n",splat_version);
4581         fprintf(fd,"Analysis of great circle path between %s and %s:\n",xmtr.name, rcvr.name);
4582         fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
4583         fprintf(fd,"Transmitter site: %s\n",xmtr.name);
4584
4585
4586         if (xmtr.lat>=0.0)
4587         {
4588                 fprintf(fd,"Site location: %.4f North / %.4f West",xmtr.lat, xmtr.lon);
4589                 fprintf(fd, " (%s N / ", dec2dms(xmtr.lat));
4590         }
4591
4592         else
4593         {
4594                 fprintf(fd,"Site location: %.4f South / %.4f West",-xmtr.lat, xmtr.lon);
4595                 fprintf(fd, " (%s S / ", dec2dms(xmtr.lat));
4596         }
4597
4598         fprintf(fd, "%s W)\n", dec2dms(xmtr.lon));
4599
4600         if (metric)
4601         {
4602                 fprintf(fd,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(xmtr));
4603                 fprintf(fd,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*xmtr.alt, METERS_PER_FOOT*(xmtr.alt+GetElevation(xmtr)));
4604         }
4605
4606         else
4607         {
4608                 fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(xmtr));
4609                 fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",xmtr.alt, xmtr.alt+GetElevation(xmtr));
4610         }
4611
4612         haavt=haat(xmtr);
4613
4614         if (haavt>-4999.0)
4615         {
4616                 if (metric)
4617                         fprintf(fd,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4618                 else
4619                         fprintf(fd,"Antenna height above average terrain: %.2f feet\n",haavt);
4620         }
4621
4622         pattern=1.0;
4623         patterndB=0.0;
4624         distance=Distance(xmtr,rcvr);
4625         azimuth=Azimuth(xmtr,rcvr);
4626         angle1=ElevationAngle(xmtr,rcvr);
4627         angle2=ElevationAngle2(xmtr,rcvr,earthradius);
4628
4629         if (got_azimuth_pattern || got_elevation_pattern)
4630         {
4631                 x=(int)rint(10.0*(10.0-angle2));
4632
4633                 if (x>=0 && x<=1000)
4634                         pattern=(double)LR.antenna_pattern[(int)rint(azimuth)][x];
4635
4636                 if (pattern!=1.0)
4637                 {
4638                         fprintf(fd,"Antenna pattern toward %s: %.3f",rcvr.name,pattern);
4639                         patterndB=20.0*log10(pattern);
4640                         fprintf(fd," (%.2f dB)\n",patterndB);
4641                 }
4642         }
4643
4644         if (metric)
4645                 fprintf(fd,"Distance to %s: %.2f kilometers\n",rcvr.name,KM_PER_MILE*distance);
4646
4647         else
4648                 fprintf(fd,"Distance to %s: %.2f miles\n",rcvr.name,distance);
4649
4650         fprintf(fd,"Azimuth to %s: %.2f degrees\n",rcvr.name,azimuth);
4651
4652         if (angle1>=0.0)
4653                 fprintf(fd,"Elevation angle to %s: %+.4f degrees\n",rcvr.name,angle1);
4654
4655         else
4656                 fprintf(fd,"Depression angle to %s: %+.4f degrees\n",rcvr.name,angle1);
4657
4658         if (angle1!=angle2)
4659         {
4660                 if (angle2<0.0)
4661                         fprintf(fd,"Depression");
4662                 else
4663                         fprintf(fd,"Elevation");
4664
4665                 fprintf(fd," angle to the first obstruction: %+.4f degrees\n",angle2);
4666         }
4667
4668         fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
4669
4670         /* Receiver */
4671
4672         fprintf(fd,"Receiver site: %s\n",rcvr.name);
4673
4674         if (rcvr.lat>=0.0)
4675         {
4676                 fprintf(fd,"Site location: %.4f North / %.4f West",rcvr.lat, rcvr.lon);
4677                 fprintf(fd, " (%s N / ", dec2dms(rcvr.lat));
4678         }
4679
4680         else
4681         {
4682                 fprintf(fd,"Site location: %.4f South / %.4f West",-rcvr.lat, rcvr.lon);
4683                 fprintf(fd, " (%s S / ", dec2dms(rcvr.lat));
4684         }
4685
4686         fprintf(fd, "%s W)\n", dec2dms(rcvr.lon));
4687
4688         if (metric)
4689         {
4690                 fprintf(fd,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(rcvr));
4691                 fprintf(fd,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*rcvr.alt, METERS_PER_FOOT*(rcvr.alt+GetElevation(rcvr)));
4692         }
4693
4694         else
4695         {
4696                 fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(rcvr));
4697                 fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",rcvr.alt, rcvr.alt+GetElevation(rcvr));
4698         }
4699
4700         haavt=haat(rcvr);
4701
4702         if (haavt>-4999.0)
4703         {
4704                 if (metric)
4705                         fprintf(fd,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4706                 else
4707                         fprintf(fd,"Antenna height above average terrain: %.2f feet\n",haavt);
4708         }
4709
4710         azimuth=Azimuth(rcvr,xmtr);
4711         angle1=ElevationAngle(rcvr,xmtr);
4712         angle2=ElevationAngle2(rcvr,xmtr,earthradius);
4713
4714         if (metric)
4715                 fprintf(fd,"Distance to %s: %.2f kilometers\n",xmtr.name,KM_PER_MILE*distance);
4716         else
4717                 fprintf(fd,"Distance to %s: %.2f miles\n",xmtr.name,distance);
4718
4719         fprintf(fd,"Azimuth to %s: %.2f degrees\n",xmtr.name,azimuth);
4720
4721         if (angle1>=0.0)
4722                 fprintf(fd,"Elevation to %s: %+.4f degrees\n",xmtr.name,angle1);
4723
4724         else
4725                 fprintf(fd,"Depression angle to %s: %+.4f degrees\n",xmtr.name,angle1);
4726
4727         if (angle1!=angle2)
4728         {
4729                 if (angle2<0.0)
4730                         fprintf(fd,"Depression");
4731                 else
4732                         fprintf(fd,"Elevation");
4733
4734                 fprintf(fd," angle to the first obstruction: %+.4f degrees\n",angle2);
4735
4736         }
4737
4738         fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
4739
4740         if (report=='y')
4741         {
4742                 /* Generate profile of the terrain.  Create the path
4743                    from transmitter to receiver because that's the
4744                    way the original los() function did it, and going
4745                    the other way can yield slightly different results. */
4746
4747                 ReadPath(xmtr,rcvr);
4748                 h_r=GetElevation(rcvr)+rcvr.alt+earthradius;
4749                 h_r_f1=h_r;
4750                 h_r_fpt6=h_r;
4751                 h_r_orig=h_r;
4752                 h_t=GetElevation(xmtr)+xmtr.alt+earthradius;
4753                 d_tx=5280.0*Distance(rcvr,xmtr);
4754                 cos_tx_angle=((h_r*h_r)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r*d_tx);
4755                 cos_tx_angle_f1=cos_tx_angle;
4756                 cos_tx_angle_fpt6=cos_tx_angle;
4757
4758                 if (f)
4759                         lambda=9.8425e8/(f*1e6);
4760
4761                 /* At each point along the path calculate the cosine
4762                    of a sort of "inverse elevation angle" at the receiver.
4763                    From the antenna, 0 deg. looks at the ground, and 90 deg.
4764                    is parallel to the ground.
4765
4766                    Start at the receiver.  If this is the lowest antenna,
4767                    then terrain obstructions will be nearest to it.  (Plus,
4768                    that's the way the original los() did it.)
4769
4770                    Calculate cosines only.  That's sufficient to compare
4771                    angles and it saves the extra computational burden of
4772                    acos().  However, note the inverted comparison: if
4773                    acos(A) > acos(B), then B > A. */
4774
4775                 for (x=path.length-1; x>0; x--)
4776                 {
4777                         site_x.lat=path.lat[x];
4778                         site_x.lon=path.lon[x];
4779                         site_x.alt=0.0;
4780
4781                         h_x=GetElevation(site_x)+earthradius;
4782                         d_x=5280.0*Distance(rcvr,site_x);
4783
4784                         /* Deal with the LOS path first. */
4785
4786                         cos_test_angle=((h_r*h_r)+(d_x*d_x)-(h_x*h_x))/(2.0*h_r*d_x);
4787
4788                         if (cos_tx_angle>cos_test_angle)
4789                         {
4790                                 if (h_r==h_r_orig)
4791                                         fprintf(fd,"SPLAT! detected obstructions at:\n\n");
4792
4793                                 if (site_x.lat>=0.0)
4794                                 {
4795                                         if (metric)
4796                                                 fprintf(fd,"\t%.4f N, %.4f W, %5.2f kilometers, %6.2f meters AMSL\n",site_x.lat, site_x.lon, KM_PER_MILE*(d_x/5280.0), METERS_PER_FOOT*(h_x-earthradius));
4797                                         else
4798                                                 fprintf(fd,"\t%.4f N, %.4f W, %5.2f miles, %6.2f feet AMSL\n",site_x.lat, site_x.lon, d_x/5280.0, h_x-earthradius);
4799                                 }
4800
4801                                 else
4802                                 {
4803                                         if (metric)
4804                                                 fprintf(fd,"\t%.4f S, %.4f W, %5.2f kilometers, %6.2f meters AMSL\n",-site_x.lat, site_x.lon, KM_PER_MILE*(d_x/5280.0), METERS_PER_FOOT*(h_x-earthradius));
4805                                         else
4806
4807                                                 fprintf(fd,"\t%.4f S, %.4f W, %5.2f miles, %6.2f feet AMSL\n",-site_x.lat, site_x.lon, d_x/5280.0, h_x-earthradius);
4808                                 }
4809                         }
4810
4811                         while (cos_tx_angle>cos_test_angle)
4812                         {
4813                                 h_r+=1;
4814                                 cos_test_angle=((h_r*h_r)+(d_x*d_x)-(h_x*h_x))/(2.0*h_r*d_x);
4815                                 cos_tx_angle=((h_r*h_r)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r*d_tx);
4816                         }
4817
4818                         if (f)
4819                         {
4820                                 /* Now clear the first Fresnel zone, but don't
4821                                    clutter the obstruction report. */
4822
4823                                 cos_tx_angle_f1=((h_r_f1*h_r_f1)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r_f1*d_tx);
4824                                 h_los=sqrt(h_r_f1*h_r_f1+d_x*d_x-2*h_r_f1*d_x*cos_tx_angle_f1);
4825                                 h_f=h_los-sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4826
4827                                 while (h_f<h_x)
4828                                 {
4829                                         h_r_f1+=1;
4830                                         cos_tx_angle_f1=((h_r_f1*h_r_f1)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r_f1*d_tx);
4831                                         h_los=sqrt(h_r_f1*h_r_f1+d_x*d_x-2*h_r_f1*d_x*cos_tx_angle_f1);
4832                                         h_f=h_los-sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4833                                 }
4834
4835                                 /* And clear the 60% F1 zone. */
4836
4837                                 cos_tx_angle_fpt6=((h_r_fpt6*h_r_fpt6)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r_fpt6*d_tx);
4838                                 h_los=sqrt(h_r_fpt6*h_r_fpt6+d_x*d_x-2*h_r_fpt6*d_x*cos_tx_angle_fpt6);
4839                                 h_f=h_los-0.6*sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4840
4841                                 while (h_f<h_x)
4842                                 {
4843                                         h_r_fpt6+=1;
4844                                         cos_tx_angle_fpt6=((h_r_fpt6*h_r_fpt6)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r_fpt6*d_tx);
4845                                         h_los=sqrt(h_r_fpt6*h_r_fpt6+d_x*d_x-2*h_r_fpt6*d_x*cos_tx_angle_fpt6);
4846                                         h_f=h_los-0.6*sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4847                                 }
4848                         }
4849                 }
4850
4851                 if (h_r>h_r_orig)
4852                 {
4853                         if (metric)
4854                                 sprintf(string,"\nAntenna at %s must be raised to at least %.2f meters AGL\nto clear all obstructions detected by SPLAT!\n",rcvr.name, METERS_PER_FOOT*(h_r-GetElevation(rcvr)-earthradius));
4855                         else
4856                                 sprintf(string,"\nAntenna at %s must be raised to at least %.2f feet AGL\nto clear all obstructions detected by SPLAT!\n",rcvr.name, h_r-GetElevation(rcvr)-earthradius);
4857                 }
4858
4859                 else
4860                         sprintf(string,"\nNo obstructions to LOS path due to terrain were detected by SPLAT!\n");
4861
4862                 if (f)
4863                 {
4864                         if (h_r_fpt6>h_r_orig)
4865                         {
4866                                 if (metric)
4867                                         sprintf(string_fpt6,"\nAntenna at %s must be raised to at least %.2f meters AGL\nto clear 60%c of the first Fresnel zone.\n",rcvr.name, METERS_PER_FOOT*(h_r_fpt6-GetElevation(rcvr)-earthradius),37);
4868
4869                                 else
4870                                         sprintf(string_fpt6,"\nAntenna at %s must be raised to at least %.2f feet AGL\nto clear 60%c of the first Fresnel zone.\n",rcvr.name, h_r_fpt6-GetElevation(rcvr)-earthradius,37);
4871                         }
4872
4873                         else
4874                                 sprintf(string_fpt6,"\n60%c of the first Fresnel zone is clear.\n",37);
4875
4876                         if (h_r_f1>h_r_orig)
4877                         {
4878                                 if (metric)
4879                                         sprintf(string_f1,"\nAntenna at %s must be raised to at least %.2f meters AGL\nto clear the first Fresnel zone.\n",rcvr.name, METERS_PER_FOOT*(h_r_f1-GetElevation(rcvr)-earthradius));
4880
4881                                 else
4882                                         sprintf(string_f1,"\nAntenna at %s must be raised to at least %.2f feet AGL\nto clear the first Fresnel zone.\n",rcvr.name, h_r_f1-GetElevation(rcvr)-earthradius);
4883
4884                         }
4885
4886                         else
4887                             sprintf(string_f1,"\nThe first Fresnel zone is clear.\n\n");
4888                 }
4889         }
4890
4891         fprintf(fd,"%s",string);
4892
4893         if (f)
4894         {
4895                 fprintf(fd,"%s",string_f1);
4896                 fprintf(fd,"%s",string_fpt6);
4897         }
4898
4899         fclose(fd);
4900
4901         /* Display report summary on terminal */
4902
4903         /* Line-of-sight status */
4904
4905         fprintf(stdout,"%s",string);
4906
4907         if (f)
4908         {
4909                 /* Fresnel zone status */
4910
4911                 fprintf(stdout,"%s",string_f1);
4912                 fprintf(stdout,"%s",string_fpt6);
4913         }
4914
4915         fprintf(stdout, "\nObstruction report written to: \"%s\"\n",report_name);
4916
4917         fflush(stdout);
4918 }
4919
4920 void SiteReport(struct site xmtr)
4921 {
4922         char    report_name[80];
4923         double  terrain;
4924         int     x, azi;
4925         FILE    *fd;
4926
4927         sprintf(report_name,"%s-site_report.txt",xmtr.name);
4928
4929         for (x=0; report_name[x]!=0; x++)
4930                 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
4931                         report_name[x]='_';
4932
4933         fd=fopen(report_name,"w");
4934
4935         fprintf(fd,"\n\t--==[ SPLAT! v%s Site Analysis Report For: %s ]==--\n\n",splat_version,xmtr.name);
4936
4937         fprintf(fd,"---------------------------------------------------------------------------\n\n");
4938
4939         if (xmtr.lat>=0.0)
4940         {
4941                 fprintf(fd,"Site location: %.4f North / %.4f West",xmtr.lat, xmtr.lon);
4942                 fprintf(fd, " (%s N / ",dec2dms(xmtr.lat));
4943         }
4944
4945         else
4946         {
4947                 fprintf(fd,"Site location: %.4f South / %.4f West",-xmtr.lat, xmtr.lon);
4948                 fprintf(fd, " (%s S / ",dec2dms(xmtr.lat));
4949         }
4950
4951         fprintf(fd, "%s W)\n",dec2dms(xmtr.lon));
4952
4953         if (metric)
4954         {
4955                 fprintf(fd,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(xmtr));
4956                 fprintf(fd,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*xmtr.alt, METERS_PER_FOOT*(xmtr.alt+GetElevation(xmtr)));
4957         }
4958
4959         else
4960         {
4961                 fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(xmtr));
4962                 fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",xmtr.alt, xmtr.alt+GetElevation(xmtr));
4963         }
4964
4965         terrain=haat(xmtr);
4966
4967         if (terrain>-4999.0)
4968         {
4969                 if (metric)
4970                         fprintf(fd,"Antenna height above average terrain: %.2f meters\n\n",METERS_PER_FOOT*terrain);
4971                 else
4972                         fprintf(fd,"Antenna height above average terrain: %.2f feet\n\n",terrain);
4973
4974                 /* Display the average terrain between 2 and 10 miles
4975                    from the transmitter site at azimuths of 0, 45, 90,
4976                    135, 180, 225, 270, and 315 degrees. */
4977
4978                 for (azi=0; azi<=315; azi+=45)
4979                 {
4980                         fprintf(fd,"Average terrain at %3d degrees azimuth: ",azi);
4981                         terrain=AverageTerrain(xmtr,(double)azi,2.0,10.0);
4982
4983                         if (terrain>-4999.0)
4984                         {
4985                                 if (metric)
4986                                         fprintf(fd,"%.2f meters AMSL\n",METERS_PER_FOOT*terrain);
4987                                 else
4988                                         fprintf(fd,"%.2f feet AMSL\n",terrain);
4989                         }
4990
4991                         else
4992                                 fprintf(fd,"No terrain\n");
4993                 }
4994         }
4995
4996         fprintf(fd,"\n---------------------------------------------------------------------------\n\n");
4997         fclose(fd);
4998         fprintf(stdout,"\nSite analysis report written to: \"%s\"\n",report_name);
4999 }
5000
5001 void LoadTopoData(int max_lon, int min_lon, int max_lat, int min_lat)
5002 {
5003         /* This function loads the SDF files required
5004            to cover the limits of the region specified. */
5005
5006         int     x, y, width, ymin, ymax;
5007
5008         width=ReduceAngle(max_lon-min_lon);
5009
5010         if ((max_lon-min_lon)<=180.0)
5011         {
5012                 for (y=0; y<=width; y++)
5013                         for (x=min_lat; x<=max_lat; x++)
5014                         {
5015                                 ymin=(int)(min_lon+(double)y);
5016
5017                                 while (ymin<0)
5018                                         ymin+=360;
5019
5020                                 while (ymin>=360)
5021                                         ymin-=360;
5022
5023                                 ymax=ymin+1;
5024
5025                                 while (ymax<0)
5026                                         ymax+=360;
5027
5028                                 while (ymax>=360)
5029                                         ymax-=360;
5030
5031                                 sprintf(string,"%d:%d:%d:%d",x, x+1, ymin, ymax);
5032                                 LoadSDF(string);
5033                         }
5034         }
5035
5036         else
5037         {
5038                 for (y=0; y<=width; y++)
5039                         for (x=min_lat; x<=max_lat; x++)
5040                         {
5041                                 ymin=max_lon+y;
5042
5043                                 while (ymin<0)
5044                                         ymin+=360;
5045
5046                                 while (ymin>=360)
5047                                         ymin-=360;
5048
5049                                 ymax=ymin+1;
5050
5051                                 while (ymax<0)
5052                                         ymax+=360;
5053
5054                                 while (ymax>=360)
5055                                         ymax-=360;
5056
5057                                 sprintf(string,"%d:%d:%d:%d",x, x+1, ymin, ymax);
5058                                 LoadSDF(string);
5059                         }
5060         }
5061 }
5062
5063 int LoadPLI(char *filename)
5064 {
5065         int     error=0, max_west, min_west, max_north, min_north;
5066         char    string[80], *pointer=NULL;
5067         double  latitude=0.0, longitude=0.0, azimuth=0.0, elevation=0.0,
5068                 loss=0.0;
5069         FILE    *fd;
5070
5071         fd=fopen(filename,"r");
5072
5073         if (fd!=NULL)
5074         {
5075                 fgets(string,78,fd);
5076                 pointer=strchr(string,';');
5077
5078                 if (pointer!=NULL)
5079                         *pointer=0;
5080
5081                 sscanf(string,"%d, %d",&max_west, &min_west);
5082
5083                 fgets(string,78,fd);
5084                 pointer=strchr(string,';');
5085
5086                 if (pointer!=NULL)
5087                         *pointer=0;
5088
5089                 sscanf(string,"%d, %d",&max_north, &min_north);
5090
5091                 fgets(string,78,fd);
5092                 pointer=strchr(string,';');
5093
5094                 if (pointer!=NULL)
5095                         *pointer=0;
5096
5097                 LoadTopoData(max_west-1, min_west, max_north-1, min_north);
5098
5099                 fprintf(stdout,"\nReading \"%s\"... ",filename);
5100                 fflush(stdout);
5101
5102                 fscanf(fd,"%lf, %lf, %lf, %lf, %lf",&latitude, &longitude, &azimuth, &elevation, &loss);
5103
5104                 while (feof(fd)==0)
5105                 {
5106                         if (loss>225.0)
5107                                 loss=225.0;
5108
5109                         if (loss<75.0)
5110                                 loss=75.0;
5111
5112                         loss-=75.0;
5113                         loss/=10.0;
5114                         loss+=1.0;
5115
5116                         if (loss<=(double)maxdB)
5117                                 OrMask(latitude,longitude,((unsigned char)(loss))<<3);
5118
5119                         fscanf(fd,"%lf, %lf, %lf, %lf, %lf",&latitude, &longitude, &azimuth, &elevation, &loss);
5120                 }
5121
5122                 fclose(fd);
5123
5124                 fprintf(stdout," Done!\n");
5125                 fflush(stdout);
5126         }
5127
5128         else
5129                 error=1;
5130
5131         return error;
5132 }
5133
5134 void WriteKML(struct site source, struct site destination)
5135 {
5136         int     x, y;
5137         char    block, report_name[80];
5138         double  distance, rx_alt, tx_alt, cos_xmtr_angle,
5139                 azimuth, cos_test_angle, test_alt;
5140         FILE    *fd=NULL;
5141
5142         ReadPath(source,destination);
5143
5144         sprintf(report_name,"%s-to-%s.kml",source.name,destination.name);
5145
5146         for (x=0; report_name[x]!=0; x++)
5147                 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
5148                         report_name[x]='_';
5149
5150         fd=fopen(report_name,"w");
5151
5152         fprintf(fd,"<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n");
5153         fprintf(fd,"<kml xmlns=\"http://earth.google.com/kml/2.0\">\n");
5154         fprintf(fd,"<!-- Generated by SPLAT! Version %s -->\n",splat_version);
5155         fprintf(fd,"<Folder>\n");
5156         fprintf(fd,"<name>SPLAT! Path</name>\n");
5157         fprintf(fd,"<open>1</open>\n");
5158         fprintf(fd,"<description>Path Between %s and %s</description>\n",source.name,destination.name);
5159
5160         fprintf(fd,"<Placemark>\n");
5161         fprintf(fd,"    <name>%s</name>\n",source.name);
5162         fprintf(fd,"    <description>\n");
5163         fprintf(fd,"       Transmit Site\n");
5164
5165         if (source.lat>=0.0)
5166                 fprintf(fd,"       <BR>%s North</BR>\n",dec2dms(source.lat));
5167         else
5168                 fprintf(fd,"       <BR>%s South</BR>\n",dec2dms(source.lat));
5169
5170         fprintf(fd,"       <BR>%s West</BR>\n",dec2dms(source.lon));
5171
5172         azimuth=Azimuth(source,destination);
5173         distance=Distance(source,destination);
5174
5175         if (metric)
5176                 fprintf(fd,"       <BR>%.2f km",distance*KM_PER_MILE);
5177         else
5178                 fprintf(fd,"       <BR>%.2f miles",distance);
5179
5180         fprintf(fd," to %s</BR>\n       <BR>toward an azimuth of %.2f%c</BR>\n",destination.name,azimuth,176);
5181
5182         fprintf(fd,"    </description>\n");
5183         fprintf(fd,"    <visibility>1</visibility>\n");
5184         fprintf(fd,"    <Style>\n");
5185         fprintf(fd,"      <IconStyle>\n");
5186         fprintf(fd,"        <Icon>\n");
5187         fprintf(fd,"          <href>root://icons/palette-5.png</href>\n");
5188         fprintf(fd,"          <x>224</x>\n");
5189         fprintf(fd,"          <y>224</y>\n");
5190         fprintf(fd,"          <w>32</w>\n");
5191         fprintf(fd,"          <h>32</h>\n");
5192         fprintf(fd,"        </Icon>\n");
5193         fprintf(fd,"      </IconStyle>\n");
5194         fprintf(fd,"    </Style>\n");
5195         fprintf(fd,"    <Point>\n");
5196         fprintf(fd,"      <extrude>1</extrude>\n");
5197         fprintf(fd,"      <altitudeMode>relativeToGround</altitudeMode>\n");
5198         fprintf(fd,"      <coordinates>%f,%f,30</coordinates>\n",(source.lon<180.0?-source.lon:360.0-source.lon),source.lat);
5199         fprintf(fd,"    </Point>\n");
5200         fprintf(fd,"</Placemark>\n");
5201
5202         fprintf(fd,"<Placemark>\n");
5203         fprintf(fd,"    <name>%s</name>\n",destination.name);
5204         fprintf(fd,"    <description>\n");
5205         fprintf(fd,"       Receive Site\n");
5206
5207         if (destination.lat>=0.0)
5208                 fprintf(fd,"       <BR>%s North</BR>\n",dec2dms(destination.lat));
5209         else
5210                 fprintf(fd,"       <BR>%s South</BR>\n",dec2dms(destination.lat));
5211
5212         fprintf(fd,"       <BR>%s West</BR>\n",dec2dms(destination.lon));
5213
5214
5215         if (metric)
5216                 fprintf(fd,"       <BR>%.2f km",distance*KM_PER_MILE);
5217         else
5218                 fprintf(fd,"       <BR>%.2f miles",distance);
5219
5220         fprintf(fd," to %s</BR>\n       <BR>toward an azimuth of %.2f%c</BR>\n",source.name,Azimuth(destination,source),176);
5221
5222         fprintf(fd,"    </description>\n");
5223         fprintf(fd,"    <visibility>1</visibility>\n");
5224         fprintf(fd,"    <Style>\n");
5225         fprintf(fd,"      <IconStyle>\n");
5226         fprintf(fd,"        <Icon>\n");
5227         fprintf(fd,"          <href>root://icons/palette-5.png</href>\n");
5228         fprintf(fd,"          <x>224</x>\n");
5229         fprintf(fd,"          <y>224</y>\n");
5230         fprintf(fd,"          <w>32</w>\n");
5231         fprintf(fd,"          <h>32</h>\n");
5232         fprintf(fd,"        </Icon>\n");
5233         fprintf(fd,"      </IconStyle>\n");
5234         fprintf(fd,"    </Style>\n");
5235         fprintf(fd,"    <Point>\n");
5236         fprintf(fd,"      <extrude>1</extrude>\n");
5237         fprintf(fd,"      <altitudeMode>relativeToGround</altitudeMode>\n");
5238         fprintf(fd,"      <coordinates>%f,%f,30</coordinates>\n",(destination.lon<180.0?-destination.lon:360.0-destination.lon),destination.lat);
5239         fprintf(fd,"    </Point>\n");
5240         fprintf(fd,"</Placemark>\n");
5241
5242         fprintf(fd,"<Placemark>\n");
5243         fprintf(fd,"<name>Point-to-Point Path</name>\n");
5244         fprintf(fd,"  <visibility>1</visibility>\n");
5245         fprintf(fd,"  <open>0</open>\n");
5246         fprintf(fd,"  <Style>\n");
5247         fprintf(fd,"    <LineStyle>\n");
5248         fprintf(fd,"      <color>7fffffff</color>\n");
5249         fprintf(fd,"    </LineStyle>\n");
5250         fprintf(fd,"    <PolyStyle>\n");
5251         fprintf(fd,"       <color>7fffffff</color>\n");
5252         fprintf(fd,"    </PolyStyle>\n");
5253         fprintf(fd,"  </Style>\n");
5254         fprintf(fd,"  <LineString>\n");
5255         fprintf(fd,"    <extrude>1</extrude>\n");
5256         fprintf(fd,"    <tessellate>1</tessellate>\n");
5257         fprintf(fd,"    <altitudeMode>relativeToGround</altitudeMode>\n");
5258         fprintf(fd,"    <coordinates>\n");
5259
5260         for (x=0; x<path.length; x++)
5261                 fprintf(fd,"      %f,%f,5\n",(path.lon[x]<180.0?-path.lon[x]:360.0-path.lon[x]),path.lat[x]);
5262
5263         fprintf(fd,"    </coordinates>\n");
5264         fprintf(fd,"   </LineString>\n");
5265         fprintf(fd,"</Placemark>\n");
5266
5267         fprintf(fd,"<Placemark>\n");
5268         fprintf(fd,"<name>Line-of-Sight Path</name>\n");
5269         fprintf(fd,"  <visibility>1</visibility>\n");
5270         fprintf(fd,"  <open>0</open>\n");
5271         fprintf(fd,"  <Style>\n");
5272         fprintf(fd,"    <LineStyle>\n");
5273         fprintf(fd,"      <color>ff00ff00</color>\n");
5274         fprintf(fd,"    </LineStyle>\n");
5275         fprintf(fd,"    <PolyStyle>\n");
5276         fprintf(fd,"       <color>7f00ff00</color>\n");
5277         fprintf(fd,"    </PolyStyle>\n");
5278         fprintf(fd,"  </Style>\n");
5279         fprintf(fd,"  <LineString>\n");
5280         fprintf(fd,"    <extrude>1</extrude>\n");
5281         fprintf(fd,"    <tessellate>1</tessellate>\n");
5282         fprintf(fd,"    <altitudeMode>relativeToGround</altitudeMode>\n");
5283         fprintf(fd,"    <coordinates>\n");
5284
5285         /* Walk across the "path", indentifying obstructions along the way */
5286
5287         for (y=0; y<path.length; y++)
5288         {
5289                 distance=5280.0*path.distance[y];
5290                 tx_alt=earthradius+source.alt+path.elevation[0];
5291                 rx_alt=earthradius+destination.alt+path.elevation[y];
5292
5293                 /* Calculate the cosine of the elevation of the
5294                    transmitter as seen at the temp rx point. */
5295
5296                 cos_xmtr_angle=((rx_alt*rx_alt)+(distance*distance)-(tx_alt*tx_alt))/(2.0*rx_alt*distance);
5297
5298                 for (x=y, block=0; x>=0 && block==0; x--)
5299                 {
5300                         distance=5280.0*(path.distance[y]-path.distance[x]);
5301                         test_alt=earthradius+path.elevation[x];
5302
5303                         cos_test_angle=((rx_alt*rx_alt)+(distance*distance)-(test_alt*test_alt))/(2.0*rx_alt*distance);
5304
5305                         /* Compare these two angles to determine if
5306                            an obstruction exists.  Since we're comparing
5307                            the cosines of these angles rather than
5308                            the angles themselves, the following "if"
5309                            statement is reversed from what it would
5310                            be if the actual angles were compared. */
5311
5312                         if (cos_xmtr_angle>cos_test_angle)
5313                                 block=1;
5314                 }
5315
5316                 if (block)
5317                         fprintf(fd,"      %f,%f,-30\n",(path.lon[y]<180.0?-path.lon[y]:360.0-path.lon[y]),path.lat[y]);
5318                 else
5319                         fprintf(fd,"      %f,%f,5\n",(path.lon[y]<180.0?-path.lon[y]:360.0-path.lon[y]),path.lat[y]);
5320         }
5321
5322         fprintf(fd,"    </coordinates>\n");
5323         fprintf(fd,"  </LineString>\n");
5324         fprintf(fd,"</Placemark>\n");
5325
5326         fprintf(fd,"    <LookAt>\n");
5327         fprintf(fd,"      <longitude>%f</longitude>\n",(source.lon<180.0?-source.lon:360.0-source.lon));
5328         fprintf(fd,"      <latitude>%f</latitude>\n",source.lat);
5329         fprintf(fd,"      <range>300.0</range>\n");
5330         fprintf(fd,"      <tilt>45.0</tilt>\n");
5331         fprintf(fd,"      <heading>%f</heading>\n",azimuth);
5332         fprintf(fd,"    </LookAt>\n");
5333
5334         fprintf(fd,"</Folder>\n");
5335         fprintf(fd,"</kml>\n");
5336
5337         fclose(fd);
5338
5339         fprintf(stdout, "KML file written to: \"%s\"\n",report_name);
5340
5341         fflush(stdout);
5342 }
5343
5344 int main(int argc, char *argv[])
5345 {
5346         int             x, y, z=0, min_lat, min_lon, max_lat, max_lon,
5347                         rxlat, rxlon, txlat, txlon, west_min, west_max,
5348                         north_min, north_max;
5349
5350         unsigned char   coverage=0, LRmap=0, ext[20], terrain_plot=0,
5351                         elevation_plot=0, height_plot=0,
5352                         longley_plot=0, cities=0, bfs=0, txsites=0,
5353                         count, report='y', norm=0, topomap=0, geo=0,
5354                         kml=0;
5355
5356         char            mapfile[255], header[80], city_file[5][255],
5357                         elevation_file[255], height_file[255],
5358                         longley_file[255], terrain_file[255],
5359                         string[255], rxfile[255], *env=NULL,
5360                         txfile[255], map=0, boundary_file[5][255],
5361                         udt_file[255], rxsite=0, plo_filename[255],
5362                         pli_filename[255], nf=0;
5363
5364         double          altitude=0.0, altitudeLR=0.0, tx_range=0.0,
5365                         rx_range=0.0, deg_range=0.0, deg_limit,
5366                         deg_range_lon, er_mult, freq=0.0;
5367
5368         struct          site tx_site[4], rx_site;
5369
5370         FILE            *fd;
5371
5372
5373         if (argc==1)
5374         {
5375                 fprintf(stdout,"\n\t\t --==[ SPLAT! v%s Available Options... ]==--\n\n",splat_version);
5376                 fprintf(stdout,"       -t txsite(s).qth (max of 4)\n");
5377                 fprintf(stdout,"       -r rxsite.qth\n");
5378                 fprintf(stdout,"       -c plot coverage of TX(s) with an RX antenna at X feet/meters AGL\n");
5379                 fprintf(stdout,"       -L plot path loss map of TX based on an RX at X feet/meters AGL\n");
5380                 fprintf(stdout,"       -s filename(s) of city/site file(s) to import (5 max)\n");
5381                 fprintf(stdout,"       -b filename(s) of cartographic boundary file(s) to import (max of 5)\n");
5382                 fprintf(stdout,"       -p filename of terrain profile graph to plot\n");
5383                 fprintf(stdout,"       -e filename of terrain elevation graph to plot\n");
5384                 fprintf(stdout,"       -h filename of terrain height graph to plot\n");
5385                 fprintf(stdout,"       -H filename of normalized terrain height graph to plot\n");
5386                 fprintf(stdout,"       -l filename of Longley-Rice graph to plot\n");
5387                 fprintf(stdout,"       -o filename of topographic map to generate (.ppm)\n");
5388                 fprintf(stdout,"       -u filename of user-defined terrain file to import\n");
5389                 fprintf(stdout,"       -d sdf file directory path (overrides path in ~/.splat_path file)\n");
5390                 fprintf(stdout,"       -n no analysis, brief report\n");
5391                 fprintf(stdout,"       -N no analysis, no report\n");
5392                 fprintf(stdout,"       -m earth radius multiplier\n");
5393                 fprintf(stdout,"       -f frequency for Fresnel zone calculation (MHz)\n");
5394                 fprintf(stdout,"       -R modify default range for -c or -L (miles/kilometers)\n");
5395                 fprintf(stdout,"      -db maximum loss contour to display on path loss maps (80-230 dB)\n");
5396                 fprintf(stdout,"      -nf do not plot Fresnel zones in height plots\n");
5397                 fprintf(stdout,"     -plo filename of path-loss output file\n");
5398                 fprintf(stdout,"     -pli filename of path-loss input file\n");
5399                 fprintf(stdout,"     -udt filename of user defined terrain input file\n");
5400                 fprintf(stdout,"     -geo generate an Xastir .geo georeference file (with .ppm output)\n");
5401                 fprintf(stdout,"     -kml generate a Google Earth .kml file (for point-to-point links)\n");
5402                 fprintf(stdout,"  -metric employ metric rather than imperial units for all user I/O\n\n");
5403
5404                 fprintf(stdout,"If that flew by too fast, consider piping the output through 'less':\n");
5405                 fprintf(stdout,"\n\tsplat | less\n\n");
5406                 fprintf(stdout,"Type 'man splat', or see the documentation for more details.\n\n");
5407                 fflush(stdout);
5408                 return 1;
5409         }
5410
5411         y=argc-1;
5412
5413         metric=0;
5414         rxfile[0]=0;
5415         txfile[0]=0;
5416         string[0]=0;
5417         mapfile[0]=0;
5418         elevation_file[0]=0;
5419         terrain_file[0]=0;
5420         sdf_path[0]=0;
5421         udt_file[0]=0;
5422         path.length=0;
5423         LR.frq_mhz=0.0;
5424         rx_site.lat=91.0;
5425         rx_site.lon=361.0;
5426         plo_filename[0]=0;
5427         pli_filename[0]=0;
5428         earthradius=EARTHRADIUS;
5429
5430         sprintf(header,"\n\t\t--==[ Welcome To SPLAT! v%s ]==--\n\n", splat_version);
5431
5432         for (x=0; x<4; x++)
5433         {
5434                 tx_site[x].lat=91.0;
5435                 tx_site[x].lon=361.0;
5436         }
5437
5438         for (x=0; x<MAXSLOTS; x++)
5439         {
5440                 dem[x].min_el=32768;
5441                 dem[x].max_el=-32768;
5442                 dem[x].min_north=90;
5443                 dem[x].max_north=-90;
5444                 dem[x].min_west=360;
5445                 dem[x].max_west=-1;
5446         }
5447
5448         /* Scan for command line arguments */
5449
5450         for (x=1; x<=y; x++)
5451         {
5452                 if (strcmp(argv[x],"-R")==0)
5453                 {
5454                         z=x+1;
5455
5456                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5457                         {
5458                                 sscanf(argv[z],"%lf",&max_range);
5459
5460                                 if (max_range<0.0)
5461                                         max_range=0.0;
5462
5463                                 if (max_range>1000.0)
5464                                         max_range=1000.0;
5465                         }
5466                 }
5467
5468                 if (strcmp(argv[x],"-m")==0)
5469                 {
5470                         z=x+1;
5471
5472                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5473                         {
5474                                 sscanf(argv[z],"%lf",&er_mult);
5475
5476                                 if (er_mult<0.1)
5477                                         er_mult=1.0;
5478
5479                                 if (er_mult>1.0e6)
5480                                         er_mult=1.0e6;
5481
5482                                 earthradius*=er_mult;
5483                         }
5484                 }
5485
5486                 if (strcmp(argv[x],"-o")==0)
5487                 {
5488                         z=x+1;
5489
5490                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5491                                 strncpy(mapfile,argv[z],253);
5492                         map=1;
5493                 }
5494
5495                 if (strcmp(argv[x],"-u")==0)
5496                 {
5497                         z=x+1;
5498
5499                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5500                                 strncpy(udt_file,argv[z],253);
5501                 }
5502
5503                 if (strcmp(argv[x],"-c")==0)
5504                 {
5505                         z=x+1;
5506
5507                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5508                         {
5509                                 sscanf(argv[z],"%lf",&altitude);
5510                                 coverage=1;
5511                         }
5512                 }
5513
5514                 if (strcmp(argv[x],"-db")==0 || strcmp(argv[x],"-dB")==0)
5515                 {
5516                         z=x+1;
5517
5518                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5519                         {
5520                                 sscanf(argv[z],"%d",&maxdB);
5521
5522                                 maxdB=abs(maxdB);
5523
5524                                 if (maxdB<80)
5525                                         maxdB=80;
5526
5527                                 if (maxdB>230)
5528                                         maxdB=230;
5529                         }
5530                 }
5531
5532                 if (strcmp(argv[x],"-p")==0)
5533                 {
5534                         z=x+1;
5535
5536                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5537                         {
5538                                 strncpy(terrain_file,argv[z],253);
5539                                 terrain_plot=1;
5540                         }
5541                 }
5542
5543                 if (strcmp(argv[x],"-e")==0)
5544                 {
5545                         z=x+1;
5546
5547                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5548                         {
5549                                 strncpy(elevation_file,argv[z],253);
5550                                 elevation_plot=1;
5551                         }
5552                 }
5553
5554                 if (strcmp(argv[x],"-h")==0 || strcmp(argv[x],"-H")==0)
5555                 {
5556                         z=x+1;
5557
5558                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5559                         {
5560                                 strncpy(height_file,argv[z],253);
5561                                 height_plot=1;
5562                         }
5563
5564                         if (strcmp(argv[x],"-H")==0)
5565                                 norm=1;
5566                         else
5567                                 norm=0;
5568                 }
5569
5570                 if (strcmp(argv[x],"-n")==0)
5571                 {
5572                         report='n';
5573                         map=1;
5574                 }
5575
5576                 if (strcmp(argv[x],"-N")==0)
5577                 {
5578                         report='N';
5579                         map=1;
5580                 }
5581
5582                 if (strcmp(argv[x],"-metric")==0)
5583                         metric=1;
5584
5585                 if (strcmp(argv[x],"-geo")==0)
5586                         geo=1;
5587
5588                 if (strcmp(argv[x],"-kml")==0)
5589                         kml=1;
5590
5591                 if (strcmp(argv[x],"-nf")==0)
5592                         nf=1;
5593
5594                 if (strcmp(argv[x],"-d")==0)
5595                 {
5596                         z=x+1;
5597
5598                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5599                                 strncpy(sdf_path,argv[z],253);
5600                 }
5601
5602                 if (strcmp(argv[x],"-t")==0)
5603                 {
5604                         /* Read Transmitter Location */
5605
5606                         z=x+1;
5607
5608                         while (z<=y && argv[z][0] && argv[z][0]!='-' && txsites<4)
5609                         {
5610                                 strncpy(txfile,argv[z],253);
5611                                 tx_site[txsites]=LoadQTH(txfile);
5612                                 txsites++;
5613                                 z++;
5614                         }
5615
5616                         z--;
5617                 }
5618
5619                 if (strcmp(argv[x],"-L")==0)
5620                 {
5621                         z=x+1;
5622
5623                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5624                         {
5625                                 sscanf(argv[z],"%lf",&altitudeLR);
5626
5627                                 if (coverage)
5628                                         fprintf(stdout,"c and L are exclusive options, ignoring L.\n");
5629                                 else
5630                                 {
5631                                         LRmap=1;
5632                                         ReadLRParm(txfile);
5633                                 }
5634                         }
5635                 }
5636
5637                 if (strcmp(argv[x],"-l")==0)
5638                 {
5639                         z=x+1;
5640
5641                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5642                         {
5643                                 strncpy(longley_file,argv[z],253);
5644                                 longley_plot=1;
5645                                 /* Doing this twice is harmless */
5646                                 ReadLRParm(txfile);
5647                         }
5648                 }
5649
5650                 if (strcmp(argv[x],"-r")==0)
5651                 {
5652                         /* Read Receiver Location */
5653
5654                         z=x+1;
5655
5656                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5657                         {
5658                                 strncpy(rxfile,argv[z],253);
5659                                 rx_site=LoadQTH(rxfile);
5660                                 rxsite=1;
5661                         }
5662                 }
5663
5664                 if (strcmp(argv[x],"-s")==0)
5665                 {
5666                         /* Read city file(s) */
5667
5668                         z=x+1;
5669
5670                         while (z<=y && argv[z][0] && argv[z][0]!='-' && cities<5)
5671                         {
5672                                 strncpy(city_file[cities],argv[z],253);
5673                                 cities++;
5674                                 z++;
5675                         }
5676
5677                         z--;
5678                 }
5679
5680                 if (strcmp(argv[x],"-b")==0)
5681                 {
5682                         /* Read Boundary File(s) */
5683
5684                         z=x+1;
5685
5686                         while (z<=y && argv[z][0] && argv[z][0]!='-' && bfs<5)
5687                         {
5688                                 strncpy(boundary_file[bfs],argv[z],253);
5689                                 bfs++;
5690                                 z++;
5691                         }
5692
5693                         z--;
5694                 }
5695
5696                 if (strcmp(argv[x],"-f")==0)
5697                 {
5698                         z=x+1;
5699
5700                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5701                         {
5702                                 sscanf(argv[z],"%lf",&freq);
5703
5704                                 if (freq<20)
5705                                         freq=20;
5706
5707                                 if (freq>20e3)
5708                                         freq=20e3;
5709                         }
5710                 }
5711
5712                 if (strcmp(argv[x],"-plo")==0)
5713                 {
5714                         z=x+1;
5715
5716                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5717                                 strncpy(plo_filename,argv[z],253);
5718                 }
5719
5720                 if (strcmp(argv[x],"-pli")==0)
5721                 {
5722                         z=x+1;
5723
5724                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5725                                 strncpy(pli_filename,argv[z],253);
5726                 }
5727         }
5728
5729         /* Perform some error checking on the arguments
5730            and switches parsed from the command-line.
5731            If an error is encountered, print a message
5732            and exit gracefully. */
5733
5734         if (txsites==0)
5735         {
5736                 fprintf(stderr,"\n%c*** ERROR: No transmitter site(s) specified!\n\n",7);
5737                 exit (-1);
5738         }
5739
5740         for (x=0, y=0; x<txsites; x++)
5741         {
5742                 if (tx_site[x].lat==91.0 && tx_site[x].lon==361.0)
5743                 {
5744                         fprintf(stderr,"\n*** ERROR: Transmitter site #%d not found!",x+1);
5745                         y++;
5746                 }
5747         }
5748
5749         if (y)
5750         {
5751                 fprintf(stderr,"%c\n\n",7);
5752                 exit (-1);
5753         }
5754
5755         if ((coverage+LRmap+pli_filename[0])==0 && rx_site.lat==91.0 && rx_site.lon==361.0)
5756         {
5757                 if (max_range!=0.0 && txsites!=0)
5758                 {
5759                         /* Plot topographic map of radius "max_range" */
5760
5761                         map=0;
5762                         topomap=1;
5763                         report='N';
5764                 }
5765
5766                 else
5767                 {
5768                         fprintf(stderr,"\n%c*** ERROR: No receiver site found or specified!\n\n",7);
5769                         exit (-1);
5770                 }
5771         }
5772
5773         /* No major errors were detected.  Whew!  :-) */
5774
5775         /* Adjust input parameters if -metric option is used */
5776
5777         if (metric)
5778         {
5779                 altitudeLR/=METERS_PER_FOOT;    /* meters --> feet */
5780                 max_range/=KM_PER_MILE;         /* kilometers --> miles */
5781                 altitude/=METERS_PER_FOOT;      /* meters --> feet */
5782         }
5783
5784         /* If no SDF path was specified on the command line (-d), check
5785            for a path specified in the $HOME/.splat_path file.  If the
5786            file is not found, then sdf_path[] remains NULL, and the
5787            current working directory is assumed to contain the SDF
5788            files. */
5789
5790         if (sdf_path[0]==0)
5791         {
5792                 env=getenv("HOME");
5793                 sprintf(string,"%s/.splat_path",env);
5794                 fd=fopen(string,"r");
5795
5796                 if (fd!=NULL)
5797                 {
5798                         fgets(string,253,fd);
5799
5800                         /* Remove <CR> and/or <LF> from string */
5801
5802                         for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0 && x<253; x++);
5803                         string[x]=0;
5804
5805                         strncpy(sdf_path,string,253);
5806
5807                         fclose(fd);
5808                 }
5809         }
5810
5811         /* Ensure a trailing '/' is present in sdf_path */
5812
5813         if (sdf_path[0])
5814         {
5815                 x=strlen(sdf_path);
5816
5817                 if (sdf_path[x-1]!='/' && x!=0)
5818                 {
5819                         sdf_path[x]='/';
5820                         sdf_path[x+1]=0;
5821                 }
5822         }
5823
5824         fprintf(stdout,"%s",header);
5825         fflush(stdout);
5826
5827         if (pli_filename[0])
5828         {
5829                 y=LoadPLI(pli_filename);
5830
5831                 for (x=0; x<txsites; x++)
5832                         PlaceMarker(tx_site[x]);
5833
5834                 if (rxsite)
5835                         PlaceMarker(rx_site);
5836
5837                 if (bfs)
5838                 {
5839                         for (x=0; x<bfs; x++)
5840                                 LoadBoundaries(boundary_file[x]);
5841                 }
5842
5843                 if (cities)
5844                 {
5845                         for (x=0; x<cities; x++)
5846                                 LoadCities(city_file[x]);
5847                 }
5848
5849                 WritePPMLR(mapfile,geo);
5850
5851                 exit(0);
5852         }
5853
5854         x=0;
5855         y=0;
5856
5857         min_lat=90;
5858         max_lat=-90;
5859
5860         min_lon=(int)floor(tx_site[0].lon);
5861         max_lon=(int)floor(tx_site[0].lon);
5862
5863         for (y=0, z=0; z<txsites; z++)
5864         {
5865                 txlat=(int)floor(tx_site[z].lat);
5866                 txlon=(int)floor(tx_site[z].lon);
5867
5868                 if (txlat<min_lat)
5869                         min_lat=txlat;
5870
5871                 if (txlat>max_lat)
5872                         max_lat=txlat;
5873
5874                 if (LonDiff(txlon,min_lon)<0.0)
5875                         min_lon=txlon;
5876
5877                 if (LonDiff(txlon,max_lon)>0.0)
5878                         max_lon=txlon;
5879         }
5880
5881         if (rxsite)
5882         {
5883                 rxlat=(int)floor(rx_site.lat);
5884                 rxlon=(int)floor(rx_site.lon);
5885
5886                 if (rxlat<min_lat)
5887                         min_lat=rxlat;
5888
5889                 if (rxlat>max_lat)
5890                         max_lat=rxlat;
5891
5892                 if (LonDiff(rxlon,min_lon)<0.0)
5893                         min_lon=rxlon;
5894
5895                 if (LonDiff(rxlon,max_lon)>0.0)
5896                         max_lon=rxlon;
5897         }
5898
5899         /* Load the required SDF files */
5900
5901         LoadTopoData(max_lon, min_lon, max_lat, min_lat);
5902
5903         if (coverage | LRmap | topomap)
5904         {
5905                 if (LRmap)
5906                         txsites=1;
5907
5908                 for (z=0; z<txsites; z++)
5909                 {
5910                         /* "Ball park" estimates used to load any additional
5911                            SDF files required to conduct this analysis. */
5912
5913                         tx_range=sqrt(1.5*(tx_site[z].alt+GetElevation(tx_site[z])));
5914
5915                         if (LRmap)
5916                                 rx_range=sqrt(1.5*altitudeLR);
5917                         else
5918                                 rx_range=sqrt(1.5*altitude);
5919
5920                         /* deg_range determines the maximum
5921                            amount of topo data we read */
5922
5923                         deg_range=(tx_range+rx_range)/69.0;
5924
5925                         /* max_range sets the maximum size of the
5926                            analysis.  A small, non-zero amount can
5927                            be used to shrink the size of the analysis
5928                            and limit the amount of topo data read by
5929                            SPLAT!  A very large number will only increase
5930                            the width of the analysis, not the size of
5931                            the map. */
5932
5933                         if (max_range==0.0)
5934                                 max_range=tx_range+rx_range;
5935
5936                         if (max_range<(tx_range+rx_range))
5937                                 deg_range=max_range/69.0;
5938
5939                         /* Prevent the demand for a really wide coverage
5940                            from allocating more slots than are available
5941                            in memory. */
5942
5943                         switch (MAXSLOTS)
5944                         {
5945                                 case 2: deg_limit=0.25;
5946                                         break;
5947
5948                                 case 4: deg_limit=0.5;
5949                                         break;
5950
5951                                 case 9: deg_limit=1.0;
5952                                         break;
5953
5954                                 case 16: deg_limit=2.0;
5955                                         break;
5956
5957                                 case 25: deg_limit=3.0;
5958                         }
5959
5960                         if (tx_site[z].lat<70.0)
5961                                 deg_range_lon=deg_range/cos(deg2rad*tx_site[z].lat);
5962                         else
5963                                 deg_range_lon=deg_range/cos(deg2rad*70.0);
5964
5965                         /* Correct for squares in degrees not being square in miles */
5966
5967                         if (deg_range>deg_limit)
5968                                 deg_range=deg_limit;
5969
5970                         if (deg_range_lon>deg_limit)
5971                                 deg_range_lon=deg_limit;
5972
5973
5974                         north_min=(int)floor(tx_site[z].lat-deg_range);
5975                         north_max=(int)floor(tx_site[z].lat+deg_range);
5976
5977                         west_min=(int)floor(tx_site[z].lon-deg_range_lon);
5978
5979                         while (west_min<0)
5980                                 west_min+=360;
5981
5982                         while (west_min>=360)
5983                                 west_min-=360;
5984
5985                         west_max=(int)floor(tx_site[z].lon+deg_range_lon);
5986
5987                         while (west_max<0)
5988                                 west_max+=360;
5989
5990                         while (west_max>=360)
5991                                 west_max-=360;
5992
5993                         if (north_min<min_lat)
5994                                 min_lat=north_min;
5995
5996                         if (north_max>max_lat)
5997                                 max_lat=north_max;
5998
5999                         if (LonDiff(west_min,min_lon)<0.0)
6000                                 min_lon=west_min;
6001
6002                         if (LonDiff(west_max,max_lon)>0.0)
6003                                 max_lon=west_max;
6004                 }
6005
6006                 /* Load any additional SDF files, if required */
6007
6008                 LoadTopoData(max_lon, min_lon, max_lat, min_lat);
6009         }
6010
6011         if (udt_file[0])
6012                 LoadUDT(udt_file);
6013
6014         if (mapfile[0] && topomap==0)
6015                 map=1;
6016
6017         if (freq==0.0 && nf==0)
6018                 freq=LR.frq_mhz;
6019
6020         else if (nf==1)
6021                 freq=0.0;
6022
6023         if (coverage | LRmap)
6024         {
6025                 for (x=0; x<txsites; x++)
6026                 {
6027                         if (coverage)
6028                                 PlotCoverage(tx_site[x],altitude);
6029
6030                         if (LRmap)
6031                                 PlotLRMap(tx_site[x],altitudeLR,plo_filename);
6032
6033                         PlaceMarker(tx_site[x]);
6034
6035                         if (report!='N')
6036                                 SiteReport(tx_site[x]);
6037                 }
6038
6039                 map=1;
6040         }
6041
6042         if (coverage==0 && LRmap==0)
6043         {
6044                 PlaceMarker(rx_site);
6045
6046                 for (x=0; x<txsites; x++)
6047                 {
6048                         PlaceMarker(tx_site[x]);
6049
6050                         if (report=='y')
6051                         {
6052                                 switch (x)
6053                                 {
6054                                         case 0:
6055                                                 PlotPath(tx_site[x],rx_site,1);
6056                                                 break;
6057
6058                                         case 1:
6059                                                 PlotPath(tx_site[x],rx_site,8);
6060                                                 break;
6061
6062                                         case 2:
6063                                                 PlotPath(tx_site[x],rx_site,16);
6064                                                 break;
6065
6066                                         case 3:
6067                                                 PlotPath(tx_site[x],rx_site,32);
6068                                 }
6069                         }
6070
6071                         if (report!='N')
6072                                 ObstructionReport(tx_site[x],rx_site,report,freq);
6073
6074                         if (kml)
6075                                 WriteKML(tx_site[x],rx_site);
6076                 }
6077         }
6078
6079         if (map | topomap)
6080         {
6081                 if (bfs)
6082                 {
6083                         for (x=0; x<bfs; x++)
6084                                 LoadBoundaries(boundary_file[x]);
6085                 }
6086
6087                 if (cities)
6088                 {
6089                         for (x=0; x<cities; x++)
6090                                 LoadCities(city_file[x]);
6091                 }
6092
6093                 if (LRmap)
6094                         WritePPMLR(mapfile,geo);
6095                 else
6096                         WritePPM(mapfile,geo);
6097         }
6098
6099         if (terrain_plot)
6100         {
6101                 if (txsites>1)
6102                 {
6103                         for (x=0; terrain_file[x]!='.' && terrain_file[x]!=0 && x<80; x++);
6104
6105                         if (terrain_file[x]=='.')  /* extension */
6106                         {
6107                                 ext[0]='.';
6108                                 for (y=1, z=x, x++; terrain_file[x]!=0 && x<253 && y<14; x++, y++)
6109                                         ext[y]=terrain_file[x];
6110
6111                                 ext[y]=0;
6112                                 terrain_file[z]=0;
6113                         }
6114
6115                         else
6116                         {
6117                                 ext[0]=0;  /* No extension */
6118                                 terrain_file[x]=0;
6119                         }
6120
6121                         for (count=0; count<txsites; count++)
6122                         {
6123                                 sprintf(string,"%s-%c%s%c",terrain_file,'1'+count,ext,0);
6124                                 GraphTerrain(tx_site[count],rx_site,string);
6125                         }
6126                 }
6127
6128                 else
6129                         GraphTerrain(tx_site[0],rx_site,terrain_file);
6130         }
6131
6132         if (elevation_plot)
6133         {
6134                 if (txsites>1)
6135                 {
6136                         for (x=0; elevation_file[x]!='.' && elevation_file[x]!=0 && x<80; x++);
6137
6138                         if (elevation_file[x]=='.')  /* extension */
6139                         {
6140                                 ext[0]='.';
6141                                 for (y=1, z=x, x++; elevation_file[x]!=0 && x<253 && y<14; x++, y++)
6142                                         ext[y]=elevation_file[x];
6143
6144                                 ext[y]=0;
6145                                 elevation_file[z]=0;
6146                         }
6147
6148                         else
6149                         {
6150                                 ext[0]=0;  /* No extension */
6151                                 elevation_file[x]=0;
6152                         }
6153
6154                         for (count=0; count<txsites; count++)
6155                         {
6156                                 sprintf(string,"%s-%c%s%c",elevation_file,'1'+count,ext,0);
6157                                 GraphElevation(tx_site[count],rx_site,string);
6158                         }
6159                 }
6160
6161                 else
6162                         GraphElevation(tx_site[0],rx_site,elevation_file);
6163         }
6164
6165         if (height_plot)
6166         {
6167                 if (txsites>1)
6168                 {
6169                         for (x=0; height_file[x]!='.' && height_file[x]!=0 && x<80; x++);
6170
6171                         if (height_file[x]=='.')  /* extension */
6172                         {
6173                                 ext[0]='.';
6174                                 for (y=1, z=x, x++; height_file[x]!=0 && x<253 && y<14; x++, y++)
6175                                         ext[y]=height_file[x];
6176
6177                                 ext[y]=0;
6178                                 height_file[z]=0;
6179                         }
6180
6181                         else
6182                         {
6183                                 ext[0]=0;  /* No extension */
6184                                 height_file[x]=0;
6185                         }
6186
6187                         for (count=0; count<txsites; count++)
6188                         {
6189                                 sprintf(string,"%s-%c%s%c",height_file,'1'+count,ext,0);
6190                                 GraphHeight(tx_site[count],rx_site,string,freq,norm);
6191                         }
6192                 }
6193
6194                 else
6195                         GraphHeight(tx_site[0],rx_site,height_file,freq,norm);
6196         }
6197
6198         if (longley_plot)
6199         {
6200                 if (txsites>1)
6201                 {
6202                         for (x=0; longley_file[x]!='.' && longley_file[x]!=0 && x<80; x++);
6203
6204                         if (longley_file[x]=='.')  /* extension */
6205                         {
6206                                 ext[0]='.';
6207                                 for (y=1, z=x, x++; longley_file[x]!=0 && x<253 && y<14; x++, y++)
6208                                         ext[y]=longley_file[x];
6209
6210                                 ext[y]=0;
6211                                 longley_file[z]=0;
6212                         }
6213
6214                         else
6215                         {
6216                                 ext[0]=0;  /* No extension */
6217                                 longley_file[x]=0;
6218                         }
6219
6220                         for (count=0; count<txsites; count++)
6221                         {
6222                                 sprintf(string,"%s-%c%s%c",longley_file,'1'+count,ext,0);
6223                                 GraphLongley(tx_site[count],rx_site,string);
6224                         }
6225                 }
6226
6227                 else
6228                         GraphLongley(tx_site[0],rx_site,longley_file);
6229         }
6230
6231         return 0;
6232 }
6233