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

   1 /****************************************************************************
   2 *      SPLAT: An RF Signal Propagation Loss and Terrain Analysis Tool       *
   3 *                         Last update: 21-Dec-2006                          *
   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.0"};
  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,
2463                 four_thirds_earth;
2464         struct  site temp;
2465
2466         ReadPath(source,destination);
2467
2468         four_thirds_earth=EARTHRADIUS*(4.0/3.0);
2469
2470         /* Copy elevations along path into the elev_l[] array. */
2471
2472         for (x=0; x<path.length; x++)
2473                 elev_l[x+2]=path.elevation[x]*METERS_PER_FOOT;
2474
2475         /* Since the only energy the Longley-Rice model considers
2476            reaching the destination is based on what is scattered
2477            or deflected from the first obstruction along the path,
2478            we first need to find the location and elevation angle
2479            of that first obstruction (if it exists).  This is done
2480            using a 4/3rds Earth radius to match the model used by
2481            Longley-Rice.  This information is required for properly
2482            integrating the antenna's elevation pattern into the
2483            calculation for overall path loss.  (Using path.length-1
2484            below avoids a Longley-Rice model error from occuring at
2485            the destination point.) */
2486
2487         for (y=2; (y<(path.length-1) && path.distance[y]<=max_range); y++)
2488         {
2489                 /* Process this point only if it
2490                    has not already been processed. */
2491
2492                 if (GetMask(path.lat[y],path.lon[y])==0)
2493                 {
2494                         distance=5280.0*path.distance[y];
2495                         source_alt=four_thirds_earth+source.alt+path.elevation[0];
2496                         dest_alt=four_thirds_earth+destination.alt+path.elevation[y];
2497                         dest_alt2=dest_alt*dest_alt;
2498                         source_alt2=source_alt*source_alt;
2499
2500                         /* Calculate the cosine of the elevation of
2501                            the receiver as seen by the transmitter. */
2502
2503                         cos_xmtr_angle=((source_alt2)+(distance*distance)-(dest_alt2))/(2.0*source_alt*distance);
2504
2505                         if (got_elevation_pattern || fd!=NULL)
2506                         {
2507                                 /* If no antenna elevation pattern is available, and
2508                                    no output file is designated, the following code
2509                                    that determines the elevation angle to the first
2510                                    obstruction along the path is bypassed. */
2511
2512                                 for (x=2, block=0; x<y && block==0; x++)
2513                                 {
2514                                         distance=5280.0*(path.distance[y]-path.distance[x]);
2515                                         test_alt=four_thirds_earth+path.elevation[x];
2516
2517                                         /* Calculate the cosine of the elevation
2518                                            angle of the terrain (test point)
2519                                            as seen by the transmitter. */
2520
2521                                         cos_test_angle=((source_alt2)+(distance*distance)-(test_alt*test_alt))/(2.0*source_alt*distance);
2522
2523                                         /* Compare these two angles to determine if
2524                                            an obstruction exists.  Since we're comparing
2525                                            the cosines of these angles rather than
2526                                            the angles themselves, the sense of the
2527                                            following "if" statement is reversed from
2528                                            what it would be if the angles themselves
2529                                            were compared. */
2530
2531                                         if (cos_xmtr_angle>cos_test_angle)
2532                                                 block=1;
2533                                 }
2534
2535                                 /* At this point, we have the elevation angle
2536                                    to the first obstruction (if it exists). */
2537                         }
2538
2539                         /* Determine attenuation for each point along the
2540                            path using Longley-Rice's point_to_point mode
2541                            starting at y=2 (number_of_points = 1), the
2542                            shortest distance terrain can play a role in
2543                            path loss. */
2544
2545                         elev_l[0]=y-1;  /* (number of points - 1) */
2546
2547                         /* Distance between elevation samples */
2548                         elev_l[1]=METERS_PER_MILE*(path.distance[y]-path.distance[y-1]);
2549
2550                         point_to_point(elev_l,source.alt*METERS_PER_FOOT,
2551                         destination.alt*METERS_PER_FOOT, LR.eps_dielect,
2552                         LR.sgm_conductivity, LR.eno_ns_surfref, LR.frq_mhz,
2553                         LR.radio_climate, LR.pol, LR.conf, LR.rel, loss,
2554                         strmode, errnum);
2555
2556                         if (block)
2557                                 elevation=((acos(cos_test_angle))/deg2rad)-90.0;
2558
2559                         else
2560                                 elevation=((acos(cos_xmtr_angle))/deg2rad)-90.0;
2561
2562                         temp.lat=path.lat[y];
2563                         temp.lon=path.lon[y];
2564
2565                         azimuth=(Azimuth(source,temp));
2566
2567                         if (fd!=NULL)
2568                         {
2569                                 /* Write path loss data to output file */
2570
2571                                 fprintf(fd,"%.7f, %.7f, %.3f, %.3f, %.2f\n",path.lat[y], path.lon[y], azimuth, elevation, loss);
2572                         }
2573
2574                         /* Integrate the antenna's radiation
2575                            pattern into the overall path loss. */
2576
2577                         x=(int)rint(10.0*(10.0-elevation));
2578
2579                         if (x>=0 && x<=1000)
2580                         {
2581                                 azimuth=rint(azimuth);
2582
2583                                 pattern=(double)LR.antenna_pattern[(int)azimuth][x];
2584
2585                                 if (pattern!=0.0)
2586                                 {
2587                                         pattern=20.0*log10(pattern);
2588                                         loss-=pattern;
2589                                 }
2590                         }
2591
2592                         if (loss>225.0)
2593                                 loss=225.0;
2594
2595                         if (loss<75.0)
2596                                 loss=75.0;
2597
2598                         loss-=75.0;
2599                         loss/=10.0;
2600                         loss+=1.0;
2601
2602                         OrMask(path.lat[y],path.lon[y],((unsigned char)(loss))<<3);
2603                 }
2604
2605                 else if (GetMask(path.lat[y],path.lon[y])==0 && path.distance[y]>max_range)
2606                         OrMask(path.lat[y],path.lon[y],1);
2607         }
2608 }
2609
2610 void PlotCoverage(struct site source, double altitude)
2611 {
2612         /* This function performs a 360 degree sweep around the
2613            transmitter site (source location), and plots the
2614            line-of-sight coverage of the transmitter on the SPLAT!
2615            generated topographic map based on a receiver located
2616            at the specified altitude (in feet AGL).  Results
2617            are stored in memory, and written out in the form
2618            of a topographic map when the WritePPM() function
2619            is later invoked. */
2620
2621         float lat, lon, one_pixel;
2622         static unsigned char mask_value;
2623         int z, count;
2624         struct site edge;
2625         unsigned char symbol[4], x;
2626
2627         /* Initialize mask_value */
2628
2629         if (mask_value!=8 && mask_value!=16 && mask_value!=32)
2630                 mask_value=1;
2631
2632         one_pixel=1.0/1200.0;
2633
2634         symbol[0]='.';
2635         symbol[1]='o';
2636         symbol[2]='O';
2637         symbol[3]='o';
2638
2639         count=0;
2640
2641         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);
2642         fflush(stdout);
2643
2644         /* 18.75=1200 pixels/degree divided by 64 loops
2645            per progress indicator symbol (.oOo) printed. */
2646
2647         z=(int)(18.75*ReduceAngle(max_west-min_west));
2648
2649         for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2650         {
2651                 if (lon>=360.0)
2652                         lon-=360.0;
2653
2654                 edge.lat=max_north;
2655                 edge.lon=lon;
2656                 edge.alt=altitude;
2657
2658                 PlotPath(source,edge,mask_value);
2659                 count++;
2660
2661                 if (count==z)
2662                 {
2663                         fprintf(stdout,"%c",symbol[x]);
2664                         fflush(stdout);
2665                         count=0;
2666
2667                         if (x==3)
2668                                 x=0;
2669                         else
2670                                 x++;
2671                 }
2672         }
2673
2674         count=0;
2675         fprintf(stdout,"\n25%c to  50%c ",37,37);
2676         fflush(stdout);
2677
2678         z=(int)(18.75*(max_north-min_north));
2679
2680         for (lat=max_north, x=0; lat>=min_north; lat-=one_pixel)
2681         {
2682                 edge.lat=lat;
2683                 edge.lon=min_west;
2684                 edge.alt=altitude;
2685
2686                 PlotPath(source,edge,mask_value);
2687                 count++;
2688
2689                 if (count==z)
2690                 {
2691                         fprintf(stdout,"%c",symbol[x]);
2692                         fflush(stdout);
2693                         count=0;
2694
2695                         if (x==3)
2696                                 x=0;
2697                         else
2698                                 x++;
2699                 }
2700         }
2701
2702         count=0;
2703         fprintf(stdout,"\n50%c to  75%c ",37,37);
2704         fflush(stdout);
2705
2706         z=(int)(18.75*ReduceAngle(max_west-min_west));
2707
2708         for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2709         {
2710                 if (lon>=360.0)
2711                         lon-=360.0;
2712
2713                 edge.lat=min_north;
2714                 edge.lon=lon;
2715                 edge.alt=altitude;
2716
2717                 PlotPath(source,edge,mask_value);
2718                 count++;
2719
2720                 if (count==z)
2721                 {
2722                         fprintf(stdout,"%c",symbol[x]);
2723                         fflush(stdout);
2724                         count=0;
2725
2726                         if (x==3)
2727                                 x=0;
2728                         else
2729                                 x++;
2730                 }
2731         }
2732
2733         count=0;
2734         fprintf(stdout,"\n75%c to 100%c ",37,37);
2735         fflush(stdout);
2736
2737         z=(int)(18.75*(max_north-min_north));
2738
2739         for (lat=min_north, x=0; lat<=max_north; lat+=one_pixel)
2740         {
2741                 edge.lat=lat;
2742                 edge.lon=max_west;
2743                 edge.alt=altitude;
2744
2745                 PlotPath(source,edge,mask_value);
2746                 count++;
2747
2748                 if (count==z)
2749                 {
2750                         fprintf(stdout,"%c",symbol[x]);
2751                         fflush(stdout);
2752                         count=0;
2753
2754                         if (x==3)
2755                                 x=0;
2756                         else
2757                                 x++;
2758                 }
2759         }
2760
2761         fprintf(stdout,"\nDone!\n");
2762         fflush(stdout);
2763
2764         /* Assign next mask value */
2765
2766         switch (mask_value)
2767         {
2768                 case 1:
2769                         mask_value=8;
2770                         break;
2771
2772                 case 8:
2773                         mask_value=16;
2774                         break;
2775
2776                 case 16:
2777                         mask_value=32;
2778         }
2779 }
2780
2781 void PlotLRMap(struct site source, double altitude, char *plo_filename)
2782 {
2783         /* This function performs a 360 degree sweep around the
2784            transmitter site (source location), and plots the
2785            Longley-Rice attenuation on the SPLAT! generated
2786            topographic map based on a receiver located at
2787            the specified altitude (in feet AGL).  Results
2788            are stored in memory, and written out in the form
2789            of a topographic map when the WritePPMLR() function
2790            is later invoked. */
2791
2792         int z, count;
2793         struct site edge;
2794         float lat, lon, one_pixel;
2795         unsigned char symbol[4], x;
2796         FILE *fd=NULL;
2797
2798         one_pixel=1.0/1200.0;
2799
2800         symbol[0]='.';
2801         symbol[1]='o';
2802         symbol[2]='O';
2803         symbol[3]='o';
2804
2805         count=0;
2806
2807         fprintf(stdout,"\nComputing Longley-Rice coverage of %s ", source.name);
2808
2809         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);
2810         fflush(stdout);
2811
2812         if (plo_filename[0]!=0)
2813                 fd=fopen(plo_filename,"wb");
2814
2815         if (fd!=NULL)
2816         {
2817                 /* Write header information to output file */
2818
2819                 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);
2820         }
2821
2822         /* 18.75=1200 pixels/degree divided by 64 loops
2823            per progress indicator symbol (.oOo) printed. */
2824
2825         z=(int)(18.75*ReduceAngle(max_west-min_west));
2826
2827         for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2828         {
2829                 if (lon>=360.0)
2830                         lon-=360.0;
2831
2832                 edge.lat=max_north;
2833                 edge.lon=lon;
2834                 edge.alt=altitude;
2835
2836                 PlotLRPath(source,edge,fd);
2837                 count++;
2838
2839                 if (count==z)
2840                 {
2841                         fprintf(stdout,"%c",symbol[x]);
2842                         fflush(stdout);
2843                         count=0;
2844
2845                         if (x==3)
2846                                 x=0;
2847                         else
2848                                 x++;
2849                 }
2850         }
2851
2852         count=0;
2853         fprintf(stdout,"\n25%c to  50%c ",37,37);
2854         fflush(stdout);
2855
2856         z=(int)(18.75*(max_north-min_north));
2857
2858         for (lat=max_north, x=0; lat>=min_north; lat-=one_pixel)
2859         {
2860                 edge.lat=lat;
2861                 edge.lon=min_west;
2862                 edge.alt=altitude;
2863
2864                 PlotLRPath(source,edge,fd);
2865                 count++;
2866
2867                 if (count==z)
2868                 {
2869                         fprintf(stdout,"%c",symbol[x]);
2870                         fflush(stdout);
2871                         count=0;
2872
2873                         if (x==3)
2874                                 x=0;
2875                         else
2876                                 x++;
2877                 }
2878         }
2879
2880         count=0;
2881         fprintf(stdout,"\n50%c to  75%c ",37,37);
2882         fflush(stdout);
2883
2884         z=(int)(18.75*ReduceAngle(max_west-min_west));
2885
2886         for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2887         {
2888                 if (lon>=360.0)
2889                         lon-=360.0;
2890
2891                 edge.lat=min_north;
2892                 edge.lon=lon;
2893                 edge.alt=altitude;
2894
2895                 PlotLRPath(source,edge,fd);
2896                 count++;
2897
2898                 if (count==z)
2899                 {
2900                         fprintf(stdout,"%c",symbol[x]);
2901                         fflush(stdout);
2902                         count=0;
2903
2904                         if (x==3)
2905                                 x=0;
2906                         else
2907                                 x++;
2908                 }
2909         }
2910
2911         count=0;
2912         fprintf(stdout,"\n75%c to 100%c ",37,37);
2913         fflush(stdout);
2914
2915         z=(int)(18.75*(max_north-min_north));
2916
2917         for (lat=min_north, x=0; lat<=max_north; lat+=one_pixel)
2918         {
2919                 edge.lat=lat;
2920                 edge.lon=max_west;
2921                 edge.alt=altitude;
2922
2923                 PlotLRPath(source,edge,fd);
2924                 count++;
2925
2926                 if (count==z)
2927                 {
2928                         fprintf(stdout,"%c",symbol[x]);
2929                         fflush(stdout);
2930                         count=0;
2931
2932                         if (x==3)
2933                                 x=0;
2934                         else
2935                                 x++;
2936                 }
2937         }
2938
2939         if (fd!=NULL)
2940                 fclose(fd);
2941
2942         fprintf(stdout,"\nDone!\n");
2943         fflush(stdout);
2944 }
2945
2946 void WritePPM(char *filename, unsigned char geo)
2947 {
2948         /* This function generates a topographic map in Portable Pix Map
2949            (PPM) format based on logarithmically scaled topology data,
2950            as well as the content of flags held in the mask[][] array.
2951            The image created is rotated counter-clockwise 90 degrees
2952            from its representation in dem[][] so that north points
2953            up and east points right in the image generated. */
2954
2955         char mapfile[255], geofile[255];
2956         unsigned char found, mask;
2957         unsigned width, height, output;
2958         int indx, x, y, x0=0, y0=0, minlat, minlon;
2959         float lat, lon, one_pixel, conversion, one_over_gamma;
2960         FILE *fd;
2961
2962         one_pixel=1.0/1200.0;
2963         one_over_gamma=1.0/GAMMA;
2964         conversion=255.0/pow((double)(max_elevation-min_elevation),one_over_gamma);
2965
2966         width=(unsigned)(1200*ReduceAngle(max_west-min_west));
2967         height=(unsigned)(1200*ReduceAngle(max_north-min_north));
2968
2969         if (filename[0]==0)
2970                 strncpy(mapfile, "map.ppm\0",8);
2971         else
2972         {
2973                 for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
2974                 {
2975                         mapfile[x]=filename[x];
2976                         geofile[x]=filename[x];
2977                 }
2978
2979                 mapfile[x]='.';
2980                 geofile[x]='.';
2981                 mapfile[x+1]='p';
2982                 geofile[x+1]='g';
2983                 mapfile[x+2]='p';
2984                 geofile[x+2]='e';
2985                 mapfile[x+3]='m';
2986                 geofile[x+3]='o';
2987                 mapfile[x+4]=0;
2988                 geofile[x+4]=0;
2989         }
2990
2991         if (geo)
2992         {
2993                 fd=fopen(geofile,"wb");
2994
2995                 fprintf(fd,"FILENAME\t%s\n",mapfile);
2996                 fprintf(fd,"#\t\tX\tY\tLong\t\tLat\n");
2997                 fprintf(fd,"TIEPOINT\t0\t0\t%d.000\t\t%d.000\n",(max_west<180?-max_west:360-max_west),max_north);
2998                 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);
2999                 fprintf(fd,"IMAGESIZE\t%u\t%u\n",width,height);
3000                 fprintf(fd,"#\n# Auto Generated by SPLAT! v%s\n#\n",splat_version);
3001
3002                 fclose(fd);
3003         }
3004
3005         fd=fopen(mapfile,"wb");
3006
3007         fprintf(fd,"P6\n%u %u\n255\n",width,height);
3008
3009         fprintf(stdout,"\nWriting \"%s\" (%ux%u pixmap image)... ",mapfile,width,height);
3010         fflush(stdout);
3011
3012         for (y=0, lat=((double)max_north)-one_pixel; y<(int)height; y++, lat-=one_pixel)
3013         {
3014                 minlat=(int)floor(lat);
3015
3016                 for (x=0, lon=((double)max_west)-one_pixel; x<(int)width; x++, lon-=one_pixel)
3017                 {
3018                         if (lon<0.0)
3019                                 lon+=360.0;
3020
3021                         minlon=(int)floor(lon);
3022
3023                         for (indx=0, found=0; indx<MAXSLOTS && found==0;)
3024                                 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
3025                                         found=1;
3026                                 else
3027                                         indx++;
3028
3029                         if (found)
3030                         {
3031                                 x0=(int)(1199.0*(lat-floor(lat)));
3032                                 y0=(int)(1199.0*(lon-floor(lon)));
3033
3034                                 mask=dem[indx].mask[x0][y0];
3035
3036                                 if (mask&2)
3037                                         /* Text Labels: Red */
3038                                         fprintf(fd,"%c%c%c",255,0,0);
3039
3040                                 else if (mask&4)
3041                                         /* County Boundaries: Light Cyan */
3042                                         fprintf(fd,"%c%c%c",128,128,255);
3043
3044                                 else switch (mask&57)
3045                                 {
3046                                         case 1:
3047                                         /* TX1: Green */
3048                                         fprintf(fd,"%c%c%c",0,255,0);
3049                                         break;
3050
3051                                         case 8:
3052                                         /* TX2: Cyan */
3053                                         fprintf(fd,"%c%c%c",0,255,255);
3054                                         break;
3055
3056                                         case 9:
3057                                         /* TX1 + TX2: Yellow */
3058                                         fprintf(fd,"%c%c%c",255,255,0);
3059                                         break;
3060
3061                                         case 16:
3062                                         /* TX3: Medium Violet */
3063                                         fprintf(fd,"%c%c%c",147,112,219);
3064                                         break;
3065
3066                                         case 17:
3067                                         /* TX1 + TX3: Pink */
3068                                         fprintf(fd,"%c%c%c",255,192,203);
3069                                         break;
3070
3071                                         case 24:
3072                                         /* TX2 + TX3: Orange */
3073                                         fprintf(fd,"%c%c%c",255,165,0);
3074                                         break;
3075
3076                                         case 25:
3077                                         /* TX1 + TX2 + TX3: Dark Green */
3078                                         fprintf(fd,"%c%c%c",0,100,0);
3079                                         break;
3080
3081                                         case 32:
3082                                         /* TX4: Sienna 1 */
3083                                         fprintf(fd,"%c%c%c",255,130,71);
3084                                         break;
3085
3086                                         case 33:
3087                                         /* TX1 + TX4: Green Yellow */
3088                                         fprintf(fd,"%c%c%c",173,255,47);
3089                                         break;
3090
3091                                         case 40:
3092                                         /* TX2 + TX4: Dark Sea Green 1 */
3093                                         fprintf(fd,"%c%c%c",193,255,193);
3094                                         break;
3095
3096                                         case 41:
3097                                         /* TX1 + TX2 + TX4: Blanched Almond */
3098                                         fprintf(fd,"%c%c%c",255,235,205);
3099                                         break;
3100
3101                                         case 48:
3102                                         /* TX3 + TX4: Dark Turquoise */
3103                                         fprintf(fd,"%c%c%c",0,206,209);
3104                                         break;
3105
3106                                         case 49:
3107                                         /* TX1 + TX3 + TX4: Medium Spring Green */
3108                                         fprintf(fd,"%c%c%c",0,250,154);
3109                                         break;
3110
3111                                         case 56:
3112                                         /* TX2 + TX3 + TX4: Tan */
3113                                         fprintf(fd,"%c%c%c",210,180,140);
3114                                         break;
3115
3116                                         case 57:
3117                                         /* TX1 + TX2 + TX3 + TX4: Gold2 */
3118                                         fprintf(fd,"%c%c%c",238,201,0);
3119                                         break;
3120
3121                                         default:
3122                                         /* Water: Medium Blue */
3123                                         if (dem[indx].data[x0][y0]==0)
3124                                                 fprintf(fd,"%c%c%c",0,0,170);
3125                                         else
3126                                         {
3127                                                 /* Elevation: Greyscale */
3128                                                 output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
3129                                                 fprintf(fd,"%c%c%c",output,output,output);
3130                                         }
3131                                 }
3132                         }
3133
3134                         else
3135                         {
3136                                 /* We should never get here, but if */
3137                                 /* we do, display the region as black */
3138
3139                                 fprintf(fd,"%c%c%c",0,0,0);
3140                         }
3141                 }
3142         }
3143
3144         fclose(fd);
3145         fprintf(stdout,"Done!\n");
3146         fflush(stdout);
3147 }
3148
3149 void WritePPMLR(char *filename, unsigned char geo)
3150 {
3151         /* This function generates a topographic map in Portable Pix Map
3152            (PPM) format based on the content of flags held in the mask[][]
3153            array (only).  The image created is rotated counter-clockwise
3154            90 degrees from its representation in dem[][] so that north
3155            points up and east points right in the image generated. */
3156
3157         char mapfile[255], geofile[255];
3158         unsigned width, height, output;
3159         unsigned char found, mask, cityorcounty;
3160         int indx, x, y, t, t2, x0, y0, minlat, minlon, loss;
3161         float lat, lon, one_pixel, conversion, one_over_gamma;
3162         FILE *fd;
3163
3164         one_pixel=1.0/1200.0;
3165         one_over_gamma=1.0/GAMMA;
3166         conversion=255.0/pow((double)(max_elevation-min_elevation),one_over_gamma);
3167
3168         width=(unsigned)(1200*ReduceAngle(max_west-min_west));
3169         height=(unsigned)(1200*ReduceAngle(max_north-min_north));
3170
3171         if (filename[0]==0)
3172                 strncpy(mapfile, "map.ppm\0",8);
3173         else
3174         {
3175                 for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
3176                 {
3177                         mapfile[x]=filename[x];
3178                         geofile[x]=filename[x];
3179                 }
3180
3181                 mapfile[x]='.';
3182                 geofile[x]='.';
3183                 mapfile[x+1]='p';
3184                 geofile[x+1]='g';
3185                 mapfile[x+2]='p';
3186                 geofile[x+2]='e';
3187                 mapfile[x+3]='m';
3188                 geofile[x+3]='o';
3189                 mapfile[x+4]=0;
3190                 geofile[x+4]=0;
3191         }
3192
3193         if (geo)
3194         {
3195                 fd=fopen(geofile,"wb");
3196
3197                 fprintf(fd,"FILENAME\t%s\n",mapfile);
3198                 fprintf(fd,"#\t\tX\tY\tLong\t\tLat\n");
3199                 fprintf(fd,"TIEPOINT\t0\t0\t%d.000\t\t%d.000\n",(max_west<180?-max_west:360-max_west),max_north);
3200                 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));
3201                 fprintf(fd,"IMAGESIZE\t%u\t%u\n",width,height+30);
3202                 fprintf(fd,"#\n# Auto Generated by SPLAT! v%s\n#\n",splat_version);
3203
3204                 fclose(fd);
3205         }
3206
3207         fd=fopen(mapfile,"wb");
3208
3209         fprintf(fd,"P6\n%u %u\n255\n",width,height+30);
3210
3211         fprintf(stdout,"\nWriting \"%s\" (%ux%u pixmap image)... ",mapfile,width,height+30);
3212         fflush(stdout);
3213
3214         for (y=0, lat=((double)max_north)-one_pixel; y<(int)height; y++, lat-=one_pixel)
3215         {
3216                 minlat=(int)floor(lat);
3217
3218                 for (x=0, lon=((double)max_west)-one_pixel; x<(int)width; x++, lon-=one_pixel)
3219                 {
3220                         if (lon<0.0)
3221                                 lon+=360.0;
3222
3223                         minlon=(int)floor(lon);
3224
3225                         for (indx=0, found=0; indx<MAXSLOTS && found==0;)
3226                                 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
3227                                         found=1;
3228                                 else
3229                                         indx++;
3230                         if (found)
3231                         {
3232                                 x0=(int)(1199.0*(lat-floor(lat)));
3233                                 y0=(int)(1199.0*(lon-floor(lon)));
3234
3235                                 mask=dem[indx].mask[x0][y0];
3236                                 loss=70+(10*(int)((mask&248)>>3));
3237                                 cityorcounty=0;
3238
3239                                 if (mask&2)
3240                                 {
3241                                         /* Text Labels - Black or Red */
3242
3243                                         if ((mask&120) && (loss<=90))
3244                                                 fprintf(fd,"%c%c%c",0,0,0);
3245                                         else
3246                                                 fprintf(fd,"%c%c%c",255,0,0);
3247
3248                                         cityorcounty=1;
3249                                 }
3250
3251                                 else if (mask&4)
3252                                 {
3253                                         /* County Boundaries: Black */
3254
3255                                         fprintf(fd,"%c%c%c",0,0,0);
3256
3257                                         cityorcounty=1;
3258                                 }
3259
3260                                 if (cityorcounty==0)
3261                                 {
3262                                         if (loss>maxdB)
3263
3264                                         { /* Display land or sea elevation */
3265
3266                                                 if (dem[indx].data[x0][y0]==0)
3267                                                         fprintf(fd,"%c%c%c",0,0,170);
3268                                                 else
3269                                                 {
3270                                                         /* Elevation: Greyscale */
3271                                                         output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
3272                                                         fprintf(fd,"%c%c%c",output,output,output);
3273                                                 }
3274                                         }
3275
3276                                         else switch (loss)
3277                                         {
3278                                                 /* Plot signal loss in color */
3279
3280                                                 case 80:
3281                                                 fprintf(fd,"%c%c%c",255,0,0);
3282                                                 break;
3283
3284                                                 case 90:
3285                                                 fprintf(fd,"%c%c%c",255,128,0);
3286                                                 break;
3287
3288                                                 case 100:
3289                                                 fprintf(fd,"%c%c%c",255,165,0);
3290                                                 break;
3291
3292                                                 case 110:
3293                                                 fprintf(fd,"%c%c%c",255,206,0);
3294                                                 break;
3295
3296                                                 case 120:
3297                                                 fprintf(fd,"%c%c%c",255,255,0);
3298                                                 break;
3299
3300                                                 case 130:
3301                                                 fprintf(fd,"%c%c%c",184,255,0);
3302                                                 break;
3303
3304                                                 case 140:
3305                                                 fprintf(fd,"%c%c%c",0,255,0);
3306                                                 break;
3307
3308                                                 case 150:
3309                                                 fprintf(fd,"%c%c%c",0,208,0);
3310                                                 break;
3311
3312                                                 case 160:
3313                                                 fprintf(fd,"%c%c%c",0,196,196);
3314                                                 break;
3315
3316                                                 case 170:
3317                                                 fprintf(fd,"%c%c%c",0,148,255);
3318                                                 break;
3319
3320                                                 case 180:
3321                                                 fprintf(fd,"%c%c%c",80,80,255);
3322                                                 break;
3323
3324                                                 case 190:
3325                                                 fprintf(fd,"%c%c%c",0,38,255);
3326                                                 break;
3327
3328                                                 case 200:
3329                                                 fprintf(fd,"%c%c%c",142,63,255);
3330                                                 break;
3331
3332                                                 case 210:
3333                                                 fprintf(fd,"%c%c%c",196,54,255);
3334                                                 break;
3335
3336                                                 case 220:
3337                                                 fprintf(fd,"%c%c%c",255,0,255);
3338                                                 break;
3339
3340                                                 case 230:
3341                                                 fprintf(fd,"%c%c%c",255,194,204);
3342                                                 break;
3343
3344                                                 default:
3345
3346                                                 if (dem[indx].data[x0][y0]==0)
3347                                                         fprintf(fd,"%c%c%c",0,0,170);
3348                                                 else
3349                                                 {
3350                                                         /* Elevation: Greyscale */
3351                                                         output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
3352                                                         fprintf(fd,"%c%c%c",output,output,output);
3353                                                 }
3354                                         }
3355                                 }
3356                         }
3357
3358                         else
3359                         {
3360                                 /* We should never get here, but if */
3361                                 /* we do, display the region as black */
3362
3363                                 fprintf(fd,"%c%c%c",0,0,0);
3364                         }
3365                 }
3366         }
3367
3368         /* Display legend along bottom of image */
3369
3370         x0=width/16;
3371
3372         for (y0=0; y0<30; y0++)
3373         {
3374                 for (indx=0; indx<16; indx++)
3375                 {
3376                         for (x=0; x<x0; x++)
3377                         {
3378                                 t=indx;
3379                                 t2=indx+8;
3380
3381                                 if (y0>=10 && y0<=18)
3382                                 {
3383                                         if (t2>9)
3384                                         {
3385                                                 if (x>=11 && x<=17)
3386                                                         if (smallfont[t2/10][y0-10][x-11])
3387                                                                 t=255;
3388                                         }
3389
3390                                         if (x>=19 && x<=25)
3391                                                 if (smallfont[t2%10][y0-10][x-19])
3392                                                         t=255;
3393
3394                                         if (x>=27 && x<=33)
3395                                                 if (smallfont[0][y0-10][x-27])
3396                                                         t=255;
3397                                 }
3398
3399                                 switch (t)
3400                                 {
3401                                         case 0:
3402                                         fprintf(fd,"%c%c%c",255,0,0);
3403                                         break;
3404
3405                                         case 1:
3406                                         fprintf(fd,"%c%c%c",255,128,0);
3407                                         break;
3408
3409                                         case 2:
3410                                         fprintf(fd,"%c%c%c",255,165,0);
3411                                         break;
3412
3413                                         case 3:
3414                                         fprintf(fd,"%c%c%c",255,206,0);
3415                                         break;
3416
3417                                         case 4:
3418                                         fprintf(fd,"%c%c%c",255,255,0);
3419                                         break;
3420
3421                                         case 5:
3422                                         fprintf(fd,"%c%c%c",184,255,0);
3423                                         break;
3424
3425                                         case 6:
3426                                         fprintf(fd,"%c%c%c",0,255,0);
3427                                         break;
3428
3429                                         case 7:
3430                                         fprintf(fd,"%c%c%c",0,208,0);
3431                                         break;
3432
3433                                         case 8:
3434                                         fprintf(fd,"%c%c%c",0,196,196);
3435                                         break;
3436
3437                                         case 9:
3438                                         fprintf(fd,"%c%c%c",0,148,255);
3439                                         break;
3440
3441                                         case 10:
3442                                         fprintf(fd,"%c%c%c",80,80,255);
3443                                         break;
3444
3445                                         case 11:
3446                                         fprintf(fd,"%c%c%c",0,38,255);
3447                                         break;
3448
3449                                         case 12:
3450                                         fprintf(fd,"%c%c%c",142,63,255);
3451                                         break;
3452
3453                                         case 13:
3454                                         fprintf(fd,"%c%c%c",196,54,255);
3455                                         break;
3456
3457                                         case 14:
3458                                         fprintf(fd,"%c%c%c",255,0,255);
3459                                         break;
3460
3461                                         case 255:
3462                                         /* Black */
3463                                         fprintf(fd,"%c%c%c",0,0,0);
3464                                         break;
3465
3466                                         default:
3467                                         fprintf(fd,"%c%c%c",255,194,204);
3468                                 }
3469                         }
3470                 }
3471         }
3472
3473         fclose(fd);
3474         fprintf(stdout,"Done!\n");
3475         fflush(stdout);
3476 }
3477
3478 void GraphTerrain(struct site source, struct site destination, char *name)
3479 {
3480         /* This function invokes gnuplot to generate an appropriate
3481            output file indicating the terrain profile between the source
3482            and destination locations.  "filename" is the name assigned
3483            to the output file generated by gnuplot.  The filename extension
3484            is used to set gnuplot's terminal setting and output file type.
3485            If no extension is found, .png is assumed.  */
3486
3487         int     x, y, z;
3488         char    filename[255], term[30], ext[15];
3489         FILE    *fd=NULL;
3490
3491         ReadPath(destination,source);
3492
3493         fd=fopen("profile.gp","wb");
3494
3495         for (x=0; x<path.length; x++)
3496         {
3497                 if (metric)
3498                         fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*path.elevation[x]);
3499
3500                 else
3501                         fprintf(fd,"%f\t%f\n",path.distance[x],path.elevation[x]);
3502         }
3503
3504         fclose(fd);
3505
3506         if (name[0]==0)
3507         {
3508                 /* Default filename and output file type */
3509
3510                 strncpy(filename,"profile\0",8);
3511                 strncpy(term,"png\0",4);
3512                 strncpy(ext,"png\0",4);
3513         }
3514
3515         else
3516         {
3517                 /* Grab extension and terminal type from "name" */
3518
3519                 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
3520                         filename[x]=name[x];
3521
3522                 if (name[x]=='.')
3523                 {
3524                         for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
3525                         {
3526                                 term[y]=tolower(name[x]);
3527                                 ext[y]=term[y];
3528                         }
3529
3530                         ext[y]=0;
3531                         term[y]=0;
3532                         filename[z]=0;
3533                 }
3534
3535                 else
3536                 {       /* No extension -- Default is png */
3537
3538                         filename[x]=0;
3539                         strncpy(term,"png\0",4);
3540                         strncpy(ext,"png\0",4);
3541                 }
3542         }
3543
3544         /* Either .ps or .postscript may be used
3545            as an extension for postscript output. */
3546
3547         if (strncmp(term,"postscript",10)==0)
3548                 strncpy(ext,"ps\0",3);
3549
3550         else if (strncmp(ext,"ps",2)==0)
3551                 strncpy(term,"postscript enhanced color\0",26);
3552
3553         fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
3554         fflush(stdout);
3555
3556         fd=fopen("splat.gp","w");
3557         fprintf(fd,"set grid\n");
3558         fprintf(fd,"set autoscale\n");
3559         fprintf(fd,"set encoding iso_8859_1\n");
3560         fprintf(fd,"set term %s\n",term);
3561         fprintf(fd,"set title \"SPLAT! Terrain Profile Between %s and %s (%.2f%c Azimuth)\"\n",destination.name, source.name, Azimuth(destination,source),176);
3562
3563         if (metric)
3564         {
3565                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(source,destination));
3566                 fprintf(fd,"set ylabel \"Ground Elevation Above Sea Level (meters)\"\n");
3567
3568
3569         }
3570
3571         else
3572         {
3573                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(source,destination));
3574                 fprintf(fd,"set ylabel \"Ground Elevation Above Sea Level (feet)\"\n");
3575         }
3576
3577         fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
3578         fprintf(fd,"plot \"profile.gp\" title \"\" with lines\n");
3579         fclose(fd);
3580
3581         x=system("gnuplot splat.gp");
3582
3583         if (x!=-1)
3584         {
3585                 unlink("splat.gp");
3586                 unlink("profile.gp");
3587                 fprintf(stdout," Done!\n");
3588                 fflush(stdout);
3589         }
3590
3591         else
3592                 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
3593 }
3594
3595 void GraphElevation(struct site source, struct site destination, char *name)
3596 {
3597         /* This function invokes gnuplot to generate an appropriate
3598            output file indicating the terrain profile between the source
3599            and destination locations.  "filename" is the name assigned
3600            to the output file generated by gnuplot.  The filename extension
3601            is used to set gnuplot's terminal setting and output file type.
3602            If no extension is found, .png is assumed. */
3603
3604         int     x, y, z;
3605         char    filename[255], term[30], ext[15];
3606         double  angle, refangle, maxangle=-90.0;
3607         struct  site remote;
3608         FILE    *fd=NULL, *fd2=NULL;
3609
3610         ReadPath(destination,source);  /* destination=RX, source=TX */
3611         refangle=ElevationAngle(destination,source);
3612
3613         fd=fopen("profile.gp","wb");
3614         fd2=fopen("reference.gp","wb");
3615
3616         for (x=1; x<path.length-1; x++)
3617         {
3618                 remote.lat=path.lat[x];
3619                 remote.lon=path.lon[x];
3620                 remote.alt=0.0;
3621                 angle=ElevationAngle(destination,remote);
3622
3623                 if (metric)
3624                 {
3625                         fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[x],angle);
3626                         fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[x],refangle);
3627                 }
3628
3629                 else
3630                 {
3631                         fprintf(fd,"%f\t%f\n",path.distance[x],angle);
3632                         fprintf(fd2,"%f\t%f\n",path.distance[x],refangle);
3633                 }
3634
3635                 if (angle>maxangle)
3636                         maxangle=angle;
3637         }
3638
3639         if (metric)
3640         {
3641                 fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],refangle);
3642                 fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],refangle);
3643         }
3644
3645         else
3646         {
3647                 fprintf(fd,"%f\t%f\n",path.distance[path.length-1],refangle);
3648                 fprintf(fd2,"%f\t%f\n",path.distance[path.length-1],refangle);
3649         }
3650
3651         fclose(fd);
3652         fclose(fd2);
3653
3654         if (name[0]==0)
3655         {
3656                 /* Default filename and output file type */
3657
3658                 strncpy(filename,"profile\0",8);
3659                 strncpy(term,"png\0",4);
3660                 strncpy(ext,"png\0",4);
3661         }
3662
3663         else
3664         {
3665                 /* Grab extension and terminal type from "name" */
3666
3667                 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
3668                         filename[x]=name[x];
3669
3670                 if (name[x]=='.')
3671                 {
3672                         for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
3673                         {
3674                                 term[y]=tolower(name[x]);
3675                                 ext[y]=term[y];
3676                         }
3677
3678                         ext[y]=0;
3679                         term[y]=0;
3680                         filename[z]=0;
3681                 }
3682
3683                 else
3684                 {       /* No extension -- Default is png */
3685
3686                         filename[x]=0;
3687                         strncpy(term,"png\0",4);
3688                         strncpy(ext,"png\0",4);
3689                 }
3690         }
3691
3692         /* Either .ps or .postscript may be used
3693            as an extension for postscript output. */
3694
3695         if (strncmp(term,"postscript",10)==0)
3696                 strncpy(ext,"ps\0",3);
3697
3698         else if (strncmp(ext,"ps",2)==0)
3699                 strncpy(term,"postscript enhanced color\0",26);
3700
3701         fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
3702         fflush(stdout);
3703
3704         fd=fopen("splat.gp","w");
3705
3706         fprintf(fd,"set grid\n");
3707         fprintf(fd,"set yrange [%2.3f to %2.3f]\n", (-fabs(refangle)-0.25), maxangle+0.25);
3708         fprintf(fd,"set encoding iso_8859_1\n");
3709         fprintf(fd,"set term %s\n",term);
3710         fprintf(fd,"set title \"SPLAT! Elevation Profile Between %s and %s (%.2f%c azimuth)\"\n",destination.name,source.name,Azimuth(destination,source),176);
3711
3712         if (metric)
3713                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(source,destination));
3714         else
3715                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(source,destination));
3716
3717
3718         fprintf(fd,"set ylabel \"Elevation Angle Along LOS Path Between %s and %s (degrees)\"\n",destination.name,source.name);
3719         fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
3720         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);
3721
3722         fclose(fd);
3723
3724         x=system("gnuplot splat.gp");
3725
3726         if (x!=-1)
3727         {
3728                 unlink("splat.gp");
3729                 unlink("profile.gp");
3730                 unlink("reference.gp");
3731
3732                 fprintf(stdout," Done!\n");
3733                 fflush(stdout);
3734         }
3735
3736         else
3737                 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
3738 }
3739
3740 void GraphHeight(struct site source, struct site destination, char *name, double f, unsigned char n)
3741 {
3742         /* This function invokes gnuplot to generate an appropriate
3743            output file indicating the terrain profile between the source
3744            and destination locations referenced to the line-of-sight path
3745            between the receive and transmit sites.  "filename" is the name
3746            assigned to the output file generated by gnuplot.  The filename
3747            extension is used to set gnuplot's terminal setting and output
3748            file type.  If no extension is found, .png is assumed. */
3749
3750         int     x, y, z;
3751         char    filename[255], term[30], ext[15];
3752         double  a, b, c, height=0.0, refangle, cangle, maxheight=-100000.0,
3753                 minheight=100000.0, lambda=0.0, f_zone=0.0, fpt6_zone=0.0,
3754                 nm=0.0, nb=0.0, ed=0.0, es=0.0, r=0.0, d=0.0, d1=0.0,
3755                 terrain, azimuth, distance, dheight=0.0, minterrain=100000.0,
3756                 minearth=100000.0, miny, maxy, min2y, max2y;
3757         struct  site remote;
3758         FILE    *fd=NULL, *fd2=NULL, *fd3=NULL, *fd4=NULL, *fd5=NULL;
3759
3760         ReadPath(destination,source);  /* destination=RX, source=TX */
3761         azimuth=Azimuth(destination,source);
3762         distance=Distance(destination,source);
3763         refangle=ElevationAngle(destination,source);
3764         b=GetElevation(destination)+destination.alt+earthradius;
3765
3766         /* Wavelength and path distance (great circle) in feet. */
3767
3768         if (f)
3769         {
3770                 lambda=9.8425e8/(f*1e6);
3771                 d=5280.0*path.distance[path.length-1];
3772         }
3773
3774         if (n)
3775         {
3776                 ed=GetElevation(destination);
3777                 es=GetElevation(source);
3778                 nb=-destination.alt-ed;
3779                 nm=(-source.alt-es-nb)/(path.distance[path.length-1]);
3780         }
3781
3782         fd=fopen("profile.gp","wb");
3783         fd2=fopen("reference.gp","wb");
3784         fd5=fopen("curvature.gp", "wb");
3785
3786         if (f)
3787         {
3788                 fd3=fopen("fresnel.gp", "wb");
3789                 fd4=fopen("fresnel_pt_6.gp", "wb");
3790         }
3791
3792         for (x=0; x<path.length; x++)
3793         {
3794                 remote.lat=path.lat[x];
3795                 remote.lon=path.lon[x];
3796                 remote.alt=0.0;
3797
3798                 terrain=GetElevation(remote);
3799
3800                 if (x==0)
3801                         terrain+=destination.alt;  /* RX antenna spike */
3802
3803                 a=terrain+earthradius;
3804                 cangle=5280.0*Distance(destination,remote)/earthradius;
3805                 c=b*sin(refangle*deg2rad+HALFPI)/sin(HALFPI-refangle*deg2rad-cangle);
3806
3807                 height=a-c;
3808
3809                 /* Per Fink and Christiansen, Electronics
3810                  * Engineers' Handbook, 1989:
3811                  *
3812                  *   H = sqrt(lamba * d1 * (d - d1)/d)
3813                  *
3814                  * where H is the distance from the LOS
3815                  * path to the first Fresnel zone boundary.
3816                  */
3817
3818                 if (f)
3819                 {
3820                         d1=5280.0*path.distance[x];
3821                         f_zone=-1*sqrt(lambda*d1*(d-d1)/d);
3822                         fpt6_zone=0.6*f_zone;
3823                 }
3824
3825                 if (n)
3826                 {
3827                         r=-(nm*path.distance[x])-nb;
3828                         height+=r;
3829
3830                         if (f>0)
3831                         {
3832                                 f_zone+=r;
3833                                 fpt6_zone+=r;
3834                         }
3835                 }
3836
3837                 else
3838                         r=0.0;
3839
3840                 if (metric)
3841                 {
3842                         fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*height);
3843                         fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*r);
3844                         fprintf(fd5,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*(height-terrain));
3845                 }
3846
3847                 else
3848                 {
3849                         fprintf(fd,"%f\t%f\n",path.distance[x],height);
3850                         fprintf(fd2,"%f\t%f\n",path.distance[x],r);
3851                         fprintf(fd5,"%f\t%f\n",path.distance[x],height-terrain);
3852                 }
3853
3854                 if (f)
3855                 {
3856                         if (metric)
3857                         {
3858                                 fprintf(fd3,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*f_zone);
3859                                 fprintf(fd4,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*fpt6_zone);
3860                         }
3861
3862                         else
3863                         {
3864                                 fprintf(fd3,"%f\t%f\n",path.distance[x],f_zone);
3865                                 fprintf(fd4,"%f\t%f\n",path.distance[x],fpt6_zone);
3866                         }
3867
3868                         if (f_zone<minheight)
3869                                 minheight=f_zone;
3870                 }
3871
3872                 if (height>maxheight)
3873                         maxheight=height;
3874
3875                 if (height<minheight)
3876                         minheight=height;
3877
3878                 if (r>maxheight)
3879                         maxheight=r;
3880
3881                 if (terrain<minterrain)
3882                         minterrain=terrain;
3883
3884                 if ((height-terrain)<minearth)
3885                         minearth=height-terrain;
3886         }
3887
3888         if (n)
3889                 r=-(nm*path.distance[path.length-1])-nb;
3890         else
3891                 r=0.0;
3892
3893         if (metric)
3894         {
3895                 fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3896                 fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3897         }
3898
3899         else
3900         {
3901                 fprintf(fd,"%f\t%f\n",path.distance[path.length-1],r);
3902                 fprintf(fd2,"%f\t%f\n",path.distance[path.length-1],r);
3903         }
3904
3905         if (f)
3906         {
3907                 if (metric)
3908                 {
3909                         fprintf(fd3,"%f\t%f\n",path.distance[path.length-1],r);
3910                         fprintf(fd4,"%f\t%f\n",path.distance[path.length-1],r);
3911                 }
3912
3913                 else
3914                 {
3915                         fprintf(fd3,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3916                         fprintf(fd4,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3917                 }
3918         }
3919
3920         if (r>maxheight)
3921                 maxheight=r;
3922
3923         if (r<minheight)
3924                 minheight=r;
3925
3926         fclose(fd);
3927         fclose(fd2);
3928         fclose(fd5);
3929
3930         if (f)
3931         {
3932                 fclose(fd3);
3933                 fclose(fd4);
3934         }
3935
3936         if (name[0]==0)
3937         {
3938                 /* Default filename and output file type */
3939
3940                 strncpy(filename,"height\0",8);
3941                 strncpy(term,"png\0",4);
3942                 strncpy(ext,"png\0",4);
3943         }
3944
3945         else
3946         {
3947                 /* Grab extension and terminal type from "name" */
3948
3949                 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
3950                         filename[x]=name[x];
3951
3952                 if (name[x]=='.')
3953                 {
3954                         for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
3955                         {
3956                                 term[y]=tolower(name[x]);
3957                                 ext[y]=term[y];
3958                         }
3959
3960                         ext[y]=0;
3961                         term[y]=0;
3962                         filename[z]=0;
3963                 }
3964
3965                 else
3966                 {       /* No extension -- Default is png */
3967
3968                         filename[x]=0;
3969                         strncpy(term,"png\0",4);
3970                         strncpy(ext,"png\0",4);
3971                 }
3972         }
3973
3974         /* Either .ps or .postscript may be used
3975            as an extension for postscript output. */
3976
3977         if (strncmp(term,"postscript",10)==0)
3978                 strncpy(ext,"ps\0",3);
3979
3980         else if (strncmp(ext,"ps",2)==0)
3981                 strncpy(term,"postscript enhanced color\0",26);
3982
3983         fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
3984         fflush(stdout);
3985
3986         fd=fopen("splat.gp","w");
3987
3988         dheight=maxheight-minheight;
3989         miny=minheight-0.15*dheight;
3990         maxy=maxheight+0.05*dheight;
3991
3992         if (maxy<20.0)
3993                 maxy=20.0;
3994
3995         dheight=maxheight-minheight;
3996         min2y=miny-minterrain+0.05*dheight;
3997
3998         if (minearth<min2y)
3999         {
4000                 miny-=min2y-minearth+0.05*dheight;
4001                 min2y=minearth-0.05*dheight;
4002         }
4003
4004         max2y=min2y+maxy-miny;
4005
4006         fprintf(fd,"set grid\n");
4007         fprintf(fd,"set yrange [%2.3f to %2.3f]\n", metric?miny*METERS_PER_FOOT:miny, metric?maxy*METERS_PER_FOOT:maxy);
4008         fprintf(fd,"set y2range [%2.3f to %2.3f]\n", metric?min2y*METERS_PER_FOOT:min2y, metric?max2y*METERS_PER_FOOT:max2y);
4009         fprintf(fd,"set xrange [-0.5 to %2.3f]\n",metric?KM_PER_MILE*rint(distance+0.5):rint(distance+0.5));
4010         fprintf(fd,"set encoding iso_8859_1\n");
4011         fprintf(fd,"set term %s\n",term);
4012
4013         if (f)
4014                 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);
4015
4016         else
4017                 fprintf(fd,"set title \"SPLAT! Height Profile Between %s and %s (%.2f%c azimuth)\"\n",destination.name, source.name, azimuth,176);
4018
4019         if (metric)
4020                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(source,destination));
4021         else
4022                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(source,destination));
4023
4024         if (n)
4025         {
4026                 if (metric)
4027                         fprintf(fd,"set ylabel \"Normalized Height Referenced To LOS Path Between\\n%s and %s (meters)\"\n",destination.name,source.name);
4028
4029                 else
4030                         fprintf(fd,"set ylabel \"Normalized Height Referenced To LOS Path Between\\n%s and %s (feet)\"\n",destination.name,source.name);
4031
4032         }
4033
4034         else
4035         {
4036                 if (metric)
4037                         fprintf(fd,"set ylabel \"Height Referenced To LOS Path Between %s and %s (meters)\"\n",destination.name,source.name);
4038
4039                 else
4040                         fprintf(fd,"set ylabel \"Height Referenced To LOS Path Between %s and %s (feet)\"\n",destination.name,source.name);
4041         }
4042
4043         fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
4044
4045         if (f)
4046                 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);
4047
4048         else
4049                 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");
4050
4051         fclose(fd);
4052
4053         x=system("gnuplot splat.gp");
4054
4055         if (x!=-1)
4056         {
4057                 unlink("splat.gp");
4058                 unlink("profile.gp");
4059                 unlink("reference.gp");
4060                 unlink("curvature.gp");
4061
4062                 if (f)
4063                 {
4064                         unlink("fresnel.gp");
4065                         unlink("fresnel_pt_6.gp");
4066                 }
4067
4068                 fprintf(stdout," Done!\n");
4069                 fflush(stdout);
4070         }
4071
4072         else
4073                 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
4074 }
4075
4076 void GraphLongley(struct site source, struct site destination, char *name)
4077 {
4078         /* This function invokes gnuplot to generate an appropriate
4079            output file indicating the Longley-Rice model loss between
4080            the source and destination locations.  "filename" is the
4081            name assigned to the output file generated by gnuplot.
4082            The filename extension is used to set gnuplot's terminal
4083            setting and output file type.  If no extension is found,
4084            .png is assumed. */
4085
4086         int     x, y, z, errnum, errflag=0;
4087         char    filename[255], term[30], ext[15], strmode[100],
4088                 report_name[80], block=0;
4089         double  maxloss=-100000.0, minloss=100000.0, loss, haavt,
4090                 angle1, angle2, azimuth, pattern=1.0, patterndB=0.0,
4091                 total_loss=0.0, cos_xmtr_angle, cos_test_angle=0.0,
4092                 source_alt, test_alt, dest_alt, source_alt2, dest_alt2,
4093                 distance, elevation, four_thirds_earth;
4094         FILE    *fd=NULL, *fd2=NULL;
4095
4096         sprintf(report_name,"%s-to-%s.lro",source.name,destination.name);
4097
4098         four_thirds_earth=EARTHRADIUS*(4.0/3.0);
4099
4100         for (x=0; report_name[x]!=0; x++)
4101                 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
4102                         report_name[x]='_';
4103
4104         fd2=fopen(report_name,"w");
4105
4106         fprintf(fd2,"\n\t--==[ SPLAT! v%s Longley-Rice Model Path Loss Report ]==--\n\n",splat_version);
4107         fprintf(fd2,"Analysis of RF path conditions between %s and %s:\n",source.name, destination.name);
4108         fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4109         fprintf(fd2,"Transmitter site: %s\n",source.name);
4110
4111         if (source.lat>=0.0)
4112         {
4113                 fprintf(fd2,"Site location: %.4f North / %.4f West",source.lat, source.lon);
4114                 fprintf(fd2, " (%s N / ", dec2dms(source.lat));
4115         }
4116
4117         else
4118         {
4119
4120                 fprintf(fd2,"Site location: %.4f South / %.4f West",-source.lat, source.lon);
4121                 fprintf(fd2, " (%s S / ", dec2dms(source.lat));
4122         }
4123
4124         fprintf(fd2, "%s W)\n", dec2dms(source.lon));
4125
4126         if (metric)
4127         {
4128                 fprintf(fd2,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(source));
4129                 fprintf(fd2,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*source.alt,METERS_PER_FOOT*(source.alt+GetElevation(source)));
4130         }
4131
4132         else
4133         {
4134                 fprintf(fd2,"Ground elevation: %.2f feet AMSL\n",GetElevation(source));
4135                 fprintf(fd2,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",source.alt, source.alt+GetElevation(source));
4136         }
4137
4138         haavt=haat(source);
4139
4140         if (haavt>-4999.0)
4141         {
4142                 if (metric)
4143                         fprintf(fd2,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4144                 else
4145                         fprintf(fd2,"Antenna height above average terrain: %.2f feet\n",haavt);
4146         }
4147
4148         azimuth=Azimuth(source,destination);
4149         angle1=ElevationAngle(source,destination);
4150         angle2=ElevationAngle2(source,destination,earthradius);
4151
4152         if (got_azimuth_pattern || got_elevation_pattern)
4153         {
4154                 x=(int)rint(10.0*(10.0-angle2));
4155
4156                 if (x>=0 && x<=1000)
4157                         pattern=(double)LR.antenna_pattern[(int)rint(azimuth)][x];
4158
4159                 patterndB=20.0*log10(pattern);
4160
4161                 fprintf(fd2,"Antenna pattern between %s and %s: %.3f (%.2f dB)\n",source.name, destination.name, pattern, patterndB);
4162
4163         }
4164
4165         if (metric)
4166                 fprintf(fd2,"Distance to %s: %.2f kilometers\n",destination.name,METERS_PER_FOOT*Distance(source,destination));
4167
4168         else
4169                 fprintf(fd2,"Distance to %s: %.2f miles.\n",destination.name,Distance(source,destination));
4170
4171         fprintf(fd2,"Azimuth to %s: %.2f degrees\n",destination.name,azimuth);
4172
4173         if (angle1>=0.0)
4174                 fprintf(fd2,"Elevation angle to %s: %+.4f degrees\n",destination.name,angle1);
4175
4176         else
4177                 fprintf(fd2,"Depression angle to %s: %+.4f degrees\n",destination.name,angle1);
4178
4179         if (angle1!=angle2)
4180         {
4181                 if (angle2<0.0)
4182                         fprintf(fd2,"Depression");
4183                 else
4184                         fprintf(fd2,"Elevation");
4185
4186                 fprintf(fd2," angle to the first obstruction: %+.4f degrees\n",angle2);
4187         }
4188
4189         fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4190
4191         /* Receiver */
4192
4193         fprintf(fd2,"Receiver site: %s\n",destination.name);
4194
4195         if (destination.lat>=0.0)
4196         {
4197                 fprintf(fd2,"Site location: %.4f North / %.4f West",destination.lat, destination.lon);
4198                 fprintf(fd2, " (%s N / ", dec2dms(destination.lat));
4199         }
4200
4201         else
4202         {
4203                 fprintf(fd2,"Site location: %.4f South / %.4f West",-destination.lat, destination.lon);
4204                 fprintf(fd2, " (%s S / ", dec2dms(destination.lat));
4205         }
4206
4207         fprintf(fd2, "%s W)\n", dec2dms(destination.lon));
4208
4209         if (metric)
4210         {
4211                 fprintf(fd2,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(destination));
4212                 fprintf(fd2,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*destination.alt, METERS_PER_FOOT*(destination.alt+GetElevation(destination)));
4213         }
4214
4215         else
4216         {
4217                 fprintf(fd2,"Ground elevation: %.2f feet AMSL\n",GetElevation(destination));
4218                 fprintf(fd2,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",destination.alt, destination.alt+GetElevation(destination));
4219         }
4220
4221         haavt=haat(destination);
4222
4223         if (haavt>-4999.0)
4224         {
4225                 if (metric)
4226                         fprintf(fd2,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4227                 else
4228                         fprintf(fd2,"Antenna height above average terrain: %.2f feet\n",haavt);
4229         }
4230
4231         if (metric)
4232                 fprintf(fd2,"Distance to %s: %.2f kilometers\n",source.name,KM_PER_MILE*Distance(source,destination));
4233
4234         else
4235                 fprintf(fd2,"Distance to %s: %.2f miles\n",source.name,Distance(source,destination));
4236
4237         azimuth=Azimuth(destination,source);
4238
4239         angle1=ElevationAngle(destination,source);
4240         angle2=ElevationAngle2(destination,source,earthradius);
4241
4242         fprintf(fd2,"Azimuth to %s: %.2f degrees.\n",source.name,azimuth);
4243
4244         if (angle1>=0.0)
4245                 fprintf(fd2,"Elevation angle to %s: %+.4f degrees\n",source.name,angle1);
4246
4247         else
4248                 fprintf(fd2,"Depression angle to %s: %+.4f degrees\n",source.name,angle1);
4249
4250         if (angle1!=angle2)
4251         {
4252                 if (angle2<0.0)
4253                         fprintf(fd2,"Depression");
4254                 else
4255                         fprintf(fd2,"Elevation");
4256
4257                 fprintf(fd2," angle to the first obstruction: %+.4f degrees\n",angle2);
4258         }
4259
4260         fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4261
4262         fprintf(fd2,"Longley-Rice path calculation parameters used in this analysis:\n\n");
4263         fprintf(fd2,"Earth's Dielectric Constant: %.3lf\n",LR.eps_dielect);
4264         fprintf(fd2,"Earth's Conductivity: %.3lf\n",LR.sgm_conductivity);
4265         fprintf(fd2,"Atmospheric Bending Constant (N): %.3lf\n",LR.eno_ns_surfref);
4266         fprintf(fd2,"Frequency: %.3lf (MHz)\n",LR.frq_mhz);
4267         fprintf(fd2,"Radio Climate: %d (",LR.radio_climate);
4268
4269         switch (LR.radio_climate)
4270         {
4271                 case 1:
4272                 fprintf(fd2,"Equatorial");
4273                 break;
4274
4275                 case 2:
4276                 fprintf(fd2,"Continental Subtropical");
4277                 break;
4278
4279                 case 3:
4280                 fprintf(fd2,"Maritime Subtropical");
4281                 break;
4282
4283                 case 4:
4284                 fprintf(fd2,"Desert");
4285                 break;
4286
4287                 case 5:
4288                 fprintf(fd2,"Continental Temperate");
4289                 break;
4290
4291                 case 6:
4292                 fprintf(fd2,"Martitime Temperate, Over Land");
4293                 break;
4294
4295                 case 7:
4296                 fprintf(fd2,"Maritime Temperate, Over Sea");
4297                 break;
4298
4299                 default:
4300                 fprintf(fd2,"Unknown");
4301         }
4302
4303         fprintf(fd2,")\nPolarization: %d (",LR.pol);
4304
4305         if (LR.pol==0)
4306                 fprintf(fd2,"Horizontal");
4307
4308         if (LR.pol==1)
4309                 fprintf(fd2,"Vertical");
4310
4311         fprintf(fd2,")\nFraction of Situations: %.1lf%c\n",LR.conf*100.0,37);
4312         fprintf(fd2,"Fraction of Time: %.1lf%c\n",LR.rel*100.0,37);
4313         fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4314
4315         fprintf(fd2,"Analysis Results:\n\n");
4316
4317         ReadPath(source, destination);  /* source=TX, destination=RX */
4318
4319         /* Copy elevations along path into the elev_l[] array. */
4320
4321         for (x=0; x<path.length; x++)
4322                 elev_l[x+2]=path.elevation[x]*METERS_PER_FOOT;
4323
4324         if (metric)
4325                 fprintf(fd2,"Distance (km)");
4326         else
4327                 fprintf(fd2,"Distance (mi)");
4328
4329         fprintf(fd2,"\tLoss (dB)\tErrnum\tComment\n\n");
4330
4331         fd=fopen("profile.gp","w");
4332
4333         azimuth=rint(Azimuth(source,destination));
4334
4335         for (y=2; y<(path.length-1); y++)  /* path.length-1 avoids LR error */
4336         {
4337                 distance=5280.0*path.distance[y];
4338                 source_alt=four_thirds_earth+source.alt+path.elevation[0];
4339                 dest_alt=four_thirds_earth+destination.alt+path.elevation[y];
4340                 dest_alt2=dest_alt*dest_alt;
4341                 source_alt2=source_alt*source_alt;
4342
4343                 /* Calculate the cosine of the elevation of
4344                    the receiver as seen by the transmitter. */
4345
4346                 cos_xmtr_angle=((source_alt2)+(distance*distance)-(dest_alt2))/(2.0*source_alt*distance);
4347
4348                 if (got_elevation_pattern)
4349                 {
4350                         /* If an antenna elevation pattern is available, the
4351                            following code determines the elevation angle to
4352                            the first obstruction along the path. */
4353
4354                         for (x=2, block=0; x<y && block==0; x++)
4355                         {
4356                                 distance=5280.0*(path.distance[y]-path.distance[x]);
4357                                 test_alt=four_thirds_earth+path.elevation[x];
4358
4359                                 /* Calculate the cosine of the elevation
4360                                    angle of the terrain (test point)
4361                                    as seen by the transmitter. */
4362
4363                                 cos_test_angle=((source_alt2)+(distance*distance)-(test_alt*test_alt))/(2.0*source_alt*distance);
4364
4365                                 /* Compare these two angles to determine if
4366                                    an obstruction exists.  Since we're comparing
4367                                    the cosines of these angles rather than
4368                                    the angles themselves, the sense of the
4369                                    following "if" statement is reversed from
4370                                    what it would be if the angles themselves
4371                                    were compared. */
4372
4373                                 if (cos_xmtr_angle>cos_test_angle)
4374                                         block=1;
4375                         }
4376
4377                         /* At this point, we have the elevation angle
4378                            to the first obstruction (if it exists). */
4379                 }
4380
4381                 /* Determine path loss for each point along the
4382                    path using Longley-Rice's point_to_point mode
4383                    starting at x=2 (number_of_points = 1), the
4384                    shortest distance terrain can play a role in
4385                    path loss. */
4386
4387                 elev_l[0]=y-1;  /* (number of points - 1) */
4388
4389                 /* Distance between elevation samples */
4390                 elev_l[1]=METERS_PER_MILE*(path.distance[y]-path.distance[y-1]);
4391
4392                 point_to_point(elev_l, source.alt*METERS_PER_FOOT,
4393                 destination.alt*METERS_PER_FOOT, LR.eps_dielect,
4394                 LR.sgm_conductivity, LR.eno_ns_surfref, LR.frq_mhz,
4395                 LR.radio_climate, LR.pol, LR.conf, LR.rel, loss,
4396                 strmode, errnum);
4397
4398                 if (block)
4399                         elevation=((acos(cos_test_angle))/deg2rad)-90.0;
4400                 else
4401                         elevation=((acos(cos_xmtr_angle))/deg2rad)-90.0;
4402
4403                 /* Integrate the antenna's radiation
4404                    pattern into the overall path loss. */
4405
4406                 x=(int)rint(10.0*(10.0-elevation));
4407
4408                 if (x>=0 && x<=1000)
4409                 {
4410                         pattern=(double)LR.antenna_pattern[(int)azimuth][x];
4411
4412                         if (pattern!=0.0)
4413                                 patterndB=20.0*log10(pattern);
4414                 }
4415
4416                 else
4417                         patterndB=0.0;
4418
4419                 total_loss=loss-patterndB;
4420
4421                 if (metric)
4422                 {
4423                         fprintf(fd,"%f\t%f\n",KM_PER_MILE*(path.distance[path.length-1]-path.distance[y]),total_loss);
4424                         fprintf(fd2,"%7.2f\t\t%7.2f\t\t  %d\t%s\n",KM_PER_MILE*path.distance[y],total_loss, errnum, strmode);
4425                 }
4426
4427                 else
4428                 {
4429                         fprintf(fd,"%f\t%f\n",path.distance[path.length-1]-path.distance[y],total_loss);
4430                         fprintf(fd2,"%7.2f\t\t%7.2f\t\t  %d\t%s\n",path.distance[y],total_loss, errnum, strmode);
4431                 }
4432
4433                 errflag|=errnum;
4434
4435                 if (total_loss>maxloss)
4436                         maxloss=total_loss;
4437
4438                 if (total_loss<minloss)
4439                         minloss=total_loss;
4440         }
4441
4442         fclose(fd);
4443
4444         fprintf(fd2,"\nLongley-Rice path loss between %s and %s is %.2f dB.\n",source.name, destination.name, loss);
4445
4446         if (patterndB!=0.0)
4447                 fprintf(fd2,"Total path loss including TX antenna pattern is %.2f dB.\n",total_loss);
4448
4449         if (errflag)
4450         {
4451                 fprintf(fd2,"\nNotes on \"errnum\":\n\n");
4452                 fprintf(fd2,"  0: No error.  :-)\n");
4453                 fprintf(fd2,"  1: Warning!  Some parameters are nearly out of range.\n");
4454                 fprintf(fd2,"     Results should be used with caution.\n");
4455                 fprintf(fd2,"  2: Note: Default parameters have been substituted for impossible ones.\n");
4456                 fprintf(fd2,"  3: Warning!  A combination of parameters is out of range.\n");
4457                 fprintf(fd2,"     Results are probably invalid.\n");
4458                 fprintf(fd2,"  Other: Warning!  Some parameters are out of range.\n");
4459                 fprintf(fd2,"    Results are probably invalid.\n\nEnd of Report\n");
4460         }
4461
4462         fprintf(stdout,"Longley-Rice path loss between %s and %s is %.2f dB.\n",source.name, destination.name, loss);
4463
4464         if (patterndB!=0.0)
4465                 fprintf(stdout,"Total path loss including TX antenna pattern is %.2f dB.\n",total_loss);
4466
4467         fprintf(stdout,"Path Loss Report written to: \"%s\"\n",report_name);
4468         fflush(stdout);
4469
4470         fclose(fd2);
4471
4472         if (name[0]==0)
4473         {
4474                 /* Default filename and output file type */
4475
4476                 strncpy(filename,"loss\0",5);
4477                 strncpy(term,"png\0",4);
4478                 strncpy(ext,"png\0",4);
4479         }
4480
4481         else
4482         {
4483                 /* Grab extension and terminal type from "name" */
4484
4485                 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
4486                         filename[x]=name[x];
4487
4488                 if (name[x]=='.')
4489                 {
4490                         for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
4491                         {
4492                                 term[y]=tolower(name[x]);
4493                                 ext[y]=term[y];
4494                         }
4495
4496                         ext[y]=0;
4497                         term[y]=0;
4498                         filename[z]=0;
4499                 }
4500
4501                 else
4502                 {       /* No extension -- Default is png */
4503
4504                         filename[x]=0;
4505                         strncpy(term,"png\0",4);
4506                         strncpy(ext,"png\0",4);
4507                 }
4508         }
4509
4510         /* Either .ps or .postscript may be used
4511            as an extension for postscript output. */
4512
4513         if (strncmp(term,"postscript",10)==0)
4514                 strncpy(ext,"ps\0",3);
4515
4516         else if (strncmp(ext,"ps",2)==0)
4517                 strncpy(term,"postscript enhanced color\0",26);
4518
4519         fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
4520         fflush(stdout);
4521
4522         fd=fopen("splat.gp","w");
4523
4524         fprintf(fd,"set grid\n");
4525         fprintf(fd,"set yrange [%2.3f to %2.3f]\n", minloss, maxloss);
4526         fprintf(fd,"set encoding iso_8859_1\n");
4527         fprintf(fd,"set term %s\n",term);
4528         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);
4529
4530         if (metric)
4531                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(destination,source));
4532         else
4533                 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(destination,source));
4534
4535         if (got_azimuth_pattern || got_elevation_pattern)
4536                 fprintf(fd,"set ylabel \"Total Path Loss (including TX antenna pattern) (dB)");
4537         else
4538                 fprintf(fd,"set ylabel \"Longley-Rice Path Loss (dB)");
4539
4540         fprintf(fd,"\"\nset output \"%s.%s\"\n",filename,ext);
4541         fprintf(fd,"plot \"profile.gp\" title \"Path Loss\" with lines\n");
4542
4543         fclose(fd);
4544
4545         x=system("gnuplot splat.gp");
4546
4547         if (x!=-1)
4548         {
4549                 unlink("splat.gp");
4550                 unlink("profile.gp");
4551                 unlink("reference.gp");
4552
4553                 fprintf(stdout," Done!\n");
4554                 fflush(stdout);
4555         }
4556
4557         else
4558                 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
4559 }
4560
4561 void ObstructionReport(struct site xmtr, struct site rcvr, char report, double f)
4562 {
4563         int     x;
4564         struct  site site_x;
4565         double  h_r, h_t, h_x, h_r_orig, cos_tx_angle, cos_test_angle,
4566                 cos_tx_angle_f1, cos_tx_angle_fpt6, haavt, d_tx, d_x,
4567                 h_r_f1, h_r_fpt6, h_f, h_los, lambda=0.0, azimuth,
4568                 pattern, patterndB, distance, angle1, angle2;
4569         char    report_name[80], string[255], string_fpt6[255],
4570                 string_f1[255];
4571         FILE    *fd;
4572
4573         sprintf(report_name,"%s-to-%s.txt",xmtr.name,rcvr.name);
4574
4575         for (x=0; report_name[x]!=0; x++)
4576                 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
4577                         report_name[x]='_';
4578
4579         fd=fopen(report_name,"w");
4580
4581         fprintf(fd,"\n\t\t--==[ SPLAT! v%s Obstruction Report ]==--\n\n",splat_version);
4582         fprintf(fd,"Analysis of great circle path between %s and %s:\n",xmtr.name, rcvr.name);
4583         fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
4584         fprintf(fd,"Transmitter site: %s\n",xmtr.name);
4585
4586
4587         if (xmtr.lat>=0.0)
4588         {
4589                 fprintf(fd,"Site location: %.4f North / %.4f West",xmtr.lat, xmtr.lon);
4590                 fprintf(fd, " (%s N / ", dec2dms(xmtr.lat));
4591         }
4592
4593         else
4594         {
4595                 fprintf(fd,"Site location: %.4f South / %.4f West",-xmtr.lat, xmtr.lon);
4596                 fprintf(fd, " (%s S / ", dec2dms(xmtr.lat));
4597         }
4598
4599         fprintf(fd, "%s W)\n", dec2dms(xmtr.lon));
4600
4601         if (metric)
4602         {
4603                 fprintf(fd,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(xmtr));
4604                 fprintf(fd,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*xmtr.alt, METERS_PER_FOOT*(xmtr.alt+GetElevation(xmtr)));
4605         }
4606
4607         else
4608         {
4609                 fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(xmtr));
4610                 fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",xmtr.alt, xmtr.alt+GetElevation(xmtr));
4611         }
4612
4613         haavt=haat(xmtr);
4614
4615         if (haavt>-4999.0)
4616         {
4617                 if (metric)
4618                         fprintf(fd,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4619                 else
4620                         fprintf(fd,"Antenna height above average terrain: %.2f feet\n",haavt);
4621         }
4622
4623         pattern=1.0;
4624         patterndB=0.0;
4625         distance=Distance(xmtr,rcvr);
4626         azimuth=Azimuth(xmtr,rcvr);
4627         angle1=ElevationAngle(xmtr,rcvr);
4628         angle2=ElevationAngle2(xmtr,rcvr,earthradius);
4629
4630         if (got_azimuth_pattern || got_elevation_pattern)
4631         {
4632                 x=(int)rint(10.0*(10.0-angle2));
4633
4634                 if (x>=0 && x<=1000)
4635                         pattern=(double)LR.antenna_pattern[(int)rint(azimuth)][x];
4636
4637                 if (pattern!=1.0)
4638                 {
4639                         fprintf(fd,"Antenna pattern toward %s: %.3f",rcvr.name,pattern);
4640                         patterndB=20.0*log10(pattern);
4641                         fprintf(fd," (%.2f dB)\n",patterndB);
4642                 }
4643         }
4644
4645         if (metric)
4646                 fprintf(fd,"Distance to %s: %.2f kilometers\n",rcvr.name,KM_PER_MILE*distance);
4647
4648         else
4649                 fprintf(fd,"Distance to %s: %.2f miles\n",rcvr.name,distance);
4650
4651         fprintf(fd,"Azimuth to %s: %.2f degrees\n",rcvr.name,azimuth);
4652
4653         if (angle1>=0.0)
4654                 fprintf(fd,"Elevation angle to %s: %+.4f degrees\n",rcvr.name,angle1);
4655
4656         else
4657                 fprintf(fd,"Depression angle to %s: %+.4f degrees\n",rcvr.name,angle1);
4658
4659         if (angle1!=angle2)
4660         {
4661                 if (angle2<0.0)
4662                         fprintf(fd,"Depression");
4663                 else
4664                         fprintf(fd,"Elevation");
4665
4666                 fprintf(fd," angle to the first obstruction: %+.4f degrees\n",angle2);
4667         }
4668
4669         fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
4670
4671         /* Receiver */
4672
4673         fprintf(fd,"Receiver site: %s\n",rcvr.name);
4674
4675         if (rcvr.lat>=0.0)
4676         {
4677                 fprintf(fd,"Site location: %.4f North / %.4f West",rcvr.lat, rcvr.lon);
4678                 fprintf(fd, " (%s N / ", dec2dms(rcvr.lat));
4679         }
4680
4681         else
4682         {
4683                 fprintf(fd,"Site location: %.4f South / %.4f West",-rcvr.lat, rcvr.lon);
4684                 fprintf(fd, " (%s S / ", dec2dms(rcvr.lat));
4685         }
4686
4687         fprintf(fd, "%s W)\n", dec2dms(rcvr.lon));
4688
4689         if (metric)
4690         {
4691                 fprintf(fd,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(rcvr));
4692                 fprintf(fd,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*rcvr.alt, METERS_PER_FOOT*(rcvr.alt+GetElevation(rcvr)));
4693         }
4694
4695         else
4696         {
4697                 fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(rcvr));
4698                 fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",rcvr.alt, rcvr.alt+GetElevation(rcvr));
4699         }
4700
4701         haavt=haat(rcvr);
4702
4703         if (haavt>-4999.0)
4704         {
4705                 if (metric)
4706                         fprintf(fd,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4707                 else
4708                         fprintf(fd,"Antenna height above average terrain: %.2f feet\n",haavt);
4709         }
4710
4711         azimuth=Azimuth(rcvr,xmtr);
4712         angle1=ElevationAngle(rcvr,xmtr);
4713         angle2=ElevationAngle2(rcvr,xmtr,earthradius);
4714
4715         if (metric)
4716                 fprintf(fd,"Distance to %s: %.2f kilometers\n",xmtr.name,KM_PER_MILE*distance);
4717         else
4718                 fprintf(fd,"Distance to %s: %.2f miles\n",xmtr.name,distance);
4719
4720         fprintf(fd,"Azimuth to %s: %.2f degrees\n",xmtr.name,azimuth);
4721
4722         if (angle1>=0.0)
4723                 fprintf(fd,"Elevation to %s: %+.4f degrees\n",xmtr.name,angle1);
4724
4725         else
4726                 fprintf(fd,"Depression angle to %s: %+.4f degrees\n",xmtr.name,angle1);
4727
4728         if (angle1!=angle2)
4729         {
4730                 if (angle2<0.0)
4731                         fprintf(fd,"Depression");
4732                 else
4733                         fprintf(fd,"Elevation");
4734
4735                 fprintf(fd," angle to the first obstruction: %+.4f degrees\n",angle2);
4736
4737         }
4738
4739         fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
4740
4741         if (report=='y')
4742         {
4743                 /* Generate profile of the terrain.  Create the path
4744                    from transmitter to receiver because that's the
4745                    way the original los() function did it, and going
4746                    the other way can yield slightly different results. */
4747
4748                 ReadPath(xmtr,rcvr);
4749                 h_r=GetElevation(rcvr)+rcvr.alt+earthradius;
4750                 h_r_f1=h_r;
4751                 h_r_fpt6=h_r;
4752                 h_r_orig=h_r;
4753                 h_t=GetElevation(xmtr)+xmtr.alt+earthradius;
4754                 d_tx=5280.0*Distance(rcvr,xmtr);
4755                 cos_tx_angle=((h_r*h_r)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r*d_tx);
4756                 cos_tx_angle_f1=cos_tx_angle;
4757                 cos_tx_angle_fpt6=cos_tx_angle;
4758
4759                 if (f)
4760                         lambda=9.8425e8/(f*1e6);
4761
4762                 /* At each point along the path calculate the cosine
4763                    of a sort of "inverse elevation angle" at the receiver.
4764                    From the antenna, 0 deg. looks at the ground, and 90 deg.
4765                    is parallel to the ground.
4766
4767                    Start at the receiver.  If this is the lowest antenna,
4768                    then terrain obstructions will be nearest to it.  (Plus,
4769                    that's the way the original los() did it.)
4770
4771                    Calculate cosines only.  That's sufficient to compare
4772                    angles and it saves the extra computational burden of
4773                    acos().  However, note the inverted comparison: if
4774                    acos(A) > acos(B), then B > A. */
4775
4776                 for (x=path.length-1; x>0; x--)
4777                 {
4778                         site_x.lat=path.lat[x];
4779                         site_x.lon=path.lon[x];
4780                         site_x.alt=0.0;
4781
4782                         h_x=GetElevation(site_x)+earthradius;
4783                         d_x=5280.0*Distance(rcvr,site_x);
4784
4785                         /* Deal with the LOS path first. */
4786
4787                         cos_test_angle=((h_r*h_r)+(d_x*d_x)-(h_x*h_x))/(2.0*h_r*d_x);
4788
4789                         if (cos_tx_angle>cos_test_angle)
4790                         {
4791                                 if (h_r==h_r_orig)
4792                                         fprintf(fd,"SPLAT! detected obstructions at:\n\n");
4793
4794                                 if (site_x.lat>=0.0)
4795                                 {
4796                                         if (metric)
4797                                                 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));
4798                                         else
4799                                                 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);
4800                                 }
4801
4802                                 else
4803                                 {
4804                                         if (metric)
4805                                                 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));
4806                                         else
4807
4808                                                 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);
4809                                 }
4810                         }
4811
4812                         while (cos_tx_angle>cos_test_angle)
4813                         {
4814                                 h_r+=1;
4815                                 cos_test_angle=((h_r*h_r)+(d_x*d_x)-(h_x*h_x))/(2.0*h_r*d_x);
4816                                 cos_tx_angle=((h_r*h_r)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r*d_tx);
4817                         }
4818
4819                         if (f)
4820                         {
4821                                 /* Now clear the first Fresnel zone, but don't
4822                                    clutter the obstruction report. */
4823
4824                                 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);
4825                                 h_los=sqrt(h_r_f1*h_r_f1+d_x*d_x-2*h_r_f1*d_x*cos_tx_angle_f1);
4826                                 h_f=h_los-sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4827
4828                                 while (h_f<h_x)
4829                                 {
4830                                         h_r_f1+=1;
4831                                         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);
4832                                         h_los=sqrt(h_r_f1*h_r_f1+d_x*d_x-2*h_r_f1*d_x*cos_tx_angle_f1);
4833                                         h_f=h_los-sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4834                                 }
4835
4836                                 /* And clear the 60% F1 zone. */
4837
4838                                 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);
4839                                 h_los=sqrt(h_r_fpt6*h_r_fpt6+d_x*d_x-2*h_r_fpt6*d_x*cos_tx_angle_fpt6);
4840                                 h_f=h_los-0.6*sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4841
4842                                 while (h_f<h_x)
4843                                 {
4844                                         h_r_fpt6+=1;
4845                                         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);
4846                                         h_los=sqrt(h_r_fpt6*h_r_fpt6+d_x*d_x-2*h_r_fpt6*d_x*cos_tx_angle_fpt6);
4847                                         h_f=h_los-0.6*sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4848                                 }
4849                         }
4850                 }
4851
4852                 if (h_r>h_r_orig)
4853                 {
4854                         if (metric)
4855                                 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));
4856                         else
4857                                 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);
4858                 }
4859
4860                 else
4861                         sprintf(string,"\nNo obstructions to LOS path due to terrain were detected by SPLAT!\n");
4862
4863                 if (f)
4864                 {
4865                         if (h_r_fpt6>h_r_orig)
4866                         {
4867                                 if (metric)
4868                                         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);
4869
4870                                 else
4871                                         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);
4872                         }
4873
4874                         else
4875                                 sprintf(string_fpt6,"\n60%c of the first Fresnel zone is clear.\n",37);
4876
4877                         if (h_r_f1>h_r_orig)
4878                         {
4879                                 if (metric)
4880                                         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));
4881
4882                                 else
4883                                         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);
4884
4885                         }
4886
4887                         else
4888                             sprintf(string_f1,"\nThe first Fresnel zone is clear.\n\n");
4889                 }
4890         }
4891
4892         fprintf(fd,"%s",string);
4893
4894         if (f)
4895         {
4896                 fprintf(fd,"%s",string_f1);
4897                 fprintf(fd,"%s",string_fpt6);
4898         }
4899
4900         fclose(fd);
4901
4902         /* Display report summary on terminal */
4903
4904         /* Line-of-sight status */
4905
4906         fprintf(stdout,"%s",string);
4907
4908         if (f)
4909         {
4910                 /* Fresnel zone status */
4911
4912                 fprintf(stdout,"%s",string_f1);
4913                 fprintf(stdout,"%s",string_fpt6);
4914         }
4915
4916         fprintf(stdout, "\nObstruction report written to: \"%s\"\n",report_name);
4917
4918         fflush(stdout);
4919 }
4920
4921 void SiteReport(struct site xmtr)
4922 {
4923         char    report_name[80];
4924         double  terrain;
4925         int     x, azi;
4926         FILE    *fd;
4927
4928         sprintf(report_name,"%s-site_report.txt",xmtr.name);
4929
4930         for (x=0; report_name[x]!=0; x++)
4931                 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
4932                         report_name[x]='_';
4933
4934         fd=fopen(report_name,"w");
4935
4936         fprintf(fd,"\n\t--==[ SPLAT! v%s Site Analysis Report For: %s ]==--\n\n",splat_version,xmtr.name);
4937
4938         fprintf(fd,"---------------------------------------------------------------------------\n\n");
4939
4940         if (xmtr.lat>=0.0)
4941         {
4942                 fprintf(fd,"Site location: %.4f North / %.4f West",xmtr.lat, xmtr.lon);
4943                 fprintf(fd, " (%s N / ",dec2dms(xmtr.lat));
4944         }
4945
4946         else
4947         {
4948                 fprintf(fd,"Site location: %.4f South / %.4f West",-xmtr.lat, xmtr.lon);
4949                 fprintf(fd, " (%s S / ",dec2dms(xmtr.lat));
4950         }
4951
4952         fprintf(fd, "%s W)\n",dec2dms(xmtr.lon));
4953
4954         if (metric)
4955         {
4956                 fprintf(fd,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(xmtr));
4957                 fprintf(fd,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*xmtr.alt, METERS_PER_FOOT*(xmtr.alt+GetElevation(xmtr)));
4958         }
4959
4960         else
4961         {
4962                 fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(xmtr));
4963                 fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",xmtr.alt, xmtr.alt+GetElevation(xmtr));
4964         }
4965
4966         terrain=haat(xmtr);
4967
4968         if (terrain>-4999.0)
4969         {
4970                 if (metric)
4971                         fprintf(fd,"Antenna height above average terrain: %.2f meters\n\n",METERS_PER_FOOT*terrain);
4972                 else
4973                         fprintf(fd,"Antenna height above average terrain: %.2f feet\n\n",terrain);
4974
4975                 /* Display the average terrain between 2 and 10 miles
4976                    from the transmitter site at azimuths of 0, 45, 90,
4977                    135, 180, 225, 270, and 315 degrees. */
4978
4979                 for (azi=0; azi<=315; azi+=45)
4980                 {
4981                         fprintf(fd,"Average terrain at %3d degrees azimuth: ",azi);
4982                         terrain=AverageTerrain(xmtr,(double)azi,2.0,10.0);
4983
4984                         if (terrain>-4999.0)
4985                         {
4986                                 if (metric)
4987                                         fprintf(fd,"%.2f meters AMSL\n",METERS_PER_FOOT*terrain);
4988                                 else
4989                                         fprintf(fd,"%.2f feet AMSL\n",terrain);
4990                         }
4991
4992                         else
4993                                 fprintf(fd,"No terrain\n");
4994                 }
4995         }
4996
4997         fprintf(fd,"\n---------------------------------------------------------------------------\n\n");
4998         fclose(fd);
4999         fprintf(stdout,"\nSite analysis report written to: \"%s\"\n",report_name);
5000 }
5001
5002 void LoadTopoData(int max_lon, int min_lon, int max_lat, int min_lat)
5003 {
5004         /* This function loads the SDF files required
5005            to cover the limits of the region specified. */
5006
5007         int     x, y, width, ymin, ymax;
5008
5009         width=ReduceAngle(max_lon-min_lon);
5010
5011         if ((max_lon-min_lon)<=180.0)
5012         {
5013                 for (y=0; y<=width; y++)
5014                         for (x=min_lat; x<=max_lat; x++)
5015                         {
5016                                 ymin=(int)(min_lon+(double)y);
5017
5018                                 while (ymin<0)
5019                                         ymin+=360;
5020
5021                                 while (ymin>=360)
5022                                         ymin-=360;
5023
5024                                 ymax=ymin+1;
5025
5026                                 while (ymax<0)
5027                                         ymax+=360;
5028
5029                                 while (ymax>=360)
5030                                         ymax-=360;
5031
5032                                 sprintf(string,"%d:%d:%d:%d",x, x+1, ymin, ymax);
5033                                 LoadSDF(string);
5034                         }
5035         }
5036
5037         else
5038         {
5039                 for (y=0; y<=width; y++)
5040                         for (x=min_lat; x<=max_lat; x++)
5041                         {
5042                                 ymin=max_lon+y;
5043
5044                                 while (ymin<0)
5045                                         ymin+=360;
5046
5047                                 while (ymin>=360)
5048                                         ymin-=360;
5049
5050                                 ymax=ymin+1;
5051
5052                                 while (ymax<0)
5053                                         ymax+=360;
5054
5055                                 while (ymax>=360)
5056                                         ymax-=360;
5057
5058                                 sprintf(string,"%d:%d:%d:%d",x, x+1, ymin, ymax);
5059                                 LoadSDF(string);
5060                         }
5061         }
5062 }
5063
5064 int LoadPLI(char *filename)
5065 {
5066         int     error=0, max_west, min_west, max_north, min_north;
5067         char    string[80], *pointer=NULL;
5068         double  latitude=0.0, longitude=0.0, azimuth=0.0, elevation=0.0,
5069                 loss=0.0;
5070         FILE    *fd;
5071
5072         fd=fopen(filename,"r");
5073
5074         if (fd!=NULL)
5075         {
5076                 fgets(string,78,fd);
5077                 pointer=strchr(string,';');
5078
5079                 if (pointer!=NULL)
5080                         *pointer=0;
5081
5082                 sscanf(string,"%d, %d",&max_west, &min_west);
5083
5084                 fgets(string,78,fd);
5085                 pointer=strchr(string,';');
5086
5087                 if (pointer!=NULL)
5088                         *pointer=0;
5089
5090                 sscanf(string,"%d, %d",&max_north, &min_north);
5091
5092                 fgets(string,78,fd);
5093                 pointer=strchr(string,';');
5094
5095                 if (pointer!=NULL)
5096                         *pointer=0;
5097
5098                 LoadTopoData(max_west-1, min_west, max_north-1, min_north);
5099
5100                 fprintf(stdout,"\nReading \"%s\"... ",filename);
5101                 fflush(stdout);
5102
5103                 fscanf(fd,"%lf, %lf, %lf, %lf, %lf",&latitude, &longitude, &azimuth, &elevation, &loss);
5104
5105                 while (feof(fd)==0)
5106                 {
5107                         if (loss>225.0)
5108                                 loss=225.0;
5109
5110                         if (loss<75.0)
5111                                 loss=75.0;
5112
5113                         loss-=75.0;
5114                         loss/=10.0;
5115                         loss+=1.0;
5116
5117                         if (loss<=(double)maxdB)
5118                                 OrMask(latitude,longitude,((unsigned char)(loss))<<3);
5119
5120                         fscanf(fd,"%lf, %lf, %lf, %lf, %lf",&latitude, &longitude, &azimuth, &elevation, &loss);
5121                 }
5122
5123                 fclose(fd);
5124
5125                 fprintf(stdout," Done!\n");
5126                 fflush(stdout);
5127         }
5128
5129         else
5130                 error=1;
5131
5132         return error;
5133 }
5134
5135 void WriteKML(struct site source, struct site destination)
5136 {
5137         int     x, y;
5138         char    block, report_name[80];
5139         double  distance, rx_alt, tx_alt, cos_xmtr_angle,
5140                 azimuth, cos_test_angle, test_alt;
5141         FILE    *fd=NULL;
5142
5143         ReadPath(source,destination);
5144
5145         sprintf(report_name,"%s-to-%s.kml",source.name,destination.name);
5146
5147         for (x=0; report_name[x]!=0; x++)
5148                 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
5149                         report_name[x]='_';
5150
5151         fd=fopen(report_name,"w");
5152
5153         fprintf(fd,"<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n");
5154         fprintf(fd,"<kml xmlns=\"http://earth.google.com/kml/2.0\">\n");
5155         fprintf(fd,"<!-- Generated by SPLAT! Version %s -->\n",splat_version);
5156         fprintf(fd,"<Folder>\n");
5157         fprintf(fd,"<name>SPLAT! Path</name>\n");
5158         fprintf(fd,"<open>1</open>\n");
5159         fprintf(fd,"<description>Path Between %s and %s</description>\n",source.name,destination.name);
5160
5161         fprintf(fd,"<Placemark>\n");
5162         fprintf(fd,"    <name>%s</name>\n",source.name);
5163         fprintf(fd,"    <description>\n");
5164         fprintf(fd,"       Transmit Site\n");
5165
5166         if (source.lat>=0.0)
5167                 fprintf(fd,"       <BR>%s North</BR>\n",dec2dms(source.lat));
5168         else
5169                 fprintf(fd,"       <BR>%s South</BR>\n",dec2dms(source.lat));
5170
5171         fprintf(fd,"       <BR>%s West</BR>\n",dec2dms(source.lon));
5172
5173         azimuth=Azimuth(source,destination);
5174         distance=Distance(source,destination);
5175
5176         if (metric)
5177                 fprintf(fd,"       <BR>%.2f km",distance*KM_PER_MILE);
5178         else
5179                 fprintf(fd,"       <BR>%.2f miles",distance);
5180
5181         fprintf(fd," to %s</BR>\n       <BR>toward an azimuth of %.2f%c</BR>\n",destination.name,azimuth,176);
5182
5183         fprintf(fd,"    </description>\n");
5184         fprintf(fd,"    <visibility>1</visibility>\n");
5185         fprintf(fd,"    <Style>\n");
5186         fprintf(fd,"      <IconStyle>\n");
5187         fprintf(fd,"        <Icon>\n");
5188         fprintf(fd,"          <href>root://icons/palette-5.png</href>\n");
5189         fprintf(fd,"          <x>224</x>\n");
5190         fprintf(fd,"          <y>224</y>\n");
5191         fprintf(fd,"          <w>32</w>\n");
5192         fprintf(fd,"          <h>32</h>\n");
5193         fprintf(fd,"        </Icon>\n");
5194         fprintf(fd,"      </IconStyle>\n");
5195         fprintf(fd,"    </Style>\n");
5196         fprintf(fd,"    <Point>\n");
5197         fprintf(fd,"      <extrude>1</extrude>\n");
5198         fprintf(fd,"      <altitudeMode>relativeToGround</altitudeMode>\n");
5199         fprintf(fd,"      <coordinates>%f,%f,30</coordinates>\n",(source.lon<180.0?-source.lon:360.0-source.lon),source.lat);
5200         fprintf(fd,"    </Point>\n");
5201         fprintf(fd,"</Placemark>\n");
5202
5203         fprintf(fd,"<Placemark>\n");
5204         fprintf(fd,"    <name>%s</name>\n",destination.name);
5205         fprintf(fd,"    <description>\n");
5206         fprintf(fd,"       Receive Site\n");
5207
5208         if (destination.lat>=0.0)
5209                 fprintf(fd,"       <BR>%s North</BR>\n",dec2dms(destination.lat));
5210         else
5211                 fprintf(fd,"       <BR>%s South</BR>\n",dec2dms(destination.lat));
5212
5213         fprintf(fd,"       <BR>%s West</BR>\n",dec2dms(destination.lon));
5214
5215
5216         if (metric)
5217                 fprintf(fd,"       <BR>%.2f km",distance*KM_PER_MILE);
5218         else
5219                 fprintf(fd,"       <BR>%.2f miles",distance);
5220
5221         fprintf(fd," to %s</BR>\n       <BR>toward an azimuth of %.2f%c</BR>\n",source.name,Azimuth(destination,source),176);
5222
5223         fprintf(fd,"    </description>\n");
5224         fprintf(fd,"    <visibility>1</visibility>\n");
5225         fprintf(fd,"    <Style>\n");
5226         fprintf(fd,"      <IconStyle>\n");
5227         fprintf(fd,"        <Icon>\n");
5228         fprintf(fd,"          <href>root://icons/palette-5.png</href>\n");
5229         fprintf(fd,"          <x>224</x>\n");
5230         fprintf(fd,"          <y>224</y>\n");
5231         fprintf(fd,"          <w>32</w>\n");
5232         fprintf(fd,"          <h>32</h>\n");
5233         fprintf(fd,"        </Icon>\n");
5234         fprintf(fd,"      </IconStyle>\n");
5235         fprintf(fd,"    </Style>\n");
5236         fprintf(fd,"    <Point>\n");
5237         fprintf(fd,"      <extrude>1</extrude>\n");
5238         fprintf(fd,"      <altitudeMode>relativeToGround</altitudeMode>\n");
5239         fprintf(fd,"      <coordinates>%f,%f,30</coordinates>\n",(destination.lon<180.0?-destination.lon:360.0-destination.lon),destination.lat);
5240         fprintf(fd,"    </Point>\n");
5241         fprintf(fd,"</Placemark>\n");
5242
5243         fprintf(fd,"<Placemark>\n");
5244         fprintf(fd,"<name>Point-to-Point Path</name>\n");
5245         fprintf(fd,"  <visibility>1</visibility>\n");
5246         fprintf(fd,"  <open>0</open>\n");
5247         fprintf(fd,"  <Style>\n");
5248         fprintf(fd,"    <LineStyle>\n");
5249         fprintf(fd,"      <color>7fffffff</color>\n");
5250         fprintf(fd,"    </LineStyle>\n");
5251         fprintf(fd,"    <PolyStyle>\n");
5252         fprintf(fd,"       <color>7fffffff</color>\n");
5253         fprintf(fd,"    </PolyStyle>\n");
5254         fprintf(fd,"  </Style>\n");
5255         fprintf(fd,"  <LineString>\n");
5256         fprintf(fd,"    <extrude>1</extrude>\n");
5257         fprintf(fd,"    <tessellate>1</tessellate>\n");
5258         fprintf(fd,"    <altitudeMode>relativeToGround</altitudeMode>\n");
5259         fprintf(fd,"    <coordinates>\n");
5260
5261         for (x=0; x<path.length; x++)
5262                 fprintf(fd,"      %f,%f,5\n",(path.lon[x]<180.0?-path.lon[x]:360.0-path.lon[x]),path.lat[x]);
5263
5264         fprintf(fd,"    </coordinates>\n");
5265         fprintf(fd,"   </LineString>\n");
5266         fprintf(fd,"</Placemark>\n");
5267
5268         fprintf(fd,"<Placemark>\n");
5269         fprintf(fd,"<name>Line-of-Sight Path</name>\n");
5270         fprintf(fd,"  <visibility>1</visibility>\n");
5271         fprintf(fd,"  <open>0</open>\n");
5272         fprintf(fd,"  <Style>\n");
5273         fprintf(fd,"    <LineStyle>\n");
5274         fprintf(fd,"      <color>ff00ff00</color>\n");
5275         fprintf(fd,"    </LineStyle>\n");
5276         fprintf(fd,"    <PolyStyle>\n");
5277         fprintf(fd,"       <color>7f00ff00</color>\n");
5278         fprintf(fd,"    </PolyStyle>\n");
5279         fprintf(fd,"  </Style>\n");
5280         fprintf(fd,"  <LineString>\n");
5281         fprintf(fd,"    <extrude>1</extrude>\n");
5282         fprintf(fd,"    <tessellate>1</tessellate>\n");
5283         fprintf(fd,"    <altitudeMode>relativeToGround</altitudeMode>\n");
5284         fprintf(fd,"    <coordinates>\n");
5285
5286         /* Walk across the "path", indentifying obstructions along the way */
5287
5288         for (y=0; y<path.length; y++)
5289         {
5290                 distance=5280.0*path.distance[y];
5291                 tx_alt=earthradius+source.alt+path.elevation[0];
5292                 rx_alt=earthradius+destination.alt+path.elevation[y];
5293
5294                 /* Calculate the cosine of the elevation of the
5295                    transmitter as seen at the temp rx point. */
5296
5297                 cos_xmtr_angle=((rx_alt*rx_alt)+(distance*distance)-(tx_alt*tx_alt))/(2.0*rx_alt*distance);
5298
5299                 for (x=y, block=0; x>=0 && block==0; x--)
5300                 {
5301                         distance=5280.0*(path.distance[y]-path.distance[x]);
5302                         test_alt=earthradius+path.elevation[x];
5303
5304                         cos_test_angle=((rx_alt*rx_alt)+(distance*distance)-(test_alt*test_alt))/(2.0*rx_alt*distance);
5305
5306                         /* Compare these two angles to determine if
5307                            an obstruction exists.  Since we're comparing
5308                            the cosines of these angles rather than
5309                            the angles themselves, the following "if"
5310                            statement is reversed from what it would
5311                            be if the actual angles were compared. */
5312
5313                         if (cos_xmtr_angle>cos_test_angle)
5314                                 block=1;
5315                 }
5316
5317                 if (block)
5318                         fprintf(fd,"      %f,%f,-30\n",(path.lon[y]<180.0?-path.lon[y]:360.0-path.lon[y]),path.lat[y]);
5319                 else
5320                         fprintf(fd,"      %f,%f,5\n",(path.lon[y]<180.0?-path.lon[y]:360.0-path.lon[y]),path.lat[y]);
5321         }
5322
5323         fprintf(fd,"    </coordinates>\n");
5324         fprintf(fd,"  </LineString>\n");
5325         fprintf(fd,"</Placemark>\n");
5326
5327         fprintf(fd,"    <LookAt>\n");
5328         fprintf(fd,"      <longitude>%f</longitude>\n",(source.lon<180.0?-source.lon:360.0-source.lon));
5329         fprintf(fd,"      <latitude>%f</latitude>\n",source.lat);
5330         fprintf(fd,"      <range>300.0</range>\n");
5331         fprintf(fd,"      <tilt>45.0</tilt>\n");
5332         fprintf(fd,"      <heading>%f</heading>\n",azimuth);
5333         fprintf(fd,"    </LookAt>\n");
5334
5335         fprintf(fd,"</Folder>\n");
5336         fprintf(fd,"</kml>\n");
5337
5338         fclose(fd);
5339
5340         fprintf(stdout, "KML file written to: \"%s\"\n",report_name);
5341
5342         fflush(stdout);
5343 }
5344
5345 int main(char argc, char *argv[])
5346 {
5347         int             x, y, z=0, min_lat, min_lon, max_lat, max_lon,
5348                         rxlat, rxlon, txlat, txlon, west_min, west_max,
5349                         north_min, north_max;
5350
5351         unsigned char   coverage=0, LRmap=0, ext[20], terrain_plot=0,
5352                         elevation_plot=0, height_plot=0,
5353                         longley_plot=0, cities=0, bfs=0, txsites=0,
5354                         count, report='y', norm=0, topomap=0, geo=0,
5355                         kml=0;
5356
5357         char            mapfile[255], header[80], city_file[5][255],
5358                         elevation_file[255], height_file[255],
5359                         longley_file[255], terrain_file[255],
5360                         string[255], rxfile[255], *env=NULL,
5361                         txfile[255], map=0, boundary_file[5][255],
5362                         udt_file[255], rxsite=0, plo_filename[255],
5363                         pli_filename[255], nf=0;
5364
5365         double          altitude=0.0, altitudeLR=0.0, tx_range=0.0,
5366                         rx_range=0.0, deg_range=0.0, deg_limit,
5367                         deg_range_lon, er_mult, freq=0.0;
5368
5369         struct          site tx_site[4], rx_site;
5370
5371         FILE            *fd;
5372
5373
5374         if (argc==1)
5375         {
5376                 fprintf(stdout,"\n\t\t --==[ SPLAT! v%s Available Options... ]==--\n\n",splat_version);
5377                 fprintf(stdout,"       -t txsite(s).qth (max of 4)\n");
5378                 fprintf(stdout,"       -r rxsite.qth\n");
5379                 fprintf(stdout,"       -c plot coverage of TX(s) with an RX antenna at X feet/meters AGL\n");
5380                 fprintf(stdout,"       -L plot path loss map of TX based on an RX at X feet/meters AGL\n");
5381                 fprintf(stdout,"       -s filename(s) of city/site file(s) to import (5 max)\n");
5382                 fprintf(stdout,"       -b filename(s) of cartographic boundary file(s) to import (max of 5)\n");
5383                 fprintf(stdout,"       -p filename of terrain profile graph to plot\n");
5384                 fprintf(stdout,"       -e filename of terrain elevation graph to plot\n");
5385                 fprintf(stdout,"       -h filename of terrain height graph to plot\n");
5386                 fprintf(stdout,"       -H filename of normalized terrain height graph to plot\n");
5387                 fprintf(stdout,"       -l filename of Longley-Rice graph to plot\n");
5388                 fprintf(stdout,"       -o filename of topographic map to generate (.ppm)\n");
5389                 fprintf(stdout,"       -u filename of user-defined terrain file to import\n");
5390                 fprintf(stdout,"       -d sdf file directory path (overrides path in ~/.splat_path file)\n");
5391                 fprintf(stdout,"       -n no analysis, brief report\n");
5392                 fprintf(stdout,"       -N no analysis, no report\n");
5393                 fprintf(stdout,"       -m earth radius multiplier\n");
5394                 fprintf(stdout,"       -f frequency for Fresnel zone calculation (MHz)\n");
5395                 fprintf(stdout,"       -R modify default range for -c or -L (miles/kilometers)\n");
5396                 fprintf(stdout,"      -db maximum loss contour to display on path loss maps (80-230 dB)\n");
5397                 fprintf(stdout,"      -nf do not plot Fresnel zones in height plots\n");
5398                 fprintf(stdout,"     -plo filename of path-loss output file\n");
5399                 fprintf(stdout,"     -pli filename of path-loss input file\n");
5400                 fprintf(stdout,"     -udt filename of user defined terrain input file\n");
5401                 fprintf(stdout,"     -geo generate an Xastir .geo georeference file (with .ppm output)\n");
5402                 fprintf(stdout,"     -kml generate a Google Earth .kml file (for point-to-point links)\n");
5403                 fprintf(stdout,"  -metric employ metric rather than imperial units for all user I/O\n\n");
5404
5405                 fprintf(stdout,"If that flew by too fast, consider piping the output through 'less':\n");
5406                 fprintf(stdout,"\n\tsplat | less\n\n");
5407                 fprintf(stdout,"Type 'man splat', or see the documentation for more details.\n\n");
5408                 fflush(stdout);
5409                 return 1;
5410         }
5411
5412         y=argc-1;
5413
5414         metric=0;
5415         rxfile[0]=0;
5416         txfile[0]=0;
5417         string[0]=0;
5418         mapfile[0]=0;
5419         elevation_file[0]=0;
5420         terrain_file[0]=0;
5421         sdf_path[0]=0;
5422         udt_file[0]=0;
5423         path.length=0;
5424         LR.frq_mhz=0.0;
5425         rx_site.lat=91.0;
5426         rx_site.lon=361.0;
5427         plo_filename[0]=0;
5428         pli_filename[0]=0;
5429         earthradius=EARTHRADIUS;
5430
5431         sprintf(header,"\n\t\t--==[ Welcome To SPLAT! v%s ]==--\n\n", splat_version);
5432
5433         for (x=0; x<4; x++)
5434         {
5435                 tx_site[x].lat=91.0;
5436                 tx_site[x].lon=361.0;
5437         }
5438
5439         for (x=0; x<MAXSLOTS; x++)
5440         {
5441                 dem[x].min_el=32768;
5442                 dem[x].max_el=-32768;
5443                 dem[x].min_north=90;
5444                 dem[x].max_north=-90;
5445                 dem[x].min_west=360;
5446                 dem[x].max_west=-1;
5447         }
5448
5449         /* Scan for command line arguments */
5450
5451         for (x=1; x<=y; x++)
5452         {
5453                 if (strcmp(argv[x],"-R")==0)
5454                 {
5455                         z=x+1;
5456
5457                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5458                         {
5459                                 sscanf(argv[z],"%lf",&max_range);
5460
5461                                 if (max_range<0.0)
5462                                         max_range=0.0;
5463
5464                                 if (max_range>1000.0)
5465                                         max_range=1000.0;
5466                         }
5467                 }
5468
5469                 if (strcmp(argv[x],"-m")==0)
5470                 {
5471                         z=x+1;
5472
5473                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5474                         {
5475                                 sscanf(argv[z],"%lf",&er_mult);
5476
5477                                 if (er_mult<0.1)
5478                                         er_mult=1.0;
5479
5480                                 if (er_mult>1.0e6)
5481                                         er_mult=1.0e6;
5482
5483                                 earthradius*=er_mult;
5484                         }
5485                 }
5486
5487                 if (strcmp(argv[x],"-o")==0)
5488                 {
5489                         z=x+1;
5490
5491                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5492                                 strncpy(mapfile,argv[z],253);
5493                         map=1;
5494                 }
5495
5496                 if (strcmp(argv[x],"-u")==0)
5497                 {
5498                         z=x+1;
5499
5500                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5501                                 strncpy(udt_file,argv[z],253);
5502                 }
5503
5504                 if (strcmp(argv[x],"-c")==0)
5505                 {
5506                         z=x+1;
5507
5508                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5509                         {
5510                                 sscanf(argv[z],"%lf",&altitude);
5511                                 coverage=1;
5512                         }
5513                 }
5514
5515                 if (strcmp(argv[x],"-db")==0 || strcmp(argv[x],"-dB")==0)
5516                 {
5517                         z=x+1;
5518
5519                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5520                         {
5521                                 sscanf(argv[z],"%d",&maxdB);
5522
5523                                 maxdB=abs(maxdB);
5524
5525                                 if (maxdB<80)
5526                                         maxdB=80;
5527
5528                                 if (maxdB>230)
5529                                         maxdB=230;
5530                         }
5531                 }
5532
5533                 if (strcmp(argv[x],"-p")==0)
5534                 {
5535                         z=x+1;
5536
5537                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5538                         {
5539                                 strncpy(terrain_file,argv[z],253);
5540                                 terrain_plot=1;
5541                         }
5542                 }
5543
5544                 if (strcmp(argv[x],"-e")==0)
5545                 {
5546                         z=x+1;
5547
5548                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5549                         {
5550                                 strncpy(elevation_file,argv[z],253);
5551                                 elevation_plot=1;
5552                         }
5553                 }
5554
5555                 if (strcmp(argv[x],"-h")==0 || strcmp(argv[x],"-H")==0)
5556                 {
5557                         z=x+1;
5558
5559                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5560                         {
5561                                 strncpy(height_file,argv[z],253);
5562                                 height_plot=1;
5563                         }
5564
5565                         if (strcmp(argv[x],"-H")==0)
5566                                 norm=1;
5567                         else
5568                                 norm=0;
5569                 }
5570
5571                 if (strcmp(argv[x],"-n")==0)
5572                 {
5573                         report='n';
5574                         map=1;
5575                 }
5576
5577                 if (strcmp(argv[x],"-N")==0)
5578                 {
5579                         report='N';
5580                         map=1;
5581                 }
5582
5583                 if (strcmp(argv[x],"-metric")==0)
5584                         metric=1;
5585
5586                 if (strcmp(argv[x],"-geo")==0)
5587                         geo=1;
5588
5589                 if (strcmp(argv[x],"-kml")==0)
5590                         kml=1;
5591
5592                 if (strcmp(argv[x],"-nf")==0)
5593                         nf=1;
5594
5595                 if (strcmp(argv[x],"-d")==0)
5596                 {
5597                         z=x+1;
5598
5599                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5600                                 strncpy(sdf_path,argv[z],253);
5601                 }
5602
5603                 if (strcmp(argv[x],"-t")==0)
5604                 {
5605                         /* Read Transmitter Location */
5606
5607                         z=x+1;
5608
5609                         while (z<=y && argv[z][0] && argv[z][0]!='-' && txsites<4)
5610                         {
5611                                 strncpy(txfile,argv[z],253);
5612                                 tx_site[txsites]=LoadQTH(txfile);
5613                                 txsites++;
5614                                 z++;
5615                         }
5616
5617                         z--;
5618                 }
5619
5620                 if (strcmp(argv[x],"-L")==0)
5621                 {
5622                         z=x+1;
5623
5624                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5625                         {
5626                                 sscanf(argv[z],"%lf",&altitudeLR);
5627
5628                                 if (coverage)
5629                                         fprintf(stdout,"c and L are exclusive options, ignoring L.\n");
5630                                 else
5631                                 {
5632                                         LRmap=1;
5633                                         ReadLRParm(txfile);
5634                                 }
5635                         }
5636                 }
5637
5638                 if (strcmp(argv[x],"-l")==0)
5639                 {
5640                         z=x+1;
5641
5642                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5643                         {
5644                                 strncpy(longley_file,argv[z],253);
5645                                 longley_plot=1;
5646                                 /* Doing this twice is harmless */
5647                                 ReadLRParm(txfile);
5648                         }
5649                 }
5650
5651                 if (strcmp(argv[x],"-r")==0)
5652                 {
5653                         /* Read Receiver Location */
5654
5655                         z=x+1;
5656
5657                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5658                         {
5659                                 strncpy(rxfile,argv[z],253);
5660                                 rx_site=LoadQTH(rxfile);
5661                                 rxsite=1;
5662                         }
5663                 }
5664
5665                 if (strcmp(argv[x],"-s")==0)
5666                 {
5667                         /* Read city file(s) */
5668
5669                         z=x+1;
5670
5671                         while (z<=y && argv[z][0] && argv[z][0]!='-' && cities<5)
5672                         {
5673                                 strncpy(city_file[cities],argv[z],253);
5674                                 cities++;
5675                                 z++;
5676                         }
5677
5678                         z--;
5679                 }
5680
5681                 if (strcmp(argv[x],"-b")==0)
5682                 {
5683                         /* Read Boundary File(s) */
5684
5685                         z=x+1;
5686
5687                         while (z<=y && argv[z][0] && argv[z][0]!='-' && bfs<5)
5688                         {
5689                                 strncpy(boundary_file[bfs],argv[z],253);
5690                                 bfs++;
5691                                 z++;
5692                         }
5693
5694                         z--;
5695                 }
5696
5697                 if (strcmp(argv[x],"-f")==0)
5698                 {
5699                         z=x+1;
5700
5701                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5702                         {
5703                                 sscanf(argv[z],"%lf",&freq);
5704
5705                                 if (freq<20)
5706                                         freq=20;
5707
5708                                 if (freq>20e3)
5709                                         freq=20e3;
5710                         }
5711                 }
5712
5713                 if (strcmp(argv[x],"-plo")==0)
5714                 {
5715                         z=x+1;
5716
5717                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5718                                 strncpy(plo_filename,argv[z],253);
5719                 }
5720
5721                 if (strcmp(argv[x],"-pli")==0)
5722                 {
5723                         z=x+1;
5724
5725                         if (z<=y && argv[z][0] && argv[z][0]!='-')
5726                                 strncpy(pli_filename,argv[z],253);
5727                 }
5728         }
5729
5730         /* Perform some error checking on the arguments
5731            and switches parsed from the command-line.
5732            If an error is encountered, print a message
5733            and exit gracefully. */
5734
5735         if (txsites==0)
5736         {
5737                 fprintf(stderr,"\n%c*** ERROR: No transmitter site(s) specified!\n\n",7);
5738                 exit (-1);
5739         }
5740
5741         for (x=0, y=0; x<txsites; x++)
5742         {
5743                 if (tx_site[x].lat==91.0 && tx_site[x].lon==361.0)
5744                 {
5745                         fprintf(stderr,"\n*** ERROR: Transmitter site #%d not found!",x+1);
5746                         y++;
5747                 }
5748         }
5749
5750         if (y)
5751         {
5752                 fprintf(stderr,"%c\n\n",7);
5753                 exit (-1);
5754         }
5755
5756         if ((coverage+LRmap+pli_filename[0])==0 && rx_site.lat==91.0 && rx_site.lon==361.0)
5757         {
5758                 if (max_range!=0.0 && txsites!=0)
5759                 {
5760                         /* Plot topographic map of radius "max_range" */
5761
5762                         map=0;
5763                         topomap=1;
5764                         report='N';
5765                 }
5766
5767                 else
5768                 {
5769                         fprintf(stderr,"\n%c*** ERROR: No receiver site found or specified!\n\n",7);
5770                         exit (-1);
5771                 }
5772         }
5773
5774         /* No major errors were detected.  Whew!  :-) */
5775
5776         /* Adjust input parameters if -metric option is used */
5777
5778         if (metric)
5779         {
5780                 altitudeLR/=METERS_PER_FOOT;    /* meters --> feet */
5781                 max_range/=KM_PER_MILE;         /* kilometers --> miles */
5782                 altitude/=METERS_PER_FOOT;      /* meters --> feet */
5783         }
5784
5785         /* If no SDF path was specified on the command line (-d), check
5786            for a path specified in the $HOME/.splat_path file.  If the
5787            file is not found, then sdf_path[] remains NULL, and the
5788            current working directory is assumed to contain the SDF
5789            files. */
5790
5791         if (sdf_path[0]==0)
5792         {
5793                 env=getenv("HOME");
5794                 sprintf(string,"%s/.splat_path",env);
5795                 fd=fopen(string,"r");
5796
5797                 if (fd!=NULL)
5798                 {
5799                         fgets(string,253,fd);
5800
5801                         /* Remove <CR> and/or <LF> from string */
5802
5803                         for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0 && x<253; x++);
5804                         string[x]=0;
5805
5806                         strncpy(sdf_path,string,253);
5807
5808                         fclose(fd);
5809                 }
5810         }
5811
5812         /* Ensure a trailing '/' is present in sdf_path */
5813
5814         if (sdf_path[0])
5815         {
5816                 x=strlen(sdf_path);
5817
5818                 if (sdf_path[x-1]!='/' && x!=0)
5819                 {
5820                         sdf_path[x]='/';
5821                         sdf_path[x+1]=0;
5822                 }
5823         }
5824
5825         fprintf(stdout,"%s",header);
5826         fflush(stdout);
5827
5828         if (pli_filename[0])
5829         {
5830                 y=LoadPLI(pli_filename);
5831
5832                 for (x=0; x<txsites; x++)
5833                         PlaceMarker(tx_site[x]);
5834
5835                 if (rxsite)
5836                         PlaceMarker(rx_site);
5837
5838                 if (bfs)
5839                 {
5840                         for (x=0; x<bfs; x++)
5841                                 LoadBoundaries(boundary_file[x]);
5842                 }
5843
5844                 if (cities)
5845                 {
5846                         for (x=0; x<cities; x++)
5847                                 LoadCities(city_file[x]);
5848                 }
5849
5850                 WritePPMLR(mapfile,geo);
5851
5852                 exit(0);
5853         }
5854
5855         x=0;
5856         y=0;
5857
5858         min_lat=90;
5859         max_lat=-90;
5860
5861         min_lon=(int)floor(tx_site[0].lon);
5862         max_lon=(int)floor(tx_site[0].lon);
5863
5864         for (y=0, z=0; z<txsites; z++)
5865         {
5866                 txlat=(int)floor(tx_site[z].lat);
5867                 txlon=(int)floor(tx_site[z].lon);
5868
5869                 if (txlat<min_lat)
5870                         min_lat=txlat;
5871
5872                 if (txlat>max_lat)
5873                         max_lat=txlat;
5874
5875                 if (LonDiff(txlon,min_lon)<0.0)
5876                         min_lon=txlon;
5877
5878                 if (LonDiff(txlon,max_lon)>0.0)
5879                         max_lon=txlon;
5880         }
5881
5882         if (rxsite)
5883         {
5884                 rxlat=(int)floor(rx_site.lat);
5885                 rxlon=(int)floor(rx_site.lon);
5886
5887                 if (rxlat<min_lat)
5888                         min_lat=rxlat;
5889
5890                 if (rxlat>max_lat)
5891                         max_lat=rxlat;
5892
5893                 if (LonDiff(rxlon,min_lon)<0.0)
5894                         min_lon=rxlon;
5895
5896                 if (LonDiff(rxlon,max_lon)>0.0)
5897                         max_lon=rxlon;
5898         }
5899
5900         /* Load the required SDF files */
5901
5902         LoadTopoData(max_lon, min_lon, max_lat, min_lat);
5903
5904         if (coverage | LRmap | topomap)
5905         {
5906                 if (LRmap)
5907                         txsites=1;
5908
5909                 for (z=0; z<txsites; z++)
5910                 {
5911                         /* "Ball park" estimates used to load any additional
5912                            SDF files required to conduct this analysis. */
5913
5914                         tx_range=sqrt(1.5*(tx_site[z].alt+GetElevation(tx_site[z])));
5915
5916                         if (LRmap)
5917                                 rx_range=sqrt(1.5*altitudeLR);
5918                         else
5919                                 rx_range=sqrt(1.5*altitude);
5920
5921                         /* deg_range determines the maximum
5922                            amount of topo data we read */
5923
5924                         deg_range=(tx_range+rx_range)/69.0;
5925
5926                         /* max_range sets the maximum size of the
5927                            analysis.  A small, non-zero amount can
5928                            be used to shrink the size of the analysis
5929                            and limit the amount of topo data read by
5930                            SPLAT!  A very large number will only increase
5931                            the width of the analysis, not the size of
5932                            the map. */
5933
5934                         if (max_range==0.0)
5935                                 max_range=tx_range+rx_range;
5936
5937                         if (max_range<(tx_range+rx_range))
5938                                 deg_range=max_range/69.0;
5939
5940                         /* Prevent the demand for a really wide coverage
5941                            from allocating more slots than are available
5942                            in memory. */
5943
5944                         switch (MAXSLOTS)
5945                         {
5946                                 case 2: deg_limit=0.25;
5947                                         break;
5948
5949                                 case 4: deg_limit=0.5;
5950                                         break;
5951
5952                                 case 9: deg_limit=1.0;
5953                                         break;
5954
5955                                 case 16: deg_limit=2.0;
5956                                         break;
5957
5958                                 case 25: deg_limit=3.0;
5959                         }
5960
5961                         if (tx_site[z].lat<70.0)
5962                                 deg_range_lon=deg_range/cos(deg2rad*tx_site[z].lat);
5963                         else
5964                                 deg_range_lon=deg_range/cos(deg2rad*70.0);
5965
5966                         /* Correct for squares in degrees not being square in miles */
5967
5968                         if (deg_range>deg_limit)
5969                                 deg_range=deg_limit;
5970
5971                         if (deg_range_lon>deg_limit)
5972                                 deg_range_lon=deg_limit;
5973
5974
5975                         north_min=(int)floor(tx_site[z].lat-deg_range);
5976                         north_max=(int)floor(tx_site[z].lat+deg_range);
5977
5978                         west_min=(int)floor(tx_site[z].lon-deg_range_lon);
5979
5980                         while (west_min<0)
5981                                 west_min+=360;
5982
5983                         while (west_min>=360)
5984                                 west_min-=360;
5985
5986                         west_max=(int)floor(tx_site[z].lon+deg_range_lon);
5987
5988                         while (west_max<0)
5989                                 west_max+=360;
5990
5991                         while (west_max>=360)
5992                                 west_max-=360;
5993
5994                         if (north_min<min_lat)
5995                                 min_lat=north_min;
5996
5997                         if (north_max>max_lat)
5998                                 max_lat=north_max;
5999
6000                         if (LonDiff(west_min,min_lon)<0.0)
6001                                 min_lon=west_min;
6002
6003                         if (LonDiff(west_max,max_lon)>0.0)
6004                                 max_lon=west_max;
6005                 }
6006
6007                 /* Load any additional SDF files, if required */
6008
6009                 LoadTopoData(max_lon, min_lon, max_lat, min_lat);
6010         }
6011
6012         if (udt_file[0])
6013                 LoadUDT(udt_file);
6014
6015         if (mapfile[0] && topomap==0)
6016                 map=1;
6017
6018         if (freq==0.0 && nf==0)
6019                 freq=LR.frq_mhz;
6020         else
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