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

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