Imported Upstream version 1.1.0
[debian/splat] / splat.cpp
1 /****************************************************************************
2 *                     SPLAT: A Terrain Analysis Program                     *
3 *              Copyright John A. Magliacane, KD2BD 1997-2004                *
4 *                        Last update: 24-Jan-2004                           *
5 *****************************************************************************
6 *                                                                           *
7 * This program is free software; you can redistribute it and/or modify it   *
8 * under the terms of the GNU General Public License as published by the     *
9 * Free Software Foundation; either version 2 of the License or any later    *
10 * version.                                                                  *
11 *                                                                           *
12 * This program is distributed in the hope that it will useful, but WITHOUT  *
13 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or     *
14 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License     *
15 * for more details.                                                         *
16 *                                                                           *
17 *****************************************************************************
18 *                                                                           *
19 *     Extensively modified by J. D. McDonald in Jan. 2004 to include        *
20 *    the Longley-Rice propagation model using C++ code from NTIA/ITS.       *
21 *                                                                           *
22 *              See: http://elbert.its.bldrdoc.gov/itm.html                  *
23 *                                                                           *
24 *****************************************************************************
25  g++ -Wall -O3 -s -lm -lbz2 -fomit-frame-pointer itm.cpp splat.cpp -o splat 
26 *****************************************************************************/
27
28 #include <stdio.h>
29 #include <math.h>
30 #include <stdlib.h>
31 #include <string.h>
32 #include <ctype.h>
33 #include <bzlib.h>
34 #include <unistd.h>
35 #include "fontdata.h"
36 #include "smallfont.h"
37
38 #define GAMMA 2.5
39 #define MAXSLOTS 9
40 #define BZBUFFER 65536
41
42 #if MAXSLOTS==4
43 #define ARRAYSIZE 4950
44 #endif
45
46 #if MAXSLOTS==9
47 #define ARRAYSIZE 10870
48 #endif
49
50 #if MAXSLOTS==16
51 #define ARRAYSIZE 19240
52 #endif
53
54 #if MAXSLOTS==25
55 #define ARRAYSIZE 30025
56 #endif
57
58 char string[255], sdf_path[255], opened=0, *splat_version={"1.1.0"};
59
60 double  TWOPI=6.283185307179586, HALFPI=1.570796326794896,
61         PI=3.141592653589793, deg2rad=1.74532925199e-02,
62         EARTHRADIUS=20902230.97, METERS_PER_MILE=1609.344,
63         METERS_PER_FOOT=0.3048, earthradius, max_range=0.0;
64
65 int     min_north=0, max_north=0, min_west=0, max_west=0,
66         max_elevation=0, min_elevation=0, bzerror;
67
68 struct site { double lat;
69               double lon;
70               double alt;
71               char name[50];
72             };
73
74 struct  { float lat[ARRAYSIZE];
75           float lon[ARRAYSIZE];
76           float elevation[ARRAYSIZE];
77           float distance[ARRAYSIZE];
78           int length;
79         } path;
80
81 struct { 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        } dem[MAXSLOTS];
90
91 struct {
92         double eps_dielect; 
93         double sgm_conductivity; 
94         double eno_ns_surfref;
95         double frq_mhz; 
96         double conf; 
97         double rel;
98         int radio_climate;  
99         int pol;
100        } LR;
101
102 double elev_l[ARRAYSIZE+10];
103
104 void point_to_point(double elev[], double tht_m, double rht_m,
105           double eps_dielect, double sgm_conductivity, double eno_ns_surfref,
106           double frq_mhz, int radio_climate, int pol, double conf,
107           double rel, double &dbloss, char *strmode, int &errnum);
108
109 double arccos(double x, double y)
110 {
111         /* This function implements the arc cosine function,
112            returning a value between 0 and TWOPI. */
113
114         double result=0.0;
115
116         if (y>0.0)
117                 result=acos(x/y);
118
119         if (y<0.0)
120                 result=PI+acos(x/y);
121
122         return result;
123 }
124
125 char *dec2dms(double decimal)
126 {
127         /* Converts decimal degrees to degrees, minutes, seconds,
128            (DMS) and returns the result as a character string. */
129  
130         int degrees, minutes, seconds;
131         double a, b, c, d;
132
133         a=floor(decimal);
134         b=60.0*(decimal-a);
135         c=floor(b);
136         d=60.0*(b-c);
137
138         degrees=(int)a;
139         minutes=(int)c;
140         seconds=(int)d;
141
142         if (seconds<0)
143                 seconds=0;
144
145         if (seconds>59)
146                 seconds=59;
147
148         string[0]=0;
149         sprintf(string,"%d%c %d\' %d\"", degrees, 176, minutes, seconds);
150         return (string);
151 }
152
153 int OrMask(double lat, double lon, int value)
154 {
155         /* Lines, text, markings, and coverage areas are stored in a
156            mask that is combined with topology data when topographic
157            maps are generated by SPLAT!.  This function sets bits in
158            the mask based on the latitude and longitude of the area
159            pointed to. */
160
161         int x, y, indx, minlat, minlon;
162         char found;
163
164         minlat=(int)floor(lat);
165         minlon=(int)floor(lon);
166
167         for (indx=0, found=0; indx<MAXSLOTS && found==0;)
168                 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
169                         found=1;
170                 else
171                         indx++;
172
173         if (found)
174         {
175                 x=(int)(1199.0*(lat-floor(lat)));
176                 y=(int)(1199.0*(lon-floor(lon)));
177
178                 dem[indx].mask[x][y]|=value;
179                 return (dem[indx].mask[x][y]);
180         }
181         else
182                 return -1;
183 }
184
185 int GetMask(double lat, double lon)
186 {
187         /* This function returns the mask bits based on the latitude
188            and longitude given. */
189
190         return (OrMask(lat,lon,0));
191 }
192
193 double GetElevation(struct site location)
194 {
195         /* This function returns the elevation (in feet) of any location
196            represented by the digital elevation model data in memory.
197            Function returns -5000.0 for locations not found in memory. */
198
199         char found;
200         int x, y, indx, minlat, minlon;
201         double elevation;
202
203         elevation=-5000.0;
204
205         minlat=(int)floor(location.lat);
206         minlon=(int)floor(location.lon);
207
208         x=(int)(1199.0*(location.lat-floor(location.lat)));
209         y=(int)(1199.0*(location.lon-floor(location.lon)));
210
211         for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
212         {
213                 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
214                 {
215                         elevation=3.28084*dem[indx].data[x][y];
216                         found=1;
217                 }
218         }
219         
220         return elevation;
221 }
222
223 double Distance(struct site site1, struct site site2)
224 {
225         /* This function returns the great circle distance
226            in miles between any two site locations. */
227
228         double lat1, lon1, lat2, lon2, distance;
229
230         lat1=site1.lat*deg2rad;
231         lon1=site1.lon*deg2rad;
232         lat2=site2.lat*deg2rad;
233         lon2=site2.lon*deg2rad;
234
235         distance=3959.0*acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos((lon1)-(lon2)));
236
237         return distance;
238 }
239
240 double Azimuth(struct site source, struct site destination)
241 {
242         /* This function returns the azimuth (in degrees) to the
243            destination as seen from the location of the source. */
244
245         double dest_lat, dest_lon, src_lat, src_lon,
246                beta, azimuth, diff, num, den, fraction;
247
248         dest_lat=destination.lat*deg2rad;
249         dest_lon=destination.lon*deg2rad;
250
251         src_lat=source.lat*deg2rad;
252         src_lon=source.lon*deg2rad;
253                 
254         /* Calculate Surface Distance */
255
256         beta=acos(sin(src_lat)*sin(dest_lat)+cos(src_lat)*cos(dest_lat)*cos(src_lon-dest_lon));
257
258         /* Calculate Azimuth */
259
260         num=sin(dest_lat)-(sin(src_lat)*cos(beta));
261         den=cos(src_lat)*sin(beta);
262         fraction=num/den;
263
264         /* Trap potential problems in acos() due to rounding */
265
266         if (fraction>=1.0)
267                 fraction=1.0;
268
269         if (fraction<=-1.0)
270                 fraction=-1.0;
271
272         /* Calculate azimuth */
273
274         azimuth=acos(fraction);
275
276         /* Reference it to True North */
277
278         diff=dest_lon-src_lon;
279
280         if (diff<=-PI)
281                 diff+=TWOPI;
282
283         if (diff>=PI)
284                 diff-=TWOPI;
285
286         if (diff>0.0)
287                 azimuth=TWOPI-azimuth;
288
289         return (azimuth/deg2rad);               
290 }
291
292 double ElevationAngle(struct site local, struct site remote)
293 {
294         /* This function returns the angle of elevation (in degrees)
295            of the remote location as seen from the local site.
296            A positive result represents an angle of elevation (uptilt),
297            while a negative result represents an angle of depression
298            (downtilt), as referenced to a normal to the center of
299            the earth. */
300            
301         register double a, b, dx;
302
303         a=GetElevation(remote)+remote.alt+earthradius;
304         b=GetElevation(local)+local.alt+earthradius;
305
306         dx=5280.0*Distance(local,remote);
307
308         /* Apply the Law of Cosines */
309
310         return ((180.0*(acos(((b*b)+(dx*dx)-(a*a))/(2.0*b*dx)))/PI)-90.0);
311 }
312
313 void ReadPath(struct site source, struct site destination)
314 {
315         /* This function generates a sequence of latitude and
316            longitude positions between a source location and
317            a destination along a great circle path, and stores
318            elevation and distance information for points along
319            that path in the "path" structure for later use. */
320
321         double azimuth, distance, lat1, lon1, beta,
322                den, num, lat2, lon2, total_distance;
323         int x1, y1, c;
324         struct site tempsite;
325
326         c=0;
327
328         lat1=source.lat*deg2rad;
329         lon1=source.lon*deg2rad;
330         azimuth=Azimuth(source,destination)*deg2rad;
331         total_distance=Distance(source,destination);
332
333         for (distance=0; distance<=total_distance; distance+=0.04)
334         {
335                 beta=distance/3959.0;
336                 lat2=asin(sin(lat1)*cos(beta)+cos(azimuth)*sin(beta)*cos(lat1));
337                 num=cos(beta)-(sin(lat1)*sin(lat2));
338                 den=cos(lat1)*cos(lat2);
339
340                 if (azimuth==0.0 && (beta>HALFPI-lat1))
341                         lon2=lon1+PI;
342
343                 else if (azimuth==HALFPI && (beta>HALFPI+lat1))
344                         lon2=lon1+PI;
345
346                 else if (fabs(num/den)>1.0)
347                         lon2=lon1;
348
349                 else
350                 {
351                         if ((PI-azimuth)>=0.0)
352                                 lon2=lon1-arccos(num,den);
353                         else
354                                 lon2=lon1+arccos(num,den);
355                 }
356         
357                 while (lon2<0.0)
358                         lon2+=TWOPI;
359
360                 while (lon2>TWOPI)
361                         lon2-=TWOPI;
362  
363                 lat2=lat2/deg2rad;
364                 lon2=lon2/deg2rad;
365
366                 x1=(int)(1199.0*(lat2-floor(lat2)));
367                 y1=(int)(1199.0*(lon2-floor(lon2)));
368
369                 path.lat[c]=lat2;
370                 path.lon[c]=lon2;
371                 tempsite.lat=lat2;
372                 tempsite.lon=lon2;
373                 path.elevation[c]=GetElevation(tempsite);
374                 path.distance[c]=distance;
375                 c++;
376         }
377
378         /* Make sure exact destination point is recorded at path.length-1 */
379
380         x1=(int)(1199.0*(destination.lat-floor(destination.lat)));
381         y1=(int)(1199.0*(destination.lon-floor(destination.lon)));
382
383         path.lat[c]=destination.lat;
384         path.lon[c]=destination.lon;
385         path.elevation[c]=GetElevation(destination);
386         path.distance[c]=total_distance;
387         c++;
388
389         path.length=c;
390 }
391
392 double AverageTerrain(struct site source, double azimuthx, double start_distance, double end_distance)
393 {
394         /* This function returns the average terrain calculated in
395            the direction of "azimuth" (degrees) between "start_distance"
396            and "end_distance" (miles) from the source location.  If
397            the terrain is all water (non-critical error), -5000.0 is
398            returned.  If not enough SDF data has been loaded into
399            memory to complete the survey (critical error), then
400            -9999.0 is returned. */
401  
402         int c, samples, endpoint;
403         double beta, lat1, lon1, lat2, lon2, num, den, azimuth, terrain=0.0;
404         struct site destination;
405
406         lat1=source.lat*deg2rad;
407         lon1=source.lon*deg2rad;
408
409         /* Generate a path of elevations between the source
410            location and the remote location provided. */
411
412         beta=end_distance/3959.0;
413
414         azimuth=deg2rad*azimuthx;
415
416         lat2=asin(sin(lat1)*cos(beta)+cos(azimuth)*sin(beta)*cos(lat1));
417         num=cos(beta)-(sin(lat1)*sin(lat2));
418         den=cos(lat1)*cos(lat2);
419
420         if (azimuth==0.0 && (beta>HALFPI-lat1))
421                 lon2=lon1+PI;
422
423         else if (azimuth==HALFPI && (beta>HALFPI+lat1))
424                 lon2=lon1+PI;
425
426         else if (fabs(num/den)>1.0)
427                 lon2=lon1;
428
429         else
430         {
431                 if ((PI-azimuth)>=0.0)
432                         lon2=lon1-arccos(num,den);
433                 else
434                         lon2=lon1+arccos(num,den);
435         }
436         
437         while (lon2<0.0)
438                 lon2+=TWOPI;
439
440         while (lon2>TWOPI)
441                 lon2-=TWOPI;
442  
443         lat2=lat2/deg2rad;
444         lon2=lon2/deg2rad;
445
446         destination.lat=lat2;
447         destination.lon=lon2;
448
449         /* If SDF data is missing for the endpoint of
450            the radial, then the average terrain cannot
451            be accurately calculated.  Return -9999.0 */
452
453         if (GetElevation(destination)<-4999.0)
454                 return (-9999.0);
455         else
456         {
457                 ReadPath(source,destination);
458
459                 endpoint=path.length;
460
461                 /* Shrink the length of the radial if the
462                    outermost portion is not over U.S. land. */
463
464                 for (c=endpoint-1; c>=0 && path.elevation[c]==0.0; c--);
465
466                 endpoint=c+1;
467
468                 for (c=0, samples=0; c<endpoint; c++)
469                 {
470                         if (path.distance[c]>=start_distance)
471                         {
472                                 terrain+=path.elevation[c];
473                                 samples++;
474                         }
475                 }
476
477                 if (samples==0)
478                         terrain=-5000.0;  /* No land */
479                 else
480                         terrain=(terrain/(double)samples);
481
482                 return terrain;
483         }
484 }
485
486 double haat(struct site antenna)
487 {
488         /* This function returns the antenna's Height Above Average
489            Terrain (HAAT) based on FCC Part 73.313(d).  If a critical
490            error occurs, such as a lack of SDF data to complete the
491            survey, -5000.0 is returned. */
492
493         int azi, c;
494         char error=0;
495         double terrain, avg_terrain, haat, sum=0.0;
496
497         /* Calculate the average terrain between 2 and 10 miles
498            from the antenna site at azimuths of 0, 45, 90, 135,
499            180, 225, 270, and 315 degrees. */
500
501         for (c=0, azi=0; azi<=315 && error==0; azi+=45)
502         {
503                 terrain=AverageTerrain(antenna, (double)azi, 2.0, 10.0);
504
505                 if (terrain<-9998.0)  /* SDF data is missing */
506                         error=1;
507
508                 if (terrain>-4999.0)  /* It's land, not water */
509                 {
510                         sum+=terrain;  /* Sum of averages */
511                         c++;
512                 }
513         }
514
515         if (error)
516                 return -5000.0;
517         else
518         {
519                 avg_terrain=(sum/(double)c);
520                 haat=(antenna.alt+GetElevation(antenna))-avg_terrain;
521                 return haat;
522         }
523 }
524
525 void PlaceMarker(struct site location)
526 {
527         /* This function places text and marker data in the mask array
528            for illustration on topographic maps generated by SPLAT!.
529            By default, SPLAT! centers text information BELOW the marker,
530            but may move it above, to the left, or to the right of the
531            marker depending on how much room is available on the map,
532            or depending on whether the area is already occupied by
533            another marker or label.  If no room or clear space is
534            available on the map to place the marker and its associated
535            text, then the marker and text are not written to the map. */
536  
537         int     a, b, c, byte;
538         char    ok2print, occupied;
539         double  x, y, lat, lon, textx=0.0, texty=0.0, xmin, xmax,
540                 ymin, ymax, p1, p3, p6, p8, p12, p16, p24, label_length;
541
542         xmin=min_north;
543         xmax=max_north;
544         ymin=min_west;
545         ymax=max_west;
546         lat=location.lat;
547         lon=location.lon;
548
549         if (lat<xmax && lat>xmin && lon<ymax && lon>ymin)
550         {
551                 p1=1.0/1200.0;
552                 p3=3.0/1200.0;
553                 p6=6.0/1200.0;
554                 p8=8.0/1200.0;
555                 p12=12.0/1200.0;
556                 p16=16.0/1200.0;
557                 p24=24.0/1200.0;
558                 ok2print=0;
559                 occupied=0;
560
561                 /* Is Marker Position Clear Of Text Or Other Markers? */
562
563                 for (x=lat-p3; (x<=xmax && x>=xmin && x<=lat+p3); x+=p1)
564                         for (y=lon-p3; (y<=ymax && y>=ymin && y<=lon+p3); y+=p1)
565                                 occupied|=(GetMask(x,y)&2);
566
567                 if (occupied==0)
568                 {
569                         /* Determine Where Text Can Be Positioned */
570
571                         /* label_length=length in pixels.
572                            Each character is 8 pixels wide. */
573
574                         label_length=p1*(double)(strlen(location.name)<<3);
575
576                         if (((lon+label_length)<=ymax) && (lon-label_length)>=ymin)
577                         {
578                                 /* Default: Centered Text */
579
580                                 texty=lon+label_length/2.0;
581
582                                 if ((lat-p8)>=p16)
583                                 {
584                                         /* Position Text Below The Marker */
585
586                                         textx=lat-p8;
587
588                                         x=textx;
589                                         y=texty;
590         
591                                         /* Is This Position Clear Of
592                                            Text Or Other Markers? */
593
594                                         for (a=0, occupied=0; a<16; a++)
595                                         {
596                                                 for (b=0; b<(int)strlen(location.name); b++)
597                                                         for (c=0; c<8; c++, y-=p1)
598                                                                 occupied|=(GetMask(x,y)&2);
599                                                 x-=p1;
600                                                 y=texty;
601                                         }
602
603                                         x=textx;
604                                         y=texty;
605
606                                         if (occupied==0)
607                                                 ok2print=1;
608                                 }
609
610                                 else
611                                 {
612                                         /* Position Text Above The Marker */
613
614                                         textx=lat+p24;
615
616                                         x=textx;
617                                         y=texty;
618         
619                                         /* Is This Position Clear Of
620                                            Text Or Other Markers? */
621
622                                         for (a=0, occupied=0; a<16; a++)
623                                         {
624                                                 for (b=0; b<(int)strlen(location.name); b++)
625                                                         for (c=0; c<8; c++, y-=p1)
626                                                                 occupied|=(GetMask(x,y)&2);
627                                                 x-=p1;
628                                                 y=texty;
629                                         }
630
631                                         x=textx;
632                                         y=texty;
633
634                                         if (occupied==0)
635                                                 ok2print=1;
636                                 }
637                         }
638
639                         if (ok2print==0)
640                         {
641                                 if ((lon-label_length)>=ymin)
642                                 {
643                                         /* Position Text To The
644                                            Right Of The Marker */
645
646                                         textx=lat+p6;
647                                         texty=lon-p12;
648
649                                         x=textx;
650                                         y=texty;
651         
652                                         /* Is This Position Clear Of
653                                            Text Or Other Markers? */
654
655                                         for (a=0, occupied=0; a<16; a++)
656                                         {
657                                                 for (b=0; b<(int)strlen(location.name); b++)
658                                                         for (c=0; c<8; c++, y-=p1)
659                                                                 occupied|=(GetMask(x,y)&2);
660                                                 x-=p1;
661                                                 y=texty;
662                                         }
663
664                                         x=textx;
665                                         y=texty;
666
667                                         if (occupied==0)
668                                                 ok2print=1;
669                                 }
670
671                                 else
672                                 {
673                                         /* Position Text To The
674                                            Left Of The Marker */
675
676                                         textx=lat+p6;
677                                         texty=lon+p8+(label_length);
678
679                                         x=textx;
680                                         y=texty;
681         
682                                         /* Is This Position Clear Of
683                                            Text Or Other Markers? */
684
685                                         for (a=0, occupied=0; a<16; a++)
686                                         {
687                                                 for (b=0; b<(int)strlen(location.name); b++)
688                                                         for (c=0; c<8; c++, y-=p1)
689                                                                 occupied|=(GetMask(x,y)&2);
690                                                 x-=p1;
691                                                 y=texty;
692                                         }
693
694                                         x=textx;
695                                         y=texty;
696
697                                         if (occupied==0)
698                                                 ok2print=1;
699                                 }
700                         }
701
702                         /* textx and texty contain the latitude and longitude
703                            coordinates that describe the placement of the text
704                            on the map. */
705         
706                         if (ok2print && textx!=0.0 && texty!=0.0)
707                         {
708                                 /* Draw Text */
709
710                                 x=textx;
711                                 y=texty;
712
713                                 for (a=0; a<16 && ok2print; a++)
714                                 {
715                                         for (b=0; b<(int)strlen(location.name); b++)
716                                         {
717                                                 byte=fontdata[16*(location.name[b])+a];
718
719                                                 for (c=128; c>0; c=c>>1, y-=p1)
720                                                         if (byte&c)
721                                                                 OrMask(x,y,2);
722                                         }
723                                         x-=p1;
724                                         y=texty;
725                                 }
726         
727                                 /* Draw Square Marker Centered
728                                    On Location Specified */
729         
730                                 for (x=lat-p3; (x<=xmax && x>=xmin && x<=lat+p3); x+=p1)
731                                         for (y=lon-p3; (y<=ymax && y>=ymin && y<=lon+p3); y+=p1)
732                                                 OrMask(x,y,2);
733                         }
734                 }
735         }
736 }
737
738 double ReadBearing(char *input)
739 {
740         /* This function takes numeric input in the form of a character
741            string, and returns an equivalent bearing in degrees as a
742            decimal number (double).  The input may either be expressed
743            in decimal format (40.139722) or degree, minute, second
744            format (40 08 23).  This function also safely handles
745            extra spaces found either leading, trailing, or
746            embedded within the numbers expressed in the
747            input string.  Decimal seconds are permitted. */
748  
749         double bearing=0.0;
750         char string[20];
751         int a, b, length, degrees, minutes;
752         double seconds;
753
754         /* Copy "input" to "string", and ignore any extra
755            spaces that might be present in the process. */
756
757         string[0]=0;
758         length=strlen(input);
759
760         for (a=0, b=0; a<length && a<18; a++)
761         {
762                 if ((input[a]!=32 && input[a]!='\n') || (input[a]==32 && input[a+1]!=32 && b!=0))
763                 {
764                         string[b]=input[a];
765                         b++;
766                 }        
767         }
768
769         string[b]=0;
770
771         /* Count number of spaces in the clean string. */
772
773         length=strlen(string);
774
775         for (a=0, b=0; a<length; a++)
776                 if (string[a]==32)
777                         b++;
778
779         if (b==0)  /* Decimal Format (40.139722) */
780                 sscanf(string,"%lf",&bearing);
781
782         if (b==2)  /* Degree, Minute, Second Format (40 08 23) */
783         {
784                 sscanf(string,"%d %d %lf",&degrees, &minutes, &seconds);
785                 bearing=(double)degrees+((double)minutes/60)+(seconds/3600);
786         }
787
788         /* Anything else returns a 0.0 */
789
790         if (bearing>360.0 || bearing<0.0)
791                 bearing=0.0;
792
793         return bearing;
794 }
795
796 struct site LoadQTH(char *filename)
797 {
798         /* This function reads SPLAT! .qth (site location) files.
799            The latitude and longitude may be expressed either in
800            decimal degrees, or in degree, minute, second format.
801            Antenna height is assumed to be expressed in feet above
802            ground level (AGL), unless followed by the letter 'M',
803            or 'm', or by the word "meters" or "Meters", in which
804            case meters is assumed, and is handled accordingly. */
805
806         int x;
807         char string[50], qthfile[255];
808         struct site tempsite;
809         FILE *fd=NULL;
810
811         for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
812                 qthfile[x]=filename[x];
813
814         qthfile[x]='.';
815         qthfile[x+1]='q';
816         qthfile[x+2]='t';
817         qthfile[x+3]='h';
818         qthfile[x+4]=0;
819
820         tempsite.lat=0.0;
821         tempsite.lon=0.0;
822         tempsite.alt=0.0;
823         tempsite.name[0]=0;
824
825         fd=fopen(qthfile,"r");
826
827         if (fd!=NULL)
828         {
829                 /* Site Name */
830                 fgets(string,49,fd);
831
832                 /* Strip <CR> and/or <LF> from end of site name */
833
834                 for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0; tempsite.name[x]=string[x], x++);
835
836                 tempsite.name[x]=0;
837
838                 /* Site Latitude */
839                 fgets(string,49,fd);
840                 tempsite.lat=ReadBearing(string);
841
842                 /* Site Longitude */
843                 fgets(string,49,fd);
844                 tempsite.lon=ReadBearing(string);
845
846                 /* Antenna Height */
847                 fgets(string,49,fd);
848                 fclose(fd);
849
850                 /* Remove <CR> and/or <LF> from antenna height string */
851                 for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0; x++);
852
853                 string[x]=0;
854
855                 /* Antenna height may either be in feet or meters.
856                    If the letter 'M' or 'm' is discovered in
857                    the string, then this is an indication that
858                    the value given is expressed in meters, and
859                    must be converted to feet before exiting. */
860
861                 for (x=0; string[x]!='M' && string[x]!='m' && string[x]!=0 && x<48; x++);
862                 if (string[x]=='M' || string[x]=='m')
863                 {
864                         string[x]=0;
865                         sscanf(string,"%lf",&tempsite.alt);
866                         tempsite.alt*=3.28084;
867                 }
868
869                 else
870                 {
871                         string[x]=0;
872                         sscanf(string,"%lf",&tempsite.alt);
873                 }
874         }
875
876         return tempsite;
877 }
878
879 int LoadSDF_SDF(char *name)
880 {
881         /* This function reads uncompressed SPLAT Data Files (.sdf)
882            containing digital elevation model data into memory.
883            Elevation data, maximum and minimum elevations, and
884            quadrangle limits are stored in the first available
885            dem[] structure. */
886
887         int x, y, data, indx, minlat, minlon, maxlat, maxlon;
888         char found, free_slot=0, line[20], sdf_file[255], path_plus_name[255];
889         FILE *fd;
890
891         for (x=0; name[x]!='.' && name[x]!=0 && x<250; x++)
892                 sdf_file[x]=name[x];
893
894         sdf_file[x]=0;
895
896         /* Parse filename for minimum latitude and longitude values */
897
898         sscanf(sdf_file,"%d:%d:%d:%d",&minlat,&maxlat,&minlon,&maxlon);
899
900         sdf_file[x]='.';
901         sdf_file[x+1]='s';
902         sdf_file[x+2]='d';
903         sdf_file[x+3]='f';
904         sdf_file[x+4]=0;
905
906         /* Is it already in memory? */
907
908         for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
909         {
910                 if (minlat!=0 && minlon!=0)
911                 {
912                         if (minlat==dem[indx].min_north && minlon==dem[indx].min_west && maxlat==dem[indx].max_north && maxlon==dem[indx].max_west)
913                                 found=1;
914                 }
915         }
916
917         /* Is room available to load it? */
918
919         if (found==0)
920         {       
921                 for (indx=0, free_slot=0; indx<MAXSLOTS && free_slot==0; indx++)
922                         if (dem[indx].max_north==0 && dem[indx].max_west==0)
923                                 free_slot=1;
924         }
925
926         indx--;
927
928         if (free_slot && found==0 && indx>=0 && indx<MAXSLOTS && minlat!=0 && minlon!=0)
929         {
930                 strncpy(path_plus_name,sdf_path,255);
931                 strncat(path_plus_name,sdf_file,255);
932
933                 fd=fopen(path_plus_name,"rb");
934
935                 if (fd!=NULL)
936                 {
937                         fprintf(stdout,"Loading \"%s\" into slot %d...",path_plus_name,indx+1);
938                         fflush(stdout);
939
940                         fgets(line,19,fd);
941                         sscanf(line,"%d",&dem[indx].max_west);
942
943                         fgets(line,19,fd);
944                         sscanf(line,"%d",&dem[indx].min_north);
945
946                         fgets(line,19,fd);
947                         sscanf(line,"%d",&dem[indx].min_west);
948
949                         fgets(line,19,fd);
950                         sscanf(line,"%d",&dem[indx].max_north);
951
952                         for (x=0; x<1200; x++)
953                                 for (y=0; y<1200; y++)
954                                 {
955                                         fgets(line,19,fd);
956                                         sscanf(line,"%d",&data);
957                                         dem[indx].data[x][y]=data;
958
959                                         if (data>dem[indx].max_el)
960                                                 dem[indx].max_el=data;
961
962                                         if (dem[indx].min_el==0)
963                                                 dem[indx].min_el=data;
964                                         else
965                                         {
966                                                  if (data<dem[indx].min_el)
967                                                         dem[indx].min_el=data;
968                                         }
969                                 }
970
971                         fclose(fd);
972
973                         if (min_elevation==0)
974                                 min_elevation=dem[indx].min_el;
975
976                         else
977                         {
978                                 if (dem[indx].min_el<min_elevation)
979                                         min_elevation=dem[indx].min_el;
980                         }
981
982                         if (dem[indx].max_el>max_elevation)
983                                 max_elevation=dem[indx].max_el;
984
985                         if (dem[indx].max_north>max_north)
986                                 max_north=dem[indx].max_north;
987
988                         if (dem[indx].max_west>max_west)
989                                 max_west=dem[indx].max_west;
990
991                         if (min_north==0)
992                                 min_north=dem[indx].min_north;
993                         else
994                         {
995                                 if (dem[indx].min_north<min_north)
996                                         min_north=dem[indx].min_north;
997                         }
998
999                         if (min_west==0)
1000                                 min_west=dem[indx].min_west;
1001                         else
1002                         {
1003                                 if (dem[indx].min_west<min_west)
1004                                         min_west=dem[indx].min_west;
1005                         }
1006
1007                         fprintf(stdout," Done!\n");
1008                         fflush(stdout);
1009                         return 1;
1010                 }
1011
1012                 else
1013                         return -1;
1014         }
1015
1016         else
1017                 return 0;
1018 }
1019
1020 char *BZfgets(BZFILE *bzfd, unsigned length)
1021 {
1022         /* This function returns at most one less than 'length' number
1023            of characters from a bz2 compressed file whose file descriptor
1024            is pointed to by *bzfd.  In operation, a buffer is filled with
1025            uncompressed data (size = BZBUFFER), which is then parsed
1026            and doled out as NULL terminated character strings every time
1027            this function is invoked.  A NULL string indicates an EOF
1028            or error condition. */
1029
1030         static int x, y, nBuf;
1031         static char buffer[BZBUFFER+1], output[BZBUFFER+1];
1032         char done=0;
1033
1034         if (opened!=1 && bzerror==BZ_OK)
1035         {
1036                 /* First time through.  Initialize everything! */
1037
1038                 x=0;
1039                 y=0;
1040                 nBuf=0;
1041                 opened=1;
1042                 output[0]=0;
1043         }
1044
1045         do
1046         {
1047                 if (x==nBuf && bzerror!=BZ_STREAM_END && bzerror==BZ_OK && opened)
1048                 {
1049                         /* Uncompress data into a static buffer */
1050
1051                         nBuf=BZ2_bzRead(&bzerror, bzfd, buffer, BZBUFFER);
1052                         buffer[nBuf]=0;
1053                         x=0;
1054                 }
1055
1056                 /* Build a string from buffer contents */
1057
1058                 output[y]=buffer[x];
1059
1060                 if (output[y]=='\n' || output[y]==0 || y==(int)length-1)
1061                 {
1062                         output[y+1]=0;
1063                         done=1;
1064                         y=0;
1065                 }
1066
1067                 else
1068                         y++;
1069                 x++;
1070
1071         } while (done==0);
1072
1073         if (output[0]==0)
1074                 opened=0;
1075
1076         return (output);
1077 }
1078
1079 int LoadSDF_BZ(char *name)
1080 {
1081         /* This function reads .bz2 compressed SPLAT Data Files containing
1082            digital elevation model data into memory.  Elevation data,
1083            maximum and minimum elevations, and quadrangle limits are
1084            stored in the first available dem[] structure. */
1085
1086         int x, y, data, indx, minlat, minlon, maxlat, maxlon;
1087         char found, free_slot=0, sdf_file[255], path_plus_name[255];
1088         FILE *fd;
1089         BZFILE *bzfd;
1090
1091         for (x=0; name[x]!='.' && name[x]!=0 && x<247; x++)
1092                 sdf_file[x]=name[x];
1093
1094         sdf_file[x]=0;
1095
1096         /* Parse sdf_file name for minimum latitude and longitude values */
1097
1098         sscanf(sdf_file,"%d:%d:%d:%d",&minlat,&maxlat,&minlon,&maxlon);
1099
1100         sdf_file[x]='.';
1101         sdf_file[x+1]='s';
1102         sdf_file[x+2]='d';
1103         sdf_file[x+3]='f';
1104         sdf_file[x+4]='.';
1105         sdf_file[x+5]='b';
1106         sdf_file[x+6]='z';
1107         sdf_file[x+7]='2';
1108         sdf_file[x+8]=0;
1109
1110         /* Is it already in memory? */
1111
1112         for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
1113         {
1114                 if (minlat!=0 && minlon!=0)
1115                 {
1116                         if (minlat==dem[indx].min_north && minlon==dem[indx].min_west && maxlat==dem[indx].max_north && maxlon==dem[indx].max_west)
1117                                 found=1;
1118                 }
1119         }
1120
1121         /* Is room available to load it? */
1122
1123         if (found==0)
1124         {       
1125                 for (indx=0, free_slot=0; indx<MAXSLOTS && free_slot==0; indx++)
1126                         if (dem[indx].max_north==0 && dem[indx].max_west==0)
1127                                 free_slot=1;
1128         }
1129
1130         indx--;
1131
1132         if (free_slot && found==0 && indx>=0 && indx<MAXSLOTS && minlat!=0 && minlon!=0)
1133         {
1134                 strncpy(path_plus_name,sdf_path,255);
1135                 strncat(path_plus_name,sdf_file,255);
1136
1137                 fd=fopen(path_plus_name,"rb");
1138                 bzfd=BZ2_bzReadOpen(&bzerror,fd,0,0,NULL,0);
1139
1140                 if (fd!=NULL && bzerror==BZ_OK)
1141                 {
1142                         fprintf(stdout,"Loading \"%s\" into slot %d...",path_plus_name,indx+1);
1143                         fflush(stdout);
1144
1145                         sscanf(BZfgets(bzfd,255),"%d",&dem[indx].max_west);
1146                         sscanf(BZfgets(bzfd,255),"%d",&dem[indx].min_north);
1147                         sscanf(BZfgets(bzfd,255),"%d",&dem[indx].min_west);
1148                         sscanf(BZfgets(bzfd,255),"%d",&dem[indx].max_north);
1149         
1150                         for (x=0; x<1200; x++)
1151                                 for (y=0; y<1200; y++)
1152                                 {
1153                                         sscanf(BZfgets(bzfd,20),"%d",&data);
1154                                         dem[indx].data[x][y]=data;
1155
1156                                         if (data>dem[indx].max_el)
1157                                                 dem[indx].max_el=data;
1158
1159                                         if (dem[indx].min_el==0)
1160                                                 dem[indx].min_el=data;
1161                                         else
1162                                         {
1163                                                  if (data<dem[indx].min_el)
1164                                                         dem[indx].min_el=data;
1165                                         }
1166                                 }
1167
1168                         fclose(fd);
1169
1170                         BZ2_bzReadClose(&bzerror,bzfd);
1171
1172                         if (min_elevation==0)
1173                                 min_elevation=dem[indx].min_el;
1174
1175                         else
1176                         {
1177                                 if (dem[indx].min_el<min_elevation)
1178                                         min_elevation=dem[indx].min_el;
1179                         }
1180         
1181                         if (dem[indx].max_el>max_elevation)
1182                                 max_elevation=dem[indx].max_el;
1183
1184                         if (dem[indx].max_north>max_north)
1185                                 max_north=dem[indx].max_north;
1186
1187                         if (dem[indx].max_west>max_west)
1188                                 max_west=dem[indx].max_west;
1189
1190                         if (min_north==0)
1191                                 min_north=dem[indx].min_north;
1192                         else
1193                         {
1194                                 if (dem[indx].min_north<min_north)
1195                                         min_north=dem[indx].min_north;
1196                         }
1197
1198                         if (min_west==0)
1199                                 min_west=dem[indx].min_west;
1200                         else
1201                         {
1202                                 if (dem[indx].min_west<min_west)
1203                                         min_west=dem[indx].min_west;
1204                         }
1205
1206                         fprintf(stdout," Done!\n");
1207                         fflush(stdout);
1208                         return 1;
1209                 }
1210                 else
1211                         return -1;
1212         }
1213         else
1214                 return 0;
1215 }
1216
1217 char LoadSDF(char *name)
1218 {
1219         /* This function loads the requested SDF file from the filesystem.
1220            It first tries to invoke the LoadSDF_SDF() function to load an
1221            uncompressed SDF file (since uncompressed files load slightly
1222            faster).  If that attempt fails, then it tries to load a
1223            compressed SDF file by invoking the LoadSDF_BZ() function.
1224            If that fails, then we can assume that no elevation data
1225            exists for the region requested, and that the region
1226            requested must be entirely over water. */
1227
1228         int x, y, indx, minlat, minlon, maxlat, maxlon;
1229         char found, free_slot=0;
1230         int  return_value=-1;
1231
1232         /* Try to load an uncompressed SDF first. */
1233
1234         return_value=LoadSDF_SDF(name);
1235
1236         /* If that fails, try loading a compressed SDF. */
1237
1238         if (return_value==0 || return_value==-1)
1239                 return_value=LoadSDF_BZ(name);
1240
1241         /* If neither format can be found, then assume the area is water. */
1242
1243         if (return_value==0 || return_value==-1)
1244         {
1245                 /* Parse SDF name for minimum latitude and longitude values */
1246
1247                 sscanf(name,"%d:%d:%d:%d",&minlat,&maxlat,&minlon,&maxlon);
1248
1249                 /* Is it already in memory? */
1250
1251                 for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
1252                 {
1253                         if (minlat!=0 && minlon!=0)
1254                         {
1255                                 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west && maxlat==dem[indx].max_north && maxlon==dem[indx].max_west)
1256                                         found=1;
1257                         }
1258                 }
1259
1260                 /* Is room available to load it? */
1261
1262                 if (found==0)
1263                 {       
1264                         for (indx=0, free_slot=0; indx<MAXSLOTS && free_slot==0; indx++)
1265                                 if (dem[indx].max_north==0 && dem[indx].max_west==0)
1266                                         free_slot=1;
1267                 }
1268
1269                 indx--;
1270
1271                 if (free_slot && found==0 && indx>=0 && indx<MAXSLOTS && minlat!=0 && minlon!=0)
1272                 {
1273                         fprintf(stdout,"Region  \"%s\" assumed as sea-level into slot %d...",name,indx+1);
1274                         fflush(stdout);
1275
1276                         dem[indx].max_west=maxlon;
1277                         dem[indx].min_north=minlat;
1278                         dem[indx].min_west=minlon;
1279                         dem[indx].max_north=maxlat;
1280
1281                         for (x=0; x<1200; x++)
1282                                 for (y=0; y<1200; y++)
1283                                 {
1284                                         dem[indx].data[x][y]=0;
1285
1286                                         if (dem[indx].min_el>0)
1287                                                 dem[indx].min_el=0;
1288                                 }
1289
1290                         if (dem[indx].min_el<min_elevation)
1291                                 min_elevation=dem[indx].min_el;
1292
1293                         if (dem[indx].max_el>max_elevation)
1294                                 max_elevation=dem[indx].max_el;
1295
1296                         if (dem[indx].max_north>max_north)
1297                                 max_north=dem[indx].max_north;
1298
1299                         if (dem[indx].max_west>max_west)
1300                                 max_west=dem[indx].max_west;
1301
1302                         if (min_north==0)
1303                                 min_north=dem[indx].min_north;
1304                         else
1305                         {
1306                                 if (dem[indx].min_north<min_north)
1307                                         min_north=dem[indx].min_north;
1308                         }
1309
1310                         if (min_west==0)
1311                                 min_west=dem[indx].min_west;
1312                         else
1313                         {
1314                                 if (dem[indx].min_west<min_west)
1315                                         min_west=dem[indx].min_west;
1316                         }
1317
1318                         fprintf(stdout," Done!\n");
1319                         fflush(stdout);
1320
1321                         return_value=1;
1322                 }
1323         }
1324
1325         return return_value;
1326 }
1327
1328 void LoadCities(char *filename)
1329 {
1330         /* This function reads SPLAT! city/site files, and plots
1331            the locations and names of the cities and site locations
1332            read on topographic maps generated by SPLAT! */
1333
1334         int x, y, z;
1335         char input[80], str[3][80];
1336         struct site city_site;
1337         FILE *fd=NULL;
1338
1339         fd=fopen(filename,"r");
1340
1341         if (fd!=NULL)
1342         {
1343                 fgets(input,78,fd);
1344
1345                 fprintf(stdout,"Reading \"%s\"... ",filename);
1346                 fflush(stdout);
1347
1348                 while (fd!=NULL && feof(fd)==0)
1349                 {
1350                         /* Parse line for name, latitude, and longitude */
1351
1352                         for (x=0, y=0, z=0; x<78 && input[x]!=0 && z<3; x++)
1353                         {
1354                                 if (input[x]!=',' && y<78)
1355                                 {
1356                                         str[z][y]=input[x];
1357                                         y++;
1358                                 }
1359
1360                                 else
1361                                 {
1362                                         str[z][y]=0;
1363                                         z++;
1364                                         y=0;
1365                                 }
1366                         }
1367
1368                         strncpy(city_site.name,str[0],49);
1369                         city_site.lat=ReadBearing(str[1]);
1370                         city_site.lon=ReadBearing(str[2]);
1371                         city_site.alt=0.0;
1372
1373                         PlaceMarker(city_site);
1374
1375                         fgets(input,78,fd);
1376                 }
1377
1378                 fclose(fd);
1379                 fprintf(stdout,"Done!\n");
1380                 fflush(stdout);
1381         }
1382         else
1383                 fprintf(stderr,"*** ERROR: \"%s\": not found!\n",filename);
1384 }
1385
1386 void LoadBoundaries(char *filename)
1387 {
1388         /* This function reads Cartographic Boundary Files available from
1389            the U.S. Census Bureau, and plots the data contained in those
1390            files on the PPM Map generated by SPLAT!.  Such files contain
1391            the coordinates that describe the boundaries of cities,
1392            counties, and states. */
1393
1394         int x;
1395         double lat0, lon0, lat1, lon1;
1396         char string[80];
1397         struct site source, destination;
1398         FILE *fd=NULL;
1399
1400         fd=fopen(filename,"r");
1401
1402         if (fd!=NULL)
1403         {
1404                 fgets(string,78,fd);
1405
1406                 fprintf(stdout,"Reading \"%s\"... ",filename);
1407                 fflush(stdout);
1408
1409                 do
1410                 {
1411                         fgets(string,78,fd);
1412                         sscanf(string,"%lf %lf", &lon0, &lat0);
1413                         fgets(string,78,fd);
1414
1415                         do
1416                         {
1417                                 sscanf(string,"%lf %lf", &lon1, &lat1);
1418                         
1419                                 lon0=fabs(lon0);
1420                                 lon1=fabs(lon1);
1421
1422                                 source.lat=lat0;
1423                                 source.lon=lon0;
1424                                 destination.lat=lat1;
1425                                 destination.lon=lon1;
1426
1427                                 ReadPath(source,destination);
1428
1429                                 for (x=0; x<path.length; x++)
1430                                         OrMask(path.lat[x],path.lon[x],4);
1431
1432                                 lat0=lat1;
1433                                 lon0=lon1;
1434
1435                                 fgets(string,78,fd);
1436
1437                         } while (strncmp(string,"END",3)!=0 && feof(fd)==0);
1438
1439                         fgets(string,78,fd);
1440
1441                 } while (strncmp(string,"END",3)!=0 && feof(fd)==0);
1442
1443                 fclose(fd);
1444
1445                 fprintf(stdout,"Done!\n");
1446                 fflush(stdout);
1447         }
1448         else
1449                 fprintf(stderr,"*** ERROR: \"%s\": not found!\n",filename);
1450 }
1451
1452 void ReadLRParm(char *txsite_filename)
1453 {
1454         /* This function reads Longley-Rice parameter data for the
1455            transmitter site.  The file name is the same as the txsite,
1456            except the filename extension is .lrp.  If the needed file
1457            is not found, then the file "splat.lrp" is read from the
1458            current working directory.  Failure to load this file will
1459            result in the default parameters hard coded into this
1460            being used and written to "splat.lrp". */
1461
1462         double din;
1463         char filename[255], lookup[256], string[80];
1464         int iin, ok=0, x;
1465         FILE *fd=NULL, *outfile=NULL;
1466
1467         /* Default parameters in case things go bad */
1468
1469         LR.eps_dielect=15.0;
1470         LR.sgm_conductivity=0.005;
1471         LR.eno_ns_surfref=301.0;
1472         LR.frq_mhz=300.0;
1473         LR.radio_climate=5;
1474         LR.pol=0;
1475         LR.conf=0.50;
1476         LR.rel=0.50;
1477
1478         /* Modify txsite filename to one with a .lrp extension. */
1479
1480         strncpy(filename,txsite_filename,255);
1481
1482         for (x=0; filename[x]!='.' && filename[x]!=0 && filename[x]!='\n' && x<249; x++);
1483
1484         filename[x]='.';
1485         filename[x+1]='l';
1486         filename[x+2]='r';
1487         filename[x+3]='p';
1488         filename[x+4]=0;
1489
1490         /* Small lookup table to parse file, pass
1491            numeric data, and ignore comments. */
1492
1493         for (x=0; x<=255; x++)
1494                 lookup[x]=0;
1495
1496         /* Valid characters */
1497
1498         for (x=48; x<=57; x++)
1499                 lookup[x]=x;
1500
1501         lookup[32]=32;
1502         lookup[46]='.';
1503
1504         fd=fopen(filename,"r");
1505
1506         if (fd==NULL)
1507         {
1508                 /* Load default "splat.lrp" file */
1509
1510                 strncpy(filename,"splat.lrp\0",10);
1511                 fd=fopen(filename,"r");
1512         }
1513
1514         if (fd!=NULL)
1515         {
1516                 fgets(string,80,fd);
1517
1518                 for (x=0; lookup[(int)string[x]] && x<20; x++)
1519                         string[x]=lookup[(int)string[x]];
1520
1521                 string[x]=0;
1522
1523                 ok=sscanf(string,"%lf", &din);
1524
1525                 if (ok)
1526                 {
1527                         LR.eps_dielect=din;
1528
1529                         fgets(string,80,fd);
1530
1531                         for (x=0; lookup[(int)string[x]] && x<20; x++)
1532                                 string[x]=lookup[(int)string[x]];
1533
1534                         string[x]=0;
1535
1536                         ok=sscanf(string,"%lf", &din);
1537                 }
1538
1539                 if (ok)
1540                 {
1541                         LR.sgm_conductivity=din;
1542
1543                         fgets(string,80,fd);
1544
1545                         for (x=0; lookup[(int)string[x]] && x<20; x++)
1546                                 string[x]=lookup[(int)string[x]];
1547
1548                         string[x]=0;
1549
1550                         ok=sscanf(string,"%lf", &din);
1551                 }
1552
1553                 if (ok)
1554                 {
1555                         LR.eno_ns_surfref=din;
1556
1557                         fgets(string,80,fd);
1558
1559                         for (x=0; lookup[(int)string[x]] && x<20; x++)
1560                                 string[x]=lookup[(int)string[x]];
1561
1562                         string[x]=0;
1563
1564                         ok=sscanf(string,"%lf", &din);
1565                 }
1566
1567                 if (ok)
1568                 {
1569                         LR.frq_mhz=din;
1570
1571                         fgets(string,80,fd);
1572
1573                         for (x=0; lookup[(int)string[x]] && x<20; x++)
1574                                 string[x]=lookup[(int)string[x]];
1575
1576                         string[x]=0;
1577
1578                         ok=sscanf(string,"%d", &iin);
1579                 }
1580
1581                 if (ok)
1582                 {
1583                         LR.radio_climate=iin;
1584
1585                         fgets(string,80,fd);
1586
1587                         for (x=0; lookup[(int)string[x]] && x<20; x++)
1588                                 string[x]=lookup[(int)string[x]];
1589
1590                         string[x]=0;
1591
1592                         ok=sscanf(string,"%d", &iin);
1593                 }
1594
1595                 if (ok)
1596                 {
1597                         LR.pol=iin;
1598
1599                         fgets(string,80,fd);
1600
1601                         for (x=0; lookup[(int)string[x]] && x<20; x++)
1602                                 string[x]=lookup[(int)string[x]];
1603
1604                         string[x]=0;
1605
1606                         ok=sscanf(string,"%lf", &din);
1607                 }
1608
1609                 if (ok)
1610                 {
1611                         LR.conf=din;
1612
1613                         fgets(string,80,fd);
1614
1615                         for (x=0; lookup[(int)string[x]] && x<20; x++)
1616                                 string[x]=lookup[(int)string[x]];
1617
1618                         string[x]=0;
1619
1620                         ok=sscanf(string,"%lf", &din);
1621                 }
1622
1623                 if (ok)
1624                         LR.rel=din;
1625
1626                 fclose(fd);
1627         } 
1628
1629         if (fd==NULL)
1630         {
1631                 /* Create a "splat.lrp" file since one
1632                    could not be successfully loaded. */
1633
1634                 outfile=fopen("splat.lrp","w");
1635
1636                 fprintf(outfile,"%.3f\t; Earth Dielectric Constant (Relative permittivity)\n",LR.eps_dielect);
1637
1638                 fprintf(outfile,"%.3f\t; Earth Conductivity (Siemens per meter)\n", LR.sgm_conductivity);
1639
1640                 fprintf(outfile,"%.3f\t; Atmospheric Bending Constant (N-Units)\n",LR.eno_ns_surfref);
1641
1642                 fprintf(outfile,"%.3f\t; Frequency in MHz (20 MHz to 20 GHz)\n", LR.frq_mhz);
1643
1644                 fprintf(outfile,"%d\t; Radio Climate\n",LR.radio_climate);
1645
1646                 fprintf(outfile,"%d\t; Polarization (0 = Horizontal, 1 = Vertical)\n", LR.pol);
1647
1648                 fprintf(outfile,"%.2f\t; Fraction of situations\n",LR.conf);
1649
1650                 fprintf(outfile, "%.2f\t; Fraction of time\n",LR.rel);
1651
1652                 fprintf(outfile,"\nPlease consult SPLAT! documentation for the meaning and use of this data.\n");
1653
1654                 fclose(outfile);
1655
1656                 fprintf(stderr,"\n%c*** There were problems reading your \"%s\" file! ***\nA \"splat.lrp\" file was written to your directory with default data.\n",7,filename);
1657         }
1658
1659         if (fd==NULL || ok==0)
1660                 fprintf(stderr,"Longley-Rice default parameters have been assumed for this analysis.\n");
1661 }
1662
1663 struct site los(struct site source, struct site destination)
1664 {
1665         /* This function determines whether a line-of-sight path
1666            unobstructed by terrain exists between source (transmitter)
1667            and destination (receiver) based on the geographical
1668            locations of the two sites, their respective antenna
1669            heights above ground, and the terrain between them.
1670            A site structure is returned upon completion.  If the
1671            first character of site.name is ' ', then a clear path
1672            exists between source and destination.  If the first
1673            character is '*', then an obstruction exists, and the
1674            site.lat and site.lon elements of the structure provide
1675            the geographical location of the obstruction. */
1676            
1677         int x;
1678         char block;
1679         struct site test, blockage;
1680         register double distance, tx_alt, rx_alt,
1681                  cos_xmtr_angle, cos_test_angle, test_alt;
1682
1683         ReadPath(source,destination);
1684
1685         distance=5280.0*Distance(source,destination);
1686         tx_alt=earthradius+source.alt+GetElevation(source);
1687         rx_alt=earthradius+destination.alt+GetElevation(destination);
1688
1689         /* Elevation angle of the xmtr (source) as seen by the rcvr */
1690
1691         cos_xmtr_angle=((rx_alt*rx_alt)+(distance*distance)-(tx_alt*tx_alt))/(2.0*rx_alt*distance);
1692
1693         /* Determine the elevation angle of each discrete location
1694            along the path between the receiver and transmitter.
1695
1696            Since obstructions are more likely due to terrain effects
1697            closest to the receiver rather than farther away, we start
1698            looking for potential obstructions from the receiver's
1699            location, and work our way towards the transmitter.
1700            This loop is broken when the first obstruction is
1701            detected.  If we can travel all the way to the transmitter
1702            without detecting an obstruction, then we have a clear
1703            unobstructed path between transmitter and receiver. */
1704
1705         for (x=path.length-1, block=0; x>0 && block==0; x--)
1706         {
1707                 /* Build a structure for each test
1708                    point along the path to be surveyed. */
1709
1710                 test.lat=path.lat[x];
1711                 test.lon=path.lon[x];
1712
1713                 /* Measure the distance between the
1714                    test point and the receiver locations */
1715
1716                 distance=5280.0*Distance(test,destination);
1717                 test_alt=earthradius+path.elevation[x];
1718
1719                 /* Determine the cosine of the elevation of the test
1720                    point as seen from the location of the receiver */
1721
1722                 cos_test_angle=((rx_alt*rx_alt)+(distance*distance)-(test_alt*test_alt))/(2.0*rx_alt*distance);
1723
1724                 /* If the elevation angle to the test point (as seen from
1725                    the receiver) is greater than the elevation angle to the
1726                    transmitter (as seen by the receiver), then we have a
1727                    path obstructed by terrain.  Note: Since we're comparing
1728                    the cosines of these angles rather than the angles
1729                    themselves (eliminating the call to acos() saves
1730                    considerable time), the following "if" statement is
1731                    reversed from what it would normally be if the angles
1732                    were compared. */
1733
1734                 if (cos_xmtr_angle>cos_test_angle)
1735                 {
1736                         block=1;
1737                         blockage.lat=path.lat[x];
1738                         blockage.lon=path.lon[x];
1739                         blockage.alt=path.elevation[x];
1740                         blockage.name[0]='*';
1741                 }
1742         }
1743
1744         if (block==0)
1745         {
1746                 blockage.lat=0.0;
1747                 blockage.lon=0.0;
1748                 blockage.alt=0.0;
1749                 blockage.name[0]=' ';
1750         }
1751
1752         return blockage;
1753 }
1754
1755 void PlotPath(struct site source, struct site destination, char mask_value)
1756 {
1757         /* This function analyzes the path between the source and
1758            destination locations.  It determines which points along
1759            the path have line-of-sight visibility to the source.
1760            Points along with path having line-of-sight visibility
1761            to the source at an AGL altitude equal to that of the
1762            destination location are stored by setting bit 1 in the
1763            mask[][] array, which are displayed in green when PPM
1764            maps are later generated by SPLAT!. */
1765
1766         char block;
1767         int x, y;
1768         register double cos_xmtr_angle, cos_test_angle, test_alt;
1769         double distance, rx_alt, tx_alt;
1770
1771         ReadPath(source,destination);
1772
1773         for (y=0; y<path.length; y++)
1774         {
1775                 /* Test this point only if it hasn't been already
1776                    tested and found to be free of obstructions. */
1777
1778                 if ((GetMask(path.lat[y],path.lon[y])&mask_value)==0)
1779                 {
1780                         distance=5280.0*path.distance[y];
1781                         tx_alt=earthradius+source.alt+path.elevation[0];
1782                         rx_alt=earthradius+destination.alt+path.elevation[y];
1783
1784                         /* Calculate the cosine of the elevation of the
1785                            transmitter as seen at the temp rx point. */
1786
1787                         cos_xmtr_angle=((rx_alt*rx_alt)+(distance*distance)-(tx_alt*tx_alt))/(2.0*rx_alt*distance);
1788
1789                         for (x=y, block=0; x>=0 && block==0; x--)
1790                         {
1791                                 distance=5280.0*(path.distance[y]-path.distance[x]);
1792                                 test_alt=earthradius+path.elevation[x];
1793
1794                                 cos_test_angle=((rx_alt*rx_alt)+(distance*distance)-(test_alt*test_alt))/(2.0*rx_alt*distance);
1795
1796                                 /* Compare these two angles to determine if
1797                                    a blockage exists.  Since we're comparing
1798                                    the cosines of these angles rather than
1799                                    the angles themselves, the following "if"
1800                                    statement is reversed from what it would
1801                                    be if the actual angles were compared. */
1802
1803                                 if (cos_xmtr_angle>cos_test_angle)
1804                                         block=1;
1805                         }
1806
1807                         if (block==0)
1808                                 OrMask(path.lat[y],path.lon[y],mask_value);
1809                 }
1810         }
1811 }
1812
1813 void PlotLRPath(struct site source, struct site destination)
1814 {
1815         /* This function plots the RF signal path loss
1816            between source and destination points based
1817            on the Longley-Rice propagation model. */
1818
1819         char strmode[100];
1820         int x, y, errnum;
1821         double loss;
1822
1823         ReadPath(source,destination);
1824         elev_l[1]=0.04*METERS_PER_MILE;
1825
1826         for (x=0; x<path.length; x++)
1827                 elev_l[x+2]=path.elevation[x]*METERS_PER_FOOT;  
1828          
1829         for (y=0; y<path.length; y++)
1830         {
1831                 /* Test this point only if it
1832                    has not already been tested. */
1833
1834                 if (GetMask(path.lat[y],path.lon[y])==0 && 0.04*y<=max_range)
1835                 {
1836                         elev_l[0]=y+1;
1837
1838                         point_to_point(elev_l,source.alt*METERS_PER_FOOT, 
1839                         destination.alt*METERS_PER_FOOT,
1840                         LR.eps_dielect, LR.sgm_conductivity, LR.eno_ns_surfref,
1841                         LR.frq_mhz, LR.radio_climate, LR.pol, LR.conf, LR.rel,
1842                         loss, strmode, errnum);
1843
1844                         /* Note: PASS BY REFERENCE ... loss and errnum are pass
1845                         by reference, only used in this file by this function */
1846
1847                         if (loss>225.0)
1848                                 loss=225.0;
1849
1850                         if (loss<75.0)
1851                                 loss=75.0;
1852
1853                         loss-=75.0;
1854                         loss/=10.0;
1855                         loss+=1.0;
1856                 
1857                         OrMask(path.lat[y],path.lon[y],((unsigned char)(loss))<<3);
1858                 }
1859
1860                 else if (GetMask(path.lat[y],path.lon[y])==0 && 0.04*y>max_range)
1861                         OrMask(path.lat[y],path.lon[y],1);
1862         }
1863 }
1864
1865 void PlotCoverage(struct site source, double altitude)
1866 {
1867         /* This function performs a 360 degree sweep around the
1868            transmitter site (source location), and plots the
1869            line-of-sight coverage of the transmitter on the SPLAT!
1870            generated topographic map based on a receiver located
1871            at the specified altitude (in feet AGL).  Results
1872            are stored in memory, and written out in the form
1873            of a topographic map when the WritePPM() function
1874            is later invoked. */
1875
1876         double lat, lon, one_pixel;
1877         static unsigned char mask_value;
1878         int z, count;
1879         struct site edge;
1880         unsigned char symbol[4], x;
1881
1882         /* Initialize mask_value */
1883
1884         if (mask_value!=8 && mask_value!=16 && mask_value!=32)
1885                 mask_value=1;
1886
1887         one_pixel=1.0/1200.0;
1888
1889         symbol[0]='.';
1890         symbol[1]='o';
1891         symbol[2]='O';
1892         symbol[3]='o';
1893
1894         count=0;        
1895
1896         fprintf(stdout,"\nComputing line-of-sight coverage of %s with an RX antenna\nat %.2f feet AGL:\n\n 0%c to  25%c ",source.name,altitude,37,37);
1897         fflush(stdout);
1898
1899         /* 18.75=1200 pixels/degree divided by 64 loops
1900            per progress indicator symbol (.oOo) printed. */
1901
1902         z=(int)(18.75*(max_west-min_west));
1903
1904         for (lon=min_west, x=0; lon<=max_west; lon+=one_pixel)
1905         {
1906                 edge.lat=max_north;
1907                 edge.lon=lon;
1908                 edge.alt=altitude;
1909
1910                 PlotPath(source,edge,mask_value);
1911                 count++;
1912
1913                 if (count==z) 
1914                 {
1915                         fprintf(stdout,"%c",symbol[x]);
1916                         fflush(stdout);
1917                         count=0;
1918
1919                         if (x==3)
1920                                 x=0;
1921                         else
1922                                 x++;
1923                 }
1924         }
1925
1926         count=0;
1927         fprintf(stdout,"\n25%c to  50%c ",37,37);
1928         fflush(stdout);
1929         
1930         z=(int)(18.75*(max_north-min_north));
1931
1932         for (lat=max_north, x=0; lat>=min_north; lat-=one_pixel)
1933         {
1934                 edge.lat=lat;
1935                 edge.lon=min_west;
1936                 edge.alt=altitude;
1937
1938                 PlotPath(source,edge,mask_value);
1939                 count++;
1940
1941                 if (count==z) 
1942                 {
1943                         fprintf(stdout,"%c",symbol[x]);
1944                         fflush(stdout);
1945                         count=0;
1946
1947                         if (x==3)
1948                                 x=0;
1949                         else
1950                                 x++;
1951                 }
1952         }
1953
1954         count=0;
1955         fprintf(stdout,"\n50%c to  75%c ",37,37);
1956         fflush(stdout);
1957
1958         z=(int)(18.75*(max_west-min_west));
1959
1960         for (lon=min_west, x=0; lon<=max_west; lon+=one_pixel)
1961         {
1962                 edge.lat=min_north;
1963                 edge.lon=lon;
1964                 edge.alt=altitude;
1965
1966                 PlotPath(source,edge,mask_value);
1967                 count++;
1968
1969                 if (count==z)
1970                 {
1971                         fprintf(stdout,"%c",symbol[x]);
1972                         fflush(stdout);
1973                         count=0;
1974
1975                         if (x==3)
1976                                 x=0;
1977                         else
1978                                 x++;
1979                 }
1980         }
1981
1982         count=0;
1983         fprintf(stdout,"\n75%c to 100%c ",37,37);
1984         fflush(stdout);
1985         
1986         z=(int)(18.75*(max_north-min_north));
1987
1988         for (lat=min_north, x=0; lat<=max_north; lat+=one_pixel)
1989         {
1990                 edge.lat=lat;
1991                 edge.lon=max_west;
1992                 edge.alt=altitude;
1993
1994                 PlotPath(source,edge,mask_value);
1995                 count++;
1996
1997                 if (count==z)
1998                 {
1999                         fprintf(stdout,"%c",symbol[x]);
2000                         fflush(stdout);
2001                         count=0;
2002
2003                         if (x==3)
2004                                 x=0;
2005                         else
2006                                 x++;
2007                 }
2008         }
2009
2010         fprintf(stdout,"\nDone!\n");
2011         fflush(stdout);
2012
2013         /* Assign next mask value */
2014
2015         switch (mask_value)
2016         {
2017                 case 1:
2018                         mask_value=8;
2019                         break;
2020
2021                 case 8:
2022                         mask_value=16;
2023                         break;
2024
2025                 case 16:
2026                         mask_value=32;
2027         }
2028 }
2029
2030 void PlotLRMap(struct site source, double altitude)
2031 {
2032         /* This function performs a 360 degree sweep around the
2033            transmitter site (source location), and plots the
2034            Longley-Rice attenuation on the SPLAT! generated
2035            topographic map based on a receiver located at
2036            the specified altitude (in feet AGL).  Results
2037            are stored in memory, and written out in the form
2038            of a topographic map when the WritePPMLR() function
2039            is later invoked. */
2040
2041         int z, count;
2042         struct site edge;
2043         double lat, lon, one_pixel;
2044         unsigned char symbol[4], x;
2045
2046         one_pixel=1.0/1200.0;
2047
2048         symbol[0]='.';
2049         symbol[1]='o';
2050         symbol[2]='O';
2051         symbol[3]='o';
2052
2053         count=0;
2054
2055         fprintf(stdout,"\nComputing Longley-Rice coverage of %s ", source.name);
2056         fprintf(stdout,"out to a radius\nof %.2f miles with an RX antenna at %.2f feet AGL:\n\n 0%c to  25%c ",max_range,altitude,37,37);
2057         fflush(stdout);
2058
2059         /* 18.75=1200 pixels/degree divided by 64 loops
2060            per progress indicator symbol (.oOo) printed. */
2061
2062         z=(int)(18.75*(max_west-min_west));
2063
2064         for (lon=min_west, x=0; lon<=max_west; lon+=one_pixel)
2065         {
2066                 edge.lat=max_north;
2067                 edge.lon=lon;
2068                 edge.alt=altitude;
2069
2070                 PlotLRPath(source,edge);
2071                 count++;
2072
2073                 if (count==z) 
2074                 {
2075                         fprintf(stdout,"%c",symbol[x]);
2076                         fflush(stdout);
2077                         count=0;
2078
2079                         if (x==3)
2080                                 x=0;
2081                         else
2082                                 x++;
2083                 }
2084         }
2085
2086         count=0;
2087         fprintf(stdout,"\n25%c to  50%c ",37,37);
2088         fflush(stdout);
2089         
2090         z=(int)(18.75*(max_north-min_north));
2091
2092         for (lat=max_north, x=0; lat>=min_north; lat-=one_pixel)
2093         {
2094                 edge.lat=lat;
2095                 edge.lon=min_west;
2096                 edge.alt=altitude;
2097
2098                 PlotLRPath(source,edge);
2099                 count++;
2100
2101                 if (count==z) 
2102                 {
2103                         fprintf(stdout,"%c",symbol[x]);
2104                         fflush(stdout);
2105                         count=0;
2106
2107                         if (x==3)
2108                                 x=0;
2109                         else
2110                                 x++;
2111                 }
2112         }
2113
2114         count=0;
2115         fprintf(stdout,"\n50%c to  75%c ",37,37);
2116         fflush(stdout);
2117
2118         z=(int)(18.75*(max_west-min_west));
2119
2120         for (lon=min_west, x=0; lon<=max_west; lon+=one_pixel)
2121         {
2122                 edge.lat=min_north;
2123                 edge.lon=lon;
2124                 edge.alt=altitude;
2125
2126                 PlotLRPath(source,edge);
2127                 count++;
2128
2129                 if (count==z)
2130                 {
2131                         fprintf(stdout,"%c",symbol[x]);
2132                         fflush(stdout);
2133                         count=0;
2134
2135                         if (x==3)
2136                                 x=0;
2137                         else
2138                                 x++;
2139                 }
2140         }
2141
2142         count=0;
2143         fprintf(stdout,"\n75%c to 100%c ",37,37);
2144         fflush(stdout);
2145         
2146         z=(int)(18.75*(max_north-min_north));
2147
2148         for (lat=min_north, x=0; lat<=max_north; lat+=one_pixel)
2149         {
2150                 edge.lat=lat;
2151                 edge.lon=max_west;
2152                 edge.alt=altitude;
2153
2154                 PlotLRPath(source,edge);
2155                 count++;
2156
2157                 if (count==z)
2158                 {
2159                         fprintf(stdout,"%c",symbol[x]);
2160                         fflush(stdout);
2161                         count=0;
2162
2163                         if (x==3)
2164                                 x=0;
2165                         else
2166                                 x++;
2167                 }
2168         }
2169
2170         fprintf(stdout,"\nDone!\n");
2171         fflush(stdout);
2172 }
2173
2174 void WritePPM(char *filename)
2175 {
2176         /* This function generates a topographic map in Portable Pix Map
2177            (PPM) format based on logarithmically scaled topology data,
2178            as well as the content of flags held in the mask[][] array.
2179            The image created is rotated counter-clockwise 90 degrees
2180            from its representation in dem[][] so that north points
2181            up and east points right in the image generated. */
2182
2183         int indx, x, x0, y0, minlat, minlon;
2184         unsigned width, height, output;
2185         unsigned char found, mask;
2186         char mapfile[255];
2187         double conversion, lat, lon, one_over_gamma, one_pixel;
2188         FILE *fd;
2189
2190         one_pixel=1.0/1200.0;
2191         one_over_gamma=1.0/GAMMA;
2192         conversion=255.0/pow((double)(max_elevation-min_elevation),one_over_gamma);
2193
2194         width=1200*(max_west-min_west);
2195         height=1200*(max_north-min_north);
2196
2197         if (filename[0]==0)
2198                 strncpy(mapfile, "map.ppm\0",8);
2199         else
2200         {
2201                 for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
2202                         mapfile[x]=filename[x];
2203
2204                 mapfile[x]='.';
2205                 mapfile[x+1]='p';
2206                 mapfile[x+2]='p';
2207                 mapfile[x+3]='m';
2208                 mapfile[x+4]=0;
2209         }
2210
2211         fd=fopen(mapfile,"wb");
2212
2213         fprintf(fd,"P6\n%u %u\n255\n",width,height);
2214
2215         fprintf(stdout,"\nWriting \"%s\" (%ux%u pixmap image)... ",mapfile,width,height);
2216         fflush(stdout);
2217
2218         for (lat=(double)max_north; lat>=(double)min_north; lat-=one_pixel)
2219         {
2220                 for (lon=(double)max_west; lon>=(double)min_west; lon-=one_pixel)
2221                 {
2222                         minlat=(int)floor(lat);
2223                         minlon=(int)floor(lon);
2224
2225                         for (indx=0, found=0; indx<MAXSLOTS && found==0;)
2226                                 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
2227                                         found=1;
2228                                 else
2229                                         indx++;
2230                         if (found)
2231                         {
2232                                 x0=(int)((1199.0*(lat-floor(lat)))+0.5);
2233                                 y0=(int)((1199.0*(lon-floor(lon)))+0.5);
2234
2235                                 mask=dem[indx].mask[x0][y0];
2236
2237                                 if (mask&2)
2238                                         /* Text Labels: Red */
2239                                         fprintf(fd,"%c%c%c",255,0,0);
2240
2241                                 else if (mask&4)
2242                                         /* County Boundaries: Light Cyan */
2243                                         fprintf(fd,"%c%c%c",128,128,255);
2244
2245                                 else switch (mask&57)
2246                                 {
2247                                         case 1:
2248                                         /* TX1: Green */
2249                                         fprintf(fd,"%c%c%c",0,255,0);
2250                                         break;
2251
2252                                         case 8:
2253                                         /* TX2: Cyan */
2254                                         fprintf(fd,"%c%c%c",0,255,255);
2255                                         break;
2256
2257                                         case 9:
2258                                         /* TX1 + TX2: Yellow */
2259                                         fprintf(fd,"%c%c%c",255,255,0);
2260                                         break;
2261
2262                                         case 16:
2263                                         /* TX3: Medium Violet */
2264                                         fprintf(fd,"%c%c%c",147,112,219);
2265                                         break;
2266
2267                                         case 17:
2268                                         /* TX1 + TX3: Pink */
2269                                         fprintf(fd,"%c%c%c",255,192,203);
2270                                         break;
2271
2272                                         case 24:
2273                                         /* TX2 + TX3: Orange */
2274                                         fprintf(fd,"%c%c%c",255,165,0);
2275                                         break;
2276
2277                                         case 25:
2278                                         /* TX1 + TX2 + TX3: Dark Green */
2279                                         fprintf(fd,"%c%c%c",0,100,0);
2280                                         break;
2281
2282                                         case 32:
2283                                         /* TX4: Sienna 1 */
2284                                         fprintf(fd,"%c%c%c",255,130,71);
2285                                         break;
2286
2287                                         case 33:
2288                                         /* TX1 + TX4: Green Yellow */
2289                                         fprintf(fd,"%c%c%c",173,255,47);
2290                                         break;
2291
2292                                         case 40:
2293                                         /* TX2 + TX4: Dark Sea Green 1 */
2294                                         fprintf(fd,"%c%c%c",193,255,193);
2295                                         break;
2296
2297                                         case 41:
2298                                         /* TX1 + TX2 + TX4: Blanched Almond */
2299                                         fprintf(fd,"%c%c%c",255,235,205);
2300                                         break;
2301
2302                                         case 48:
2303                                         /* TX3 + TX4: Dark Turquoise */
2304                                         fprintf(fd,"%c%c%c",0,206,209);
2305                                         break;
2306
2307                                         case 49:
2308                                         /* TX1 + TX3 + TX4: Medium Spring Green */
2309                                         fprintf(fd,"%c%c%c",0,250,154);
2310                                         break;
2311
2312                                         case 56:
2313                                         /* TX2 + TX3 + TX4: Tan */
2314                                         fprintf(fd,"%c%c%c",210,180,140);
2315                                         break;
2316
2317                                         case 57:
2318                                         /* TX1 + TX2 + TX3 + TX4: Gold2 */
2319                                         fprintf(fd,"%c%c%c",238,201,0);
2320                                         break;
2321
2322                                         default:
2323                                         /* Water: Medium Blue */
2324                                         if (dem[indx].data[x0][y0]==0)
2325                                                 fprintf(fd,"%c%c%c",0,0,170);
2326                                         else
2327                                         {
2328                                                 /* Elevation: Greyscale */
2329                                                 output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
2330                                                 fprintf(fd,"%c%c%c",output,output,output);
2331                                         }
2332                                 }
2333                         }
2334                 }
2335         }
2336
2337         fclose(fd);
2338         fprintf(stdout,"Done!\n");
2339         fflush(stdout);
2340 }
2341
2342 void WritePPMLR(char *filename)
2343 {
2344         /* This function generates a topographic map in Portable Pix Map
2345            (PPM) format based on the content of flags held in the mask[][] 
2346            array (only).  The image created is rotated counter-clockwise
2347            90 degrees from its representation in dem[][] so that north
2348            points up and east points right in the image generated. */
2349
2350         int indx, x, t, t2, x0, y0, minlat, minlon;
2351         unsigned width, height, output;
2352         unsigned char found, mask;
2353         char mapfile[255];
2354         double conversion, lat, lon, one_over_gamma, one_pixel;
2355         FILE *fd;
2356
2357         one_pixel=1.0/1200.0;
2358         one_over_gamma=1.0/GAMMA;
2359         conversion=255.0/pow((double)(max_elevation-min_elevation),one_over_gamma);
2360
2361         width=1200*(max_west-min_west);
2362         height=1200*(max_north-min_north);
2363
2364         if (filename[0]==0)
2365                 strncpy(mapfile, "map.ppm\0",8);
2366         else
2367         {
2368                 for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
2369                         mapfile[x]=filename[x];
2370
2371                 mapfile[x]='.';
2372                 mapfile[x+1]='p';
2373                 mapfile[x+2]='p';
2374                 mapfile[x+3]='m';
2375                 mapfile[x+4]=0;
2376         }
2377
2378         fd=fopen(mapfile,"wb");
2379
2380         fprintf(fd,"P6\n%u %u\n255\n",width,height+30);
2381
2382         fprintf(stdout,"\nWriting \"%s\" (%ux%u pixmap image)... ",mapfile,width,height+30);
2383         fflush(stdout);
2384
2385         for (lat=(double)max_north; lat>=(double)min_north; lat-=one_pixel)
2386         {
2387                 for (lon=(double)max_west; lon>=(double)min_west; lon-=one_pixel)
2388                 {
2389                         minlat=(int)floor(lat);
2390                         minlon=(int)floor(lon);
2391
2392                         for (indx=0, found=0; indx<MAXSLOTS && found==0;)
2393                                 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
2394                                         found=1;
2395                                 else
2396                                         indx++;
2397                         if (found)
2398                         {
2399                                 x0=(int)((1199.0*(lat-floor(lat)))+0.5);
2400                                 y0=(int)((1199.0*(lon-floor(lon)))+0.5);
2401
2402                                 mask=dem[indx].mask[x0][y0];
2403
2404                                 if (mask&2)
2405                                 {
2406                                         /* Text Labels - Black or Red */
2407                                         if (mask&120)
2408                                                 fprintf(fd,"%c%c%c",0,0,0);
2409                                         else
2410                                                 fprintf(fd,"%c%c%c",255,0,0);
2411                                 }
2412
2413                                 else if (mask&4)
2414                                         /* County Boundaries: Black */
2415                                         fprintf(fd,"%c%c%c",0,0,0);
2416
2417                                 else if (mask&1 && !((mask&248)==192))
2418                                 {
2419                                         /* Outside Analysis Range */
2420                                         /* Display Greyscale / Sea Level */
2421
2422                                         if (dem[indx].data[x0][y0]==0)
2423                                                 fprintf(fd,"%c%c%c",0,0,170);
2424                                         else
2425                                         {
2426                                                 output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
2427                                                 fprintf(fd,"%c%c%c",output,output,output);
2428                                         }
2429                                 }
2430
2431                                 else switch ((mask&248)>>3)
2432                                 {
2433                                         case 0:
2434                                         /* Inside range, but no coverage.
2435                                            Display Sea Level / Terrain */
2436
2437                                         if (dem[indx].data[x0][y0]==0)
2438                                                 fprintf(fd,"%c%c%c",0,0,170);
2439                                         else
2440                                         {
2441                                                 /* Elevation: Greyscale */
2442                                                 output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
2443                                                 fprintf(fd,"%c%c%c",output,output,output);
2444                                         }
2445
2446                                         break;
2447
2448                                         case 1:
2449                                         /* Green */
2450                                         fprintf(fd,"%c%c%c",0,255,0);
2451                                         break;
2452
2453                                         case 2:
2454                                         /* Pink */
2455                                         fprintf(fd,"%c%c%c",255,192,203);
2456                                         break;
2457
2458                                         case 3:
2459                                         /* Cyan */
2460                                         fprintf(fd,"%c%c%c",0,255,255);
2461                                         break;
2462
2463                                         case 4:
2464                                         /* Yellow */
2465                                         fprintf(fd,"%c%c%c",255,255,0);
2466                                         break;
2467
2468                                         case 5:
2469                                         /* Medium Violet */
2470                                         fprintf(fd,"%c%c%c",161,131,224);
2471                                         break;
2472
2473                                         case 6:
2474                                         /* Orange */
2475                                         fprintf(fd,"%c%c%c",255,165,0);
2476                                         break;
2477
2478                                         case 7:
2479                                         /* Light Green */
2480                                         fprintf(fd,"%c%c%c",193,255,193);
2481                                         break;
2482
2483                                         case 8:
2484                                         /* Red Pink */
2485                                         fprintf(fd,"%c%c%c",255,108,108);
2486                                         break;
2487
2488                                         case 9:
2489                                         /* TX1 + TX4: Green Yellow */
2490                                         fprintf(fd,"%c%c%c",173,255,47);
2491                                         break;
2492
2493                                         case 10:
2494                                         /* Blanched Almond */
2495                                         fprintf(fd,"%c%c%c",255,235,184);
2496                                         break;
2497
2498                                         case 11:
2499                                         /* Dark Turquoise */
2500                                         fprintf(fd,"%c%c%c",0,206,209);
2501                                         break;
2502
2503                                         case 12:
2504                                         /* Tan */
2505                                         fprintf(fd,"%c%c%c",210,180,140);
2506                                         break;
2507
2508                                         case 13:
2509                                         /* Magenta 1 */
2510                                         fprintf(fd,"%c%c%c",243,110,205);
2511                                         break;
2512
2513                                         case 14:
2514                                         /* Gold2 */
2515                                         fprintf(fd,"%c%c%c",238,201,0);
2516                                         break;
2517
2518                                         case 15:
2519                                         /* Medium Spring Green */
2520                                         fprintf(fd,"%c%c%c",0,250,154);
2521                                         break;
2522
2523                                         case 16:
2524                                         /* Very light Blue */
2525                                         fprintf(fd,"%c%c%c",244,244,255);
2526                                         break;
2527
2528                                         default:
2529                                         /* Land / Sea */
2530                                         if (dem[indx].data[x0][y0]==0)
2531                                                 fprintf(fd,"%c%c%c",0,0,170);
2532                                         else
2533                                         {
2534                                                 /* Elevation: Greyscale */
2535                                                 output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
2536                                                 fprintf(fd,"%c%c%c",output,output,output);
2537                                         }
2538                                 }
2539                         }
2540                 }
2541         }
2542
2543         x0=width/16;
2544
2545         for (y0=0; y0<30; y0++)
2546         {
2547                 for (indx=0; indx<16; indx++)
2548                 {
2549                         for (x=0; x<x0; x++)
2550                         {
2551                                 t=indx;  
2552                                 t2=indx+8;
2553
2554                                 if (y0>=10 && y0<=18)
2555                                 {  
2556                                         if (t2>9)
2557                                         {
2558                                                 if (x>=11 && x<=17)     
2559                                                         if (smallfont[t2/10][y0-10][x-11])
2560                                                                 t=255; 
2561                                         }
2562
2563                                         if (x>=19 && x<=25)     
2564                                                 if (smallfont[t2%10][y0-10][x-19])
2565                                                         t=255;
2566  
2567                                         if (x>=27 && x<=33)
2568                                                 if (smallfont[0][y0-10][x-27])
2569                                                         t=255; 
2570                                 }
2571
2572                                 switch (t)
2573                                 {
2574                                         case 0:
2575                                         /* Green */
2576                                         fprintf(fd,"%c%c%c",0,255,0);
2577                                         break;
2578
2579                                         case 1:
2580                                         /* Pink */
2581                                         fprintf(fd,"%c%c%c",255,192,203);
2582                                         break;
2583
2584                                         case 2:
2585                                         /* Cyan */
2586                                         fprintf(fd,"%c%c%c",0,255,255);
2587                                         break;
2588
2589                                         case 3:
2590                                         /* Yellow */
2591                                         fprintf(fd,"%c%c%c",255,255,0);
2592                                         break;
2593
2594                                         case 4:
2595                                         /* Medium Violet */
2596                                         fprintf(fd,"%c%c%c",161,131,224);
2597                                         break;
2598
2599                                         case 5:
2600                                         /* Orange */
2601                                         fprintf(fd,"%c%c%c",255,165,0);
2602                                         break;
2603
2604                                         case 6:
2605                                         /* Light Green */
2606                                         fprintf(fd,"%c%c%c",193,255,193);
2607                                         break;
2608
2609                                         case 7:
2610                                         /* Red Pink */
2611                                         fprintf(fd,"%c%c%c",255,108,108);
2612                                         break;
2613
2614                                         case 8:
2615                                         /* Green Yellow */
2616                                         fprintf(fd,"%c%c%c",173,255,47);
2617                                         break;
2618
2619                                         case 9:
2620                                         /* Blanched Almond */
2621                                         fprintf(fd,"%c%c%c",255,235,184);
2622                                         break;
2623
2624                                         case 10:
2625                                         /* Dark Turquoise */
2626                                         fprintf(fd,"%c%c%c",0,206,209);
2627                                         break;
2628
2629                                         case 11:
2630                                         /* Tan */
2631                                         fprintf(fd,"%c%c%c",210,180,140);
2632                                         break;
2633
2634                                         case 12:
2635                                         /* Magenta 1 */
2636                                         fprintf(fd,"%c%c%c",243,110,205);
2637                                         break;
2638
2639                                         case 13:
2640                                         /* Gold2 */
2641                                         fprintf(fd,"%c%c%c",238,201,0);
2642                                         break;
2643
2644                                         case 14:
2645                                         /* Medium Spring Green */
2646                                         fprintf(fd,"%c%c%c",0,250,154);
2647                                         break;
2648
2649                                         case 255:
2650                                         /* Black */
2651                                         fprintf(fd,"%c%c%c",0,0,0);
2652                                         break;
2653
2654                                         default:
2655                                         /* Very Light Blue */
2656                                         fprintf(fd,"%c%c%c",244,244,255);
2657                                 }
2658                         } 
2659                 }
2660         }
2661
2662         fclose(fd);
2663         fprintf(stdout,"Done!\n");
2664         fflush(stdout);
2665 }
2666
2667 void GraphTerrain(struct site source, struct site destination, char *name)
2668 {
2669         /* This function invokes gnuplot to generate an appropriate
2670            output file indicating the terrain profile between the source
2671            and destination locations.  "filename" is the name assigned
2672            to the output file generated by gnuplot.  The filename extension
2673            is used to set gnuplot's terminal setting and output file type.
2674            If no extension is found, .gif is assumed.  */
2675
2676         int x, y, z;
2677         char filename[255], term[15], ext[15];
2678         FILE *fd=NULL;
2679
2680         ReadPath(destination,source);
2681
2682         fd=fopen("profile.gp","wb");
2683
2684         for (x=0; x<path.length; x++)
2685                 fprintf(fd,"%f\t%f\n",path.distance[x],path.elevation[x]);
2686
2687         fclose(fd);
2688
2689         if (name[0]==0)
2690         {
2691                 /* Default filename and output file type */
2692
2693                 strncpy(filename,"profile\0",8);
2694                 strncpy(term,"gif\0",4);
2695                 strncpy(ext,"gif\0",4);
2696         }
2697
2698         else
2699         {
2700                 /* Grab extension and terminal type from "name" */
2701
2702                 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
2703                         filename[x]=name[x];
2704
2705                 if (name[x]=='.')
2706                 {
2707                         for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
2708                         {
2709                                 term[y]=tolower(name[x]);
2710                                 ext[y]=term[y];
2711                         }
2712
2713                         ext[y]=0;
2714                         term[y]=0;
2715                         filename[z]=0;
2716                 }
2717
2718                 else
2719                 {       /* No extension -- Default is gif */
2720
2721                         filename[x]=0;
2722                         strncpy(term,"gif\0",4);
2723                         strncpy(ext,"gif\0",4);
2724                 }
2725         }
2726
2727         /* Either .ps or .postscript may be used
2728            as an extension for postscript output. */
2729
2730         if (strncmp(term,"postscript",10)==0)
2731                 strncpy(ext,"ps\0",3);
2732
2733         else if (strncmp(ext,"ps",2)==0)
2734                 strncpy(term,"postscript\0",11);
2735
2736         fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
2737         fflush(stdout);
2738
2739         fd=fopen("splat.gp","w");
2740         fprintf(fd,"set grid\n");
2741         fprintf(fd,"set autoscale\n");
2742         fprintf(fd,"set term %s\n",term);
2743         fprintf(fd,"set title \"SPLAT! Terrain Profile\"\n");
2744         fprintf(fd,"set xlabel \"Distance Between %s and %s (miles)\"\n",destination.name,source.name);
2745         fprintf(fd,"set ylabel \"Ground Elevation Above Sea Level (feet)\"\n");
2746         fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
2747         fprintf(fd,"plot \"profile.gp\" title \"\" with lines\n");
2748         fclose(fd);
2749                         
2750         x=system("gnuplot splat.gp");
2751
2752         if (x!=-1)
2753         {
2754                 unlink("splat.gp");
2755                 unlink("profile.gp");
2756                 fprintf(stdout," Done!\n");
2757                 fflush(stdout);
2758         }
2759
2760         else
2761                 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
2762 }
2763
2764 void GraphElevation(struct site source, struct site destination, char *name)
2765 {
2766         /* This function invokes gnuplot to generate an appropriate
2767            output file indicating the terrain profile between the source
2768            and destination locations.  "filename" is the name assigned
2769            to the output file generated by gnuplot.  The filename extension
2770            is used to set gnuplot's terminal setting and output file type.
2771            If no extension is found, .gif is assumed. */
2772
2773         int x, y, z;
2774         char filename[255], term[15], ext[15];
2775         double angle, refangle, maxangle=-90.0;
2776         struct site remote;
2777         FILE *fd=NULL, *fd2=NULL;
2778
2779         ReadPath(destination,source);  /* destination=RX, source=TX */
2780         refangle=ElevationAngle(destination,source);
2781
2782         fd=fopen("profile.gp","wb");
2783         fd2=fopen("reference.gp","wb");
2784
2785         for (x=1; x<path.length-1; x++)
2786         {
2787                 remote.lat=path.lat[x];
2788                 remote.lon=path.lon[x];
2789                 remote.alt=0.0;
2790                 angle=ElevationAngle(destination,remote);
2791                 fprintf(fd,"%f\t%f\n",path.distance[x],angle);
2792                 fprintf(fd2,"%f\t%f\n",path.distance[x],refangle);
2793
2794                 if (angle>maxangle)
2795                         maxangle=angle;
2796         }
2797
2798         fprintf(fd,"%f\t%f\n",path.distance[path.length-1],refangle);
2799         fprintf(fd2,"%f\t%f\n",path.distance[path.length-1],refangle);
2800
2801         fclose(fd);
2802         fclose(fd2);
2803
2804         if (name[0]==0)
2805         {
2806                 /* Default filename and output file type */
2807
2808                 strncpy(filename,"profile\0",8);
2809                 strncpy(term,"gif\0",4);
2810                 strncpy(ext,"gif\0",4);
2811         }
2812
2813         else
2814         {
2815                 /* Grab extension and terminal type from "name" */
2816
2817                 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
2818                         filename[x]=name[x];
2819
2820                 if (name[x]=='.')
2821                 {
2822                         for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
2823                         {
2824                                 term[y]=tolower(name[x]);
2825                                 ext[y]=term[y];
2826                         }
2827
2828                         ext[y]=0;
2829                         term[y]=0;
2830                         filename[z]=0;
2831                 }
2832
2833                 else
2834                 {       /* No extension -- Default is gif */
2835
2836                         filename[x]=0;
2837                         strncpy(term,"gif\0",4);
2838                         strncpy(ext,"gif\0",4);
2839                 }
2840         }
2841
2842         /* Either .ps or .postscript may be used
2843            as an extension for postscript output. */
2844
2845         if (strncmp(term,"postscript",10)==0)
2846                 strncpy(ext,"ps\0",3);
2847
2848         else if (strncmp(ext,"ps",2)==0)
2849                 strncpy(term,"postscript\0",11);
2850
2851         fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
2852         fflush(stdout);
2853
2854         fd=fopen("splat.gp","w");
2855
2856         fprintf(fd,"set grid\n");
2857         fprintf(fd,"set yrange [%2.3f to %2.3f]\n", (-fabs(refangle)-0.25), maxangle+0.25);
2858         fprintf(fd,"set term %s\n",term);
2859         fprintf(fd,"set title \"SPLAT! Elevation Profile\"\n");
2860         fprintf(fd,"set xlabel \"Distance Between %s and %s (miles)\"\n",destination.name,source.name);
2861         fprintf(fd,"set ylabel \"Elevation Angle Along Path Between %s and %s (degrees)\"\n",destination.name,source.name);
2862         fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
2863         fprintf(fd,"plot \"profile.gp\" title \"Real Earth Profile\" with lines, \"reference.gp\" title \"Line Of Sight Path\" with lines\n");
2864
2865         fclose(fd);
2866                         
2867         x=system("gnuplot splat.gp");
2868
2869         if (x!=-1)
2870         {
2871                 unlink("splat.gp");
2872                 unlink("profile.gp");
2873                 unlink("reference.gp"); 
2874
2875                 fprintf(stdout," Done!\n");
2876                 fflush(stdout);
2877         }
2878
2879         else
2880                 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
2881 }
2882
2883 void GraphHeight(struct site source, struct site destination, char *name)
2884 {
2885         /* This function invokes gnuplot to generate an appropriate
2886            output file indicating the terrain profile between the source
2887            and destination locations.  What is plotted is the height of land
2888            above or below a straight line between the receibe and transmit
2889            sites. "filename" is the name assigned to the output file
2890            generated by gnuplot.  The filename extension is used to
2891            set gnuplot's terminal setting and output file type.
2892            If no extension is found, .gif is assumed. */
2893
2894         int x, y, z;
2895         char filename[255], term[15], ext[15];
2896         double a, b, c, height, refangle, cangle, maxheight=-100000.0,
2897                minheight=100000.0;
2898         struct site remote;
2899         FILE *fd=NULL, *fd2=NULL;
2900
2901         ReadPath(destination,source);  /* destination=RX, source=TX */
2902         refangle=ElevationAngle(destination,source);
2903         b=GetElevation(destination)+destination.alt+earthradius;
2904
2905         fd=fopen("profile.gp","wb");
2906         fd2=fopen("reference.gp","wb");
2907
2908         for (x=1; x<path.length-1; x++)
2909         {
2910                 remote.lat=path.lat[x];
2911                 remote.lon=path.lon[x];
2912                 remote.alt=0.0;
2913
2914                 a=GetElevation(remote)+earthradius;
2915
2916                 cangle=5280.0*Distance(destination,remote)/earthradius;
2917
2918                 c=b*sin(refangle*deg2rad+HALFPI)/sin(HALFPI-refangle*deg2rad-cangle);
2919
2920                 height=a-c;
2921
2922                 fprintf(fd,"%f\t%f\n",path.distance[x],height);
2923                 fprintf(fd2,"%f\t%f\n",path.distance[x],0.);
2924
2925                 if (height>maxheight)
2926                         maxheight=height;
2927
2928                 if (height<minheight)
2929                         minheight=height;
2930         }
2931
2932         fprintf(fd,"%f\t%f\n",path.distance[path.length-1],0.0);
2933         fprintf(fd2,"%f\t%f\n",path.distance[path.length-1],0.0);
2934
2935         fclose(fd);
2936         fclose(fd2);
2937
2938         if (name[0]==0)
2939         {
2940                 /* Default filename and output file type */
2941
2942                 strncpy(filename,"height\0",8);
2943                 strncpy(term,"gif\0",4);
2944                 strncpy(ext,"gif\0",4);
2945         }
2946
2947         else
2948         {
2949                 /* Grab extension and terminal type from "name" */
2950
2951                 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
2952                         filename[x]=name[x];
2953
2954                 if (name[x]=='.')
2955                 {
2956                         for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
2957                         {
2958                                 term[y]=tolower(name[x]);
2959                                 ext[y]=term[y];
2960                         }
2961
2962                         ext[y]=0;
2963                         term[y]=0;
2964                         filename[z]=0;
2965                 }
2966
2967                 else
2968                 {       /* No extension -- Default is gif */
2969
2970                         filename[x]=0;
2971                         strncpy(term,"gif\0",4);
2972                         strncpy(ext,"gif\0",4);
2973                 }
2974         }
2975
2976         /* Either .ps or .postscript may be used
2977            as an extension for postscript output. */
2978
2979         if (strncmp(term,"postscript",10)==0)
2980                 strncpy(ext,"ps\0",3);
2981
2982         else if (strncmp(ext,"ps",2)==0)
2983                 strncpy(term,"postscript\0",11);
2984
2985         fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
2986         fflush(stdout);
2987
2988         fd=fopen("splat.gp","w");
2989
2990         minheight-=20.0;
2991         maxheight+=20.0;
2992
2993         if (maxheight<20.0)
2994                 maxheight=20.0;
2995
2996         fprintf(fd,"set grid\n");
2997         fprintf(fd,"set yrange [%2.3f to %2.3f]\n", minheight, maxheight);
2998         fprintf(fd,"set term %s\n",term);
2999         fprintf(fd,"set title \"SPLAT! Height Profile\"\n");
3000         fprintf(fd,"set xlabel \"Distance Between %s and %s (miles)\"\n",destination.name,source.name);
3001         fprintf(fd,"set ylabel \"Ground Height Above Path Between %s and %s (feet)\"\n",destination.name,source.name);
3002         fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
3003         fprintf(fd,"plot \"profile.gp\" title \"Real Earth Profile\" with lines, \"reference.gp\" title \"Line Of Sight Path\" with lines\n");
3004
3005         fclose(fd);
3006                         
3007         x=system("gnuplot splat.gp");
3008
3009         if (x!=-1)
3010         {
3011                 unlink("splat.gp");
3012                 unlink("profile.gp");
3013                 unlink("reference.gp"); 
3014                 fprintf(stdout," Done!\n");
3015                 fflush(stdout);
3016         }
3017
3018         else
3019                 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
3020 }
3021
3022 void GraphLongley(struct site source, struct site destination, char *name)
3023 {
3024         /* This function invokes gnuplot to generate an appropriate
3025            output file indicating the Longley-Rice model loss between
3026            the source and destination locations.   "filename" is the
3027            name assigned to the output file generated by gnuplot.
3028            The filename extension is used to set gnuplot's terminal
3029            setting and output file type.  If no extension is found,
3030            .gif is assumed.  */
3031
3032         int x, y, z, errnum, errflag=0;
3033         char filename[255], term[15], ext[15], strmode[100], report_name[80];
3034         double maxloss=-100000.0, minloss=100000.0, loss, haavt, angle;
3035         FILE *fd=NULL, *fd2=NULL;
3036
3037         sprintf(report_name,"%s-to-%s.lro",source.name,destination.name);
3038
3039         for (x=0; report_name[x]!=0; x++)
3040                 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
3041                         report_name[x]='_';     
3042
3043         fd2=fopen(report_name,"w");
3044
3045         fprintf(fd2,"\n\t--==[ SPLAT! v%s Longley-Rice Model Path Loss Report ]==--\n\n",splat_version);
3046         fprintf(fd2,"Analysis of RF path conditions between %s and %s:\n",source.name, destination.name);
3047         fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
3048         fprintf(fd2,"Transmitter site: %s\n",source.name);
3049         fprintf(fd2,"Site location: %.4f North / %.4f West",source.lat, source.lon);
3050         fprintf(fd2, " (%s N / ", dec2dms(source.lat));
3051         fprintf(fd2, "%s W)\n", dec2dms(source.lon));
3052         fprintf(fd2,"Ground elevation: %.2f feet AMSL\n",GetElevation(source));
3053         fprintf(fd2,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",source.alt, source.alt+GetElevation(source));
3054
3055         haavt=haat(source);
3056
3057         if (haavt>-4999.0)
3058                 fprintf(fd2,"Antenna height above average terrain: %.2f feet\n",haavt);
3059
3060         fprintf(fd2,"Distance to %s: %.2f miles.\n",destination.name,Distance(source,destination));
3061         fprintf(fd2,"Azimuth to %s: %.2f degrees.\n",destination.name,Azimuth(source,destination));
3062
3063         angle=ElevationAngle(source,destination);
3064
3065         if (angle>=0.0)
3066                 fprintf(fd2,"Angle of elevation between %s and %s: %+.4f degrees.\n",source.name,destination.name,angle);
3067
3068         if (angle<0.0)
3069                 fprintf(fd2,"Angle of depression between %s and %s: %+.4f degrees.\n",source.name,destination.name,angle);
3070
3071         fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
3072
3073         /* Receiver */
3074
3075         fprintf(fd2,"Receiver site: %s\n",destination.name);
3076         fprintf(fd2,"Site location: %.4f North / %.4f West",destination.lat, destination.lon);
3077         fprintf(fd2, " (%s N / ", dec2dms(destination.lat));
3078         fprintf(fd2, "%s W)\n", dec2dms(destination.lon));
3079         fprintf(fd2,"Ground elevation: %.2f feet AMSL\n",GetElevation(destination));
3080         fprintf(fd2,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",destination.alt, destination.alt+GetElevation(destination));
3081
3082         haavt=haat(destination);
3083
3084         if (haavt>-4999.0)
3085                 fprintf(fd2,"Antenna height above average terrain: %.2f feet\n",haavt);
3086
3087         fprintf(fd2,"Distance to %s: %.2f miles.\n",source.name,Distance(source,destination));
3088         fprintf(fd2,"Azimuth to %s: %.2f degrees.\n",source.name,Azimuth(destination,source));
3089
3090         angle=ElevationAngle(destination,source);
3091
3092         if (angle>=0.0)
3093                 fprintf(fd2,"Angle of elevation between %s and %s: %+.4f degrees.\n",destination.name,source.name,angle);
3094
3095         if (angle<0.0)
3096                 fprintf(fd2,"Angle of depression between %s and %s: %+.4f degrees.\n",destination.name,source.name,angle);
3097
3098         fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
3099
3100         fprintf(fd2,"Longley-Rice path calculation parameters used in this analysis:\n\n");
3101         fprintf(fd2,"Earth's Dielectric Constant: %.3lf\n",LR.eps_dielect);
3102         fprintf(fd2,"Earth's Conductivity: %.3lf\n",LR.sgm_conductivity);
3103         fprintf(fd2,"Atmospheric Bending Constant (N): %.3lf\n",LR.eno_ns_surfref);
3104         fprintf(fd2,"Frequency: %.3lf (MHz)\n",LR.frq_mhz);
3105         fprintf(fd2,"Radio Climate: %d (",LR.radio_climate);
3106
3107         switch (LR.radio_climate)
3108         {
3109                 case 1:
3110                 fprintf(fd2,"Equatorial");
3111                 break;
3112
3113                 case 2:
3114                 fprintf(fd2,"Continental Subtropical");
3115                 break;
3116
3117                 case 3:
3118                 fprintf(fd2,"Maritime Subtropical");
3119                 break;
3120
3121                 case 4:
3122                 fprintf(fd2,"Desert");
3123                 break;
3124
3125                 case 5:
3126                 fprintf(fd2,"Continental Temperate");
3127                 break;
3128
3129                 case 6:
3130                 fprintf(fd2,"Martitime Temperate, Over Land");
3131                 break;
3132
3133                 case 7:
3134                 fprintf(fd2,"Maritime Temperate, Over Sea");
3135                 break;
3136
3137                 default:
3138                 fprintf(fd2,"Unknown");
3139         }
3140
3141         fprintf(fd2,")\nPolarization: %d (",LR.pol);
3142
3143         if (LR.pol==0)
3144                 fprintf(fd2,"Horizontal");
3145
3146         if (LR.pol==1)
3147                 fprintf(fd2,"Vertical");
3148
3149         fprintf(fd2,")\nFraction of Situations: %.1lf%c\n",LR.conf*100.0,37);
3150         fprintf(fd2,"Fraction of Time: %.1lf%c\n",LR.rel*100.0,37);
3151
3152         fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
3153
3154         fprintf(fd2,"Analysis Results:\n\n");
3155
3156         ReadPath(source, destination);  /* destination=RX, source=TX */
3157
3158         elev_l[1]=0.04*METERS_PER_MILE;
3159
3160         for (x=1; x<path.length; x++)
3161                 elev_l[x+1]=path.elevation[x]*METERS_PER_FOOT;  
3162
3163         fprintf(fd2,"Distance (mi)\tLoss (dB)\tErrnum\tComment\n\n"); 
3164
3165         fd=fopen("profile.gp","w");
3166
3167         for (x=2; x<path.length; x++)
3168         {
3169                 elev_l[0]=x-1;
3170
3171                 point_to_point(elev_l, source.alt*METERS_PER_FOOT, 
3172                         destination.alt*METERS_PER_FOOT,
3173                         LR.eps_dielect, LR.sgm_conductivity, LR.eno_ns_surfref,
3174                         LR.frq_mhz, LR.radio_climate, LR.pol, LR.conf, LR.rel,
3175                         loss, strmode, errnum);
3176
3177                         /* Note: PASS BY REFERENCE ... loss and errnum are pass
3178                         by reference, only used in this file by this function */
3179
3180
3181                 fprintf(fd,"%f\t%f\n",path.distance[path.length-1]-path.distance[x],loss);
3182                 fprintf(fd2,"%7.2f\t\t%7.2f\t\t  %d\t%s\n",path.distance[x],loss, errnum, strmode); 
3183                 errflag|=errnum;
3184                   
3185                 if (loss>maxloss)
3186                         maxloss=loss;
3187
3188                 if (loss<minloss)
3189                         minloss=loss;
3190         }
3191
3192         fclose(fd);
3193
3194         if (errflag)
3195         {
3196                 fprintf(fd2,"\nNotes on \"errnum\"...\n\n");
3197                 fprintf(fd2,"  0: No error.  :-)\n");
3198                 fprintf(fd2,"  1: Warning!  Some parameters are nearly out of range.\n");
3199                 fprintf(fd2,"     Results should be used with caution.\n");
3200                 fprintf(fd2,"  2: Note: Default parameters have been substituted for impossible ones.\n");
3201                 fprintf(fd2,"  3: Warning!  A combination of parameters is out of range.\n");
3202                 fprintf(fd2,"     Results are probably invalid.\n");
3203                 fprintf(fd2,"  Other: Warning!  Some parameters are out of range.\n");
3204                 fprintf(fd2,"    Results are probably invalid.\n\nEnd of Report\n");
3205         }
3206
3207         fprintf(stdout,"Longley-Rice Path Loss between %s and %s is %.2f db\n",source.name, destination.name, loss);
3208         fprintf(stdout,"Path Loss Report written to: \"%s\"\n",report_name);
3209         fflush(stdout);
3210
3211         fclose(fd2);
3212
3213         if (name[0]==0)
3214         {
3215                 /* Default filename and output file type */
3216
3217                 strncpy(filename,"loss\0",5);
3218                 strncpy(term,"gif\0",4);
3219                 strncpy(ext,"gif\0",4);
3220         }
3221
3222         else
3223         {
3224                 /* Grab extension and terminal type from "name" */
3225
3226                 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
3227                         filename[x]=name[x];
3228
3229                 if (name[x]=='.')
3230                 {
3231                         for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
3232                         {
3233                                 term[y]=tolower(name[x]);
3234                                 ext[y]=term[y];
3235                         }
3236
3237                         ext[y]=0;
3238                         term[y]=0;
3239                         filename[z]=0;
3240                 }
3241
3242                 else
3243                 {       /* No extension -- Default is gif */
3244
3245                         filename[x]=0;
3246                         strncpy(term,"gif\0",4);
3247                         strncpy(ext,"gif\0",4);
3248                 }
3249         }
3250
3251         /* Either .ps or .postscript may be used
3252            as an extension for postscript output. */
3253
3254         if (strncmp(term,"postscript",10)==0)
3255                 strncpy(ext,"ps\0",3);
3256
3257         else if (strncmp(ext,"ps",2)==0)
3258                 strncpy(term,"postscript\0",11);
3259
3260         fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
3261         fflush(stdout);
3262
3263         fd=fopen("splat.gp","w");
3264
3265         fprintf(fd,"set grid\n");
3266         fprintf(fd,"set yrange [%2.3f to %2.3f]\n", minloss, maxloss);
3267         fprintf(fd,"set term %s\n",term);
3268         fprintf(fd,"set title \"SPLAT! Loss Profile\"\n");
3269         fprintf(fd,"set xlabel \"Distance Between %s and %s (miles)\"\n",destination.name,source.name);
3270         fprintf(fd,"set ylabel \"Longley-Rice Loss (dB)\"\n");
3271         fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
3272         fprintf(fd,"plot \"profile.gp\" title \"Longley-Rice Loss\" with lines\n");
3273
3274         fclose(fd);
3275                         
3276         x=system("gnuplot splat.gp");
3277
3278         if (x!=-1)
3279         {
3280                 unlink("splat.gp");
3281                 unlink("profile.gp");
3282                 unlink("reference.gp"); 
3283
3284                 fprintf(stdout," Done!\n");
3285                 fflush(stdout);
3286         }
3287
3288         else
3289                 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
3290 }
3291
3292 void ObstructionReport(struct site xmtr, struct site rcvr, char report)
3293 {
3294         struct site result, result2, new_site;
3295         double angle, haavt;
3296         unsigned char block;
3297         char report_name[80], string[255];
3298         int x;
3299         FILE *fd;
3300
3301         sprintf(report_name,"%s-to-%s.txt",xmtr.name,rcvr.name);
3302
3303         for (x=0; report_name[x]!=0; x++)
3304                 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
3305                         report_name[x]='_';     
3306
3307         fd=fopen(report_name,"w");
3308
3309         fprintf(fd,"\n\t\t--==[ SPLAT! v%s Obstruction Report ]==--\n\n",splat_version);
3310         fprintf(fd,"Analysis of line-of-sight path conditions between %s and %s:\n",xmtr.name, rcvr.name);
3311         fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
3312         fprintf(fd,"Transmitter site: %s\n",xmtr.name);
3313         fprintf(fd,"Site location: %.4f North / %.4f West",xmtr.lat, xmtr.lon);
3314         fprintf(fd, " (%s N / ", dec2dms(xmtr.lat));
3315         fprintf(fd, "%s W)\n", dec2dms(xmtr.lon));
3316         fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(xmtr));
3317         fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",xmtr.alt, xmtr.alt+GetElevation(xmtr));
3318
3319         haavt=haat(xmtr);
3320
3321         if (haavt>-4999.0)
3322                 fprintf(fd,"Antenna height above average terrain: %.2f feet\n",haavt);
3323
3324         fprintf(fd,"Distance to %s: %.2f miles.\n",rcvr.name,Distance(xmtr,rcvr));
3325         fprintf(fd,"Azimuth to %s: %.2f degrees.\n",rcvr.name,Azimuth(xmtr,rcvr));
3326
3327         angle=ElevationAngle(xmtr,rcvr);
3328
3329         if (angle>=0.0)
3330                 fprintf(fd,"Angle of elevation between %s and %s: %+.4f degrees.\n",xmtr.name,rcvr.name,angle);
3331
3332         if (angle<0.0)
3333                 fprintf(fd,"Angle of depression between %s and %s: %+.4f degrees.\n",xmtr.name,rcvr.name,angle);
3334
3335         fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
3336
3337         /* Receiver */
3338
3339         fprintf(fd,"Receiver site: %s\n",rcvr.name);
3340         fprintf(fd,"Site location: %.4f North / %.4f West",rcvr.lat, rcvr.lon);
3341         fprintf(fd, " (%s N / ", dec2dms(rcvr.lat));
3342         fprintf(fd, "%s W)\n", dec2dms(rcvr.lon));
3343         fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(rcvr));
3344         fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",rcvr.alt, rcvr.alt+GetElevation(rcvr));
3345
3346         haavt=haat(rcvr);
3347
3348         if (haavt>-4999.0)
3349                 fprintf(fd,"Antenna height above average terrain: %.2f feet\n",haavt);
3350
3351         fprintf(fd,"Distance to %s: %.2f miles.\n",xmtr.name,Distance(xmtr,rcvr));
3352         fprintf(fd,"Azimuth to %s: %.2f degrees.\n",xmtr.name,Azimuth(rcvr,xmtr));
3353
3354         angle=ElevationAngle(rcvr,xmtr);
3355
3356         if (angle>=0.0)
3357                 fprintf(fd,"Angle of elevation between %s and %s: %+.4f degrees.\n",rcvr.name,xmtr.name,angle);
3358
3359         if (angle<0.0)
3360                 fprintf(fd,"Angle of depression between %s and %s: %+.4f degrees.\n",rcvr.name,xmtr.name,angle);
3361
3362         fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
3363
3364         if (report=='y')
3365         {
3366                 /* Write an Obstruction Report */
3367
3368                 new_site=rcvr;
3369                 result=los(xmtr,rcvr);
3370                 result2=result;
3371                 result2.alt-=1;
3372                 block=result.name[0];
3373
3374                 if (block=='*')
3375                         fprintf(fd,"SPLAT! detected obstructions at:\n\n");
3376
3377                 while (block=='*')
3378                 {
3379                         if (result.lat!=result2.lat || result.lon!=result2.lon || result.alt!=result2.alt)
3380                                 fprintf(fd,"\t%.4f N, %.4f W, %5.2f miles, %6.2f feet AMSL.\n",result.lat, result.lon, Distance(rcvr,result), result.alt);
3381
3382                         result2=result;
3383                         new_site.alt+=1.0;
3384
3385                         /* Can you hear me now? :-) */
3386
3387                         result=los(xmtr,new_site);
3388                         block=result.name[0];
3389                 }
3390
3391                 if (new_site.alt!=rcvr.alt)
3392                         sprintf(string,"\nAntenna at %s must be raised to at least %.2f feet AGL\nto clear all obstructions detected by SPLAT!\n\n",rcvr.name, new_site.alt);
3393                 else
3394                         sprintf(string,"\nNo obstructions due to terrain were detected by SPLAT!\n\n");
3395         }
3396
3397         fprintf(fd,"%s",string);
3398
3399         fclose(fd);
3400
3401         /* Display LOS status to terminal */
3402
3403         fprintf(stdout,"%sObstruction report written to: \"%s\"\n",string,report_name);
3404         fflush(stdout);
3405 }
3406
3407 void SiteReport(struct site xmtr)
3408 {
3409         char report_name[80];
3410         double terrain;
3411         int x, azi;
3412         FILE *fd;
3413
3414         sprintf(report_name,"%s-site_report.txt",xmtr.name);
3415
3416         for (x=0; report_name[x]!=0; x++)
3417                 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
3418                         report_name[x]='_';     
3419
3420         fd=fopen(report_name,"w");
3421
3422         fprintf(fd,"\n\t--==[ SPLAT! v%s Site Analysis Report For: %s ]==--\n\n",splat_version,xmtr.name);
3423
3424         fprintf(fd,"---------------------------------------------------------------------------\n\n");
3425         fprintf(fd,"Site location: %.4f North / %.4f West",xmtr.lat, xmtr.lon);
3426         fprintf(fd, " (%s N / ",dec2dms(xmtr.lat));
3427         fprintf(fd, "%s W)\n",dec2dms(xmtr.lon));
3428         fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(xmtr));
3429         fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",xmtr.alt, xmtr.alt+GetElevation(xmtr));
3430
3431         terrain=haat(xmtr);
3432
3433         if (terrain>-4999.0)
3434         {
3435                 fprintf(fd,"Antenna height above average terrain: %.2f feet\n\n",terrain);
3436
3437                 /* Display the average terrain between 2 and 10 miles
3438                    from the transmitter site at azimuths of 0, 45, 90,
3439                    135, 180, 225, 270, and 315 degrees. */
3440
3441                 for (azi=0; azi<=315; azi+=45)
3442                 {
3443                         fprintf(fd,"Average terrain at %3d degrees azimuth: ",azi);
3444                         terrain=AverageTerrain(xmtr,(double)azi,2.0,10.0);
3445
3446                         if (terrain>-4999.0)
3447                                 fprintf(fd,"%.2f feet AMSL\n",terrain);
3448                         else
3449                                 fprintf(fd,"No terrain\n");
3450                 }
3451         }
3452
3453         fprintf(fd,"\n---------------------------------------------------------------------------\n\n");
3454         fclose(fd);
3455         fprintf(stdout,"\nSite analysis report written to: \"%s\"\n",report_name);
3456 }
3457
3458 int main(char argc, char *argv[])
3459 {
3460
3461         int             x, y, z=0;
3462         unsigned char   rxlat, rxlon, txlat, txlon, min_lat,
3463                         min_lon, max_lat, max_lon, 
3464                         coverage=0, LRmap=0,
3465                         ext[20], terrain_plot=0,
3466                         elevation_plot=0, height_plot=0, 
3467                         longley_plot=0, cities=0, bfs=0, txsites=0,
3468                         count, west_min, west_max, north_min, north_max,
3469                         report='y';
3470  
3471         char            mapfile[255], header[80], city_file[5][255], 
3472                         elevation_file[255], height_file[255], 
3473                         longley_file[255], terrain_file[255],
3474                         string[255], rxfile[255],
3475                         txfile[255], map=0, boundary_file[5][255];
3476
3477         double          altitude=0.0, altitudeLR=0.0, tx_range=0.0,
3478                         rx_range=0.0, deg_range=0.0, deg_limit,
3479                         deg_range_lon, er_mult;
3480         struct          site tx_site[4], rx_site;
3481         FILE            *fd;
3482
3483         sprintf(header,"\n  --==[ SPLAT! v%s Terrain Analysis Software (c) 1997-2004 KD2BD ]==--\n\n", splat_version);
3484
3485         if (argc==1)
3486         {
3487                 fprintf(stdout, "%sAvailable Options...\n\n\t-t txsite(s).qth (max of 4)\n\t-r rxsite.qth\n",header);
3488                 fprintf(stdout,"\t-c plot coverage area(s) of TX(s) based on an RX antenna at X feet AGL\n");
3489                 fprintf(stdout,"\t-L plot path loss map of TX based on an RX antenna at X feet AGL\n");
3490                 fprintf(stdout,"\t-s filename(s) of city/site file(s) to import (max of 5)\n");
3491                 fprintf(stdout,"\t-b filename(s) of cartographic boundary file(s) to import (max of 5)\n");
3492                 fprintf(stdout,"\t-p filename of terrain profile graph to plot\n");
3493                 fprintf(stdout,"\t-e filename of terrain elevation graph to plot\n");
3494                 fprintf(stdout,"\t-h filename of terrain height graph to plot\n");
3495                 fprintf(stdout,"\t-l filename of Longley-Rice graph to plot\n");
3496                 fprintf(stdout,"\t-o filename of topographic map to generate (.ppm)\n");
3497                 fprintf(stdout,"\t-d sdf file directory path (overrides path in ~/.splat_path file)\n");
3498                 fprintf(stdout,"\t-n no analysis, brief report\n\t-N no analysis, no report\n");
3499                 fprintf(stdout,"\t-m earth radius multiplier\n");
3500                 fprintf(stdout,"\t-R modify default range for -c or -L (miles)\n\n");
3501
3502                 fprintf(stdout,"Type 'man splat', or see the documentation for more details.\n\n");
3503                 fflush(stdout);
3504                 return 1;
3505         }
3506
3507         y=argc-1;
3508
3509         rxfile[0]=0;
3510         txfile[0]=0;
3511         string[0]=0;
3512         mapfile[0]=0;
3513         elevation_file[0]=0;
3514         terrain_file[0]=0;
3515         sdf_path[0]=0;
3516         path.length=0;
3517         rx_site.lat=0.0;
3518         rx_site.lon=0.0;
3519         earthradius=EARTHRADIUS;
3520
3521         for (x=0; x<MAXSLOTS; x++)
3522         {
3523                 dem[x].min_el=0;
3524                 dem[x].max_el=0;
3525                 dem[x].min_north=0;
3526                 dem[x].max_north=0;
3527                 dem[x].min_west=0;
3528                 dem[x].max_west=0;
3529         }
3530
3531         /* Scan for command line arguments */
3532
3533         for (x=1; x<=y; x++)
3534         {
3535                 if (strcmp(argv[x],"-R")==0)
3536                 {
3537                         z=x+1;
3538
3539                         if (z<=y && argv[z][0] && argv[z][0]!='-')
3540                         {
3541                                 sscanf(argv[z],"%lf",&max_range);
3542
3543                                 if (max_range<0.0)
3544                                         max_range=0.0;
3545
3546                                 if (max_range>1000.0)
3547                                         max_range=1000.0;
3548                         }                        
3549                 }
3550
3551                 if (strcmp(argv[x],"-m")==0)
3552                 {
3553                         z=x+1;
3554
3555                         if (z<=y && argv[z][0] && argv[z][0]!='-')
3556                         {
3557                                 sscanf(argv[z],"%lf",&er_mult);
3558
3559                                 if (er_mult<0.1)
3560                                         er_mult=1.0;
3561
3562                                 if (er_mult>1.0e6)
3563                                         er_mult=1.0e6;
3564
3565                                 earthradius*=er_mult;
3566                         }                        
3567                 }
3568
3569                 if (strcmp(argv[x],"-o")==0)
3570                 {
3571                         z=x+1;
3572
3573                         if (z<=y && argv[z][0] && argv[z][0]!='-')
3574                                 strncpy(mapfile,argv[z],253);
3575                         map=1;
3576                 }
3577
3578                 if (strcmp(argv[x],"-c")==0)
3579                 {
3580                         z=x+1;
3581
3582                         if (z<=y && argv[z][0] && argv[z][0]!='-')
3583                         {
3584                                 sscanf(argv[z],"%lf",&altitude);
3585                                 coverage=1;
3586                         }
3587                 }
3588
3589                 if (strcmp(argv[x],"-p")==0)
3590                 { 
3591                         z=x+1;
3592
3593                         if (z<=y && argv[z][0] && argv[z][0]!='-')
3594                         {
3595                                 strncpy(terrain_file,argv[z],253);
3596                                 terrain_plot=1;
3597                         }
3598                 }
3599
3600                 if (strcmp(argv[x],"-e")==0)
3601                 {
3602                         z=x+1;
3603
3604                         if (z<=y && argv[z][0] && argv[z][0]!='-')
3605                         {
3606                                 strncpy(elevation_file,argv[z],253);
3607                                 elevation_plot=1;
3608                         }
3609                 }
3610
3611                 if (strcmp(argv[x],"-h")==0)
3612                 {
3613                         z=x+1;
3614
3615                         if (z<=y && argv[z][0] && argv[z][0]!='-')
3616                         {
3617                                 strncpy(height_file,argv[z],253);
3618                                 height_plot=1;
3619                         }
3620                 }
3621
3622                 if (strcmp(argv[x],"-n")==0)
3623                 {
3624                         if (z<=y && argv[z][0] && argv[z][0]!='-')
3625                         {
3626                                 report='n';
3627                                 map=1;
3628                         }
3629                 }
3630
3631                 if (strcmp(argv[x],"-N")==0)
3632                 {
3633                         if (z<=y && argv[z][0] && argv[z][0]!='-');
3634                         {
3635                                 report='N';
3636                                 map=1;
3637                         }
3638                 }
3639
3640                 if (strcmp(argv[x],"-d")==0)
3641                 {
3642                         z=x+1;
3643
3644                         if (z<=y && argv[z][0] && argv[z][0]!='-')
3645                                 strncpy(sdf_path,argv[z],253);
3646                 }
3647
3648                 if (strcmp(argv[x],"-t")==0)
3649                 {
3650                         /* Read Transmitter Location */
3651
3652                         z=x+1;
3653
3654                         while (z<=y && argv[z][0] && argv[z][0]!='-' && txsites<4)
3655                         {
3656                                 strncpy(txfile,argv[z],253);
3657                                 tx_site[txsites]=LoadQTH(txfile);
3658                                 txsites++;
3659                                 z++;
3660                         }
3661                         z--;
3662                 }
3663
3664                 if (strcmp(argv[x],"-L")==0)
3665                 {
3666                         z=x+1;
3667
3668                         if (z<=y && argv[z][0] && argv[z][0]!='-')
3669                         {
3670                                 sscanf(argv[z],"%lf",&altitudeLR);
3671
3672                                 if (coverage)
3673                                         fprintf(stdout,"c and L are exclusive options, ignoring L.\n");
3674                                 else
3675                                 {
3676                                         LRmap=1;
3677                                         ReadLRParm(txfile);
3678                                 } 
3679                         }
3680                 }
3681
3682                 if (strcmp(argv[x],"-l")==0)
3683                 {
3684                         z=x+1;
3685
3686                         if (z<=y && argv[z][0] && argv[z][0]!='-')
3687                         {
3688                                 strncpy(longley_file,argv[z],253);
3689                                 longley_plot=1;
3690                                 /* Doing this twice is harmless */
3691                                 ReadLRParm(txfile);
3692                         }
3693                 }
3694
3695                 if (strcmp(argv[x],"-r")==0)
3696                 {
3697                         /* Read Receiver Location */
3698
3699                         z=x+1;
3700
3701                         if (z<=y && argv[z][0] && argv[z][0]!='-')
3702                         {
3703                                 strncpy(rxfile,argv[z],253);
3704                                 rx_site=LoadQTH(rxfile);
3705                         }
3706                 }
3707
3708                 if (strcmp(argv[x],"-s")==0)
3709                 {
3710                         /* Read city file(s) */
3711
3712                         z=x+1;
3713
3714                         while (z<=y && argv[z][0] && argv[z][0]!='-' && cities<5)
3715                         {
3716                                 strncpy(city_file[cities],argv[z],253);
3717                                 cities++;
3718                                 z++;
3719                         }
3720                         z--;
3721                 }
3722
3723                 if (strcmp(argv[x],"-b")==0)
3724                 {
3725                         /* Read Boundary File(s) */
3726
3727                         z=x+1;
3728
3729                         while (z<=y && argv[z][0] && argv[z][0]!='-' && bfs<5)
3730                         {
3731                                 strncpy(boundary_file[bfs],argv[z],253);
3732                                 bfs++;
3733                                 z++;
3734                         }
3735                         z--;
3736                 }
3737         }
3738
3739         /* Perform some error checking on the arguments
3740            and switches parsed from the command-line.
3741            If an error is encountered, print a message
3742            and exit gracefully. */
3743
3744         if (txsites==0)
3745         {
3746                 fprintf(stderr,"\n%c*** ERROR: No transmitter site(s) specified!\n\n",7);
3747                 exit (-1);
3748         }
3749
3750         for (x=0, y=0; x<txsites; x++)
3751         {
3752                 if (tx_site[x].lat==0.0 && tx_site[x].lon==0.0)
3753                 {
3754                         fprintf(stderr,"\n*** ERROR: Transmitter site #%d not found!",x+1);
3755                         y++;
3756                 }
3757         }
3758
3759         if (y)
3760         {
3761                 fprintf(stderr,"%c\n\n",7);
3762                 exit (-1);
3763         }
3764
3765         if ((coverage+LRmap)==0 && rx_site.lat==0.0 && rx_site.lon==0.0)
3766         {
3767                 fprintf(stderr,"\n%c*** ERROR: No receiver site found or specified!\n\n",7);
3768                 exit (-1);
3769         }
3770
3771         /* No errors were detected.  Whew!  :-) */
3772
3773         /* If no SDF path was specified on the command line (-d), check
3774            for a path specified in the $HOME/.splat_path file.  If the
3775            file is not found, then sdf_path[] remains NULL, and the
3776            current working directory is assumed to contain the SDF
3777            files. */
3778
3779         if (sdf_path[0]==0)
3780         {
3781                 sprintf(string,"%s/.splat_path",getenv("HOME"));
3782                 fd=fopen(string,"r");
3783
3784                 if (fd!=NULL)
3785                 {
3786                         fgets(string,253,fd);
3787
3788                         /* Remove <CR> and/or <LF> from string */
3789
3790                         for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0 && x<253; x++);
3791                         string[x]=0;
3792
3793                         strncpy(sdf_path,string,253);
3794
3795                         fclose(fd);
3796                 }
3797         }
3798
3799         /* Ensure a trailing '/' is present in sdf_path */
3800
3801         if (sdf_path[0])
3802         {
3803                 x=strlen(sdf_path);
3804
3805                 if (sdf_path[x-1]!='/' && x!=0)
3806                 {
3807                         sdf_path[x]='/';
3808                         sdf_path[x+1]=0;
3809                 }
3810         }
3811
3812         fprintf(stdout,"%s",header);
3813         fflush(stdout);
3814
3815         x=0;
3816         y=0;
3817
3818         min_lat=0;
3819         max_lat=0;
3820         min_lon=0;
3821         max_lon=0;
3822
3823         rxlat=(unsigned char)floor(rx_site.lat);
3824         rxlon=(unsigned char)floor(rx_site.lon);
3825
3826         if (rxlat!=0)
3827         {
3828                 if (min_lat==0)
3829                         min_lat=rxlat;
3830
3831                 else if (rxlat<min_lat)
3832                         min_lat=rxlat;
3833         }
3834
3835         if (rxlon!=0)
3836         {
3837                 if (min_lon==0)
3838                         min_lon=rxlon;
3839
3840                 else if (rxlon<min_lon)
3841                         min_lon=rxlon;
3842         }
3843
3844         if (rxlat>max_lat)
3845                 max_lat=rxlat;
3846
3847         if (rxlon>max_lon)
3848                 max_lon=rxlon;
3849
3850         for (y=0, z=0; z<txsites; z++)
3851         {
3852                 txlat=(unsigned char)floor(tx_site[z].lat);
3853                 txlon=(unsigned char)floor(tx_site[z].lon);
3854
3855                 if (txlat!=0)
3856                 {
3857                         if (min_lat==0)
3858                                 min_lat=txlat;
3859
3860                         else if (txlat<min_lat)
3861                                 min_lat=txlat;
3862                 }
3863
3864                 if (txlon!=0)
3865                 {
3866                         if (min_lon==0)
3867                                 min_lon=txlon;
3868
3869                         else if (txlon<min_lon)
3870                                 min_lon=txlon;
3871                 }
3872
3873                 if (txlat>max_lat)
3874                         max_lat=txlat;
3875
3876                 if (txlon>max_lon)
3877                         max_lon=txlon;
3878         }
3879
3880         if (min_lat!=0 && min_lon!=0 && max_lat!=0 && max_lon!=0)
3881         {
3882                 for (y=min_lon; y<=max_lon; y++)
3883                         for (x=min_lat; x<=max_lat; x++)
3884                         {
3885                                 sprintf(string,"%u:%u:%u:%u",x, x+1, y, y+1);
3886                                 LoadSDF(string);
3887                         }
3888         }
3889
3890         if (coverage)
3891         {
3892                 for (z=0; z<txsites; z++)
3893                 {
3894                         /* "Ball park" estimates used to load any additional
3895                            SDF files required to conduct this analysis. */
3896
3897                         tx_range=sqrt(1.5*(tx_site[z].alt+GetElevation(tx_site[z])));
3898                         rx_range=sqrt(1.5*altitude);
3899
3900                         /* deg_range determines the maximum
3901                            amount of topo data we read */
3902
3903                         deg_range=(tx_range+rx_range)/69.0;
3904
3905                         /* max_range sets the maximum size of the
3906                            analysis.  A small, non-zero amount can
3907                            be used to shrink the size of the analysis
3908                            and limit the amount of topo data read by
3909                            SPLAT!  A very large number will only increase
3910                            the width of the analysis, not the size of
3911                            the map. */
3912
3913                         if (max_range==0.0)
3914                                 max_range=tx_range+rx_range;
3915
3916                         if (max_range<(tx_range+rx_range))
3917                                 deg_range=max_range/69.0;
3918
3919                         /* Prevent the demand for a really wide coverage
3920                            from allocating more slots than are available
3921                            in memory. */
3922
3923                         switch (MAXSLOTS)
3924                         {
3925                                 case 2: deg_limit=0.25;
3926                                         break;
3927
3928                                 case 4: deg_limit=0.5;
3929                                         break;
3930
3931                                 case 9: deg_limit=1.0;
3932                                         break;
3933
3934                                 case 16: deg_limit=2.0;
3935                                         break;
3936
3937                                 case 25: deg_limit=3.0;
3938                         }
3939
3940                         if (tx_site[z].lat<70.0)
3941                                 deg_range_lon=deg_range/cos(deg2rad*tx_site[z].lat);
3942                         else
3943                                 deg_range_lon=deg_range/cos(deg2rad*70.0);
3944
3945                         /* Correct for squares in degrees not being square in miles */  
3946
3947                         if (deg_range>deg_limit)
3948                                 deg_range=deg_limit;
3949
3950                         if (deg_range_lon>deg_limit)
3951                                 deg_range_lon=deg_limit;
3952
3953                         north_min=(unsigned char)floor(tx_site[z].lat-deg_range);
3954                         north_max=(unsigned char)floor(tx_site[z].lat+deg_range);
3955                         west_min=(unsigned char)floor(tx_site[z].lon-deg_range_lon);
3956                         west_max=(unsigned char)floor(tx_site[z].lon+deg_range_lon);
3957
3958                         if (min_lat==0)
3959                                 min_lat=north_min;
3960
3961                         else if (north_min<min_lat)
3962                                 min_lat=north_min;
3963
3964                         if (min_lon==0)
3965                                 min_lon=west_min;
3966
3967                         else if (west_min<min_lon)
3968                                 min_lon=west_min;
3969
3970                         if (north_max>max_lat)
3971                                 max_lat=north_max;
3972
3973                         if (west_max>max_lon)
3974                                 max_lon=west_max;
3975                 }
3976
3977                 if (min_lat!=0 && min_lon!=0 && max_lat!=0 && max_lon!=0)
3978                 {
3979                         for (y=min_lon; y<=max_lon; y++)
3980                                 for (x=min_lat; x<=max_lat; x++)
3981                                 {
3982                                         sprintf(string,"%u:%u:%u:%u",x, x+1, y, y+1);
3983                                         LoadSDF(string);
3984                                 }
3985                 }
3986         }
3987
3988         if (LRmap)
3989         {
3990                 /* "Ball park" estimates used to load any additional
3991                    SDF files required to conduct this analysis. */
3992
3993                 tx_range=sqrt(1.5*(tx_site[0].alt+GetElevation(tx_site[0])));
3994                 rx_range=sqrt(1.5*altitudeLR);
3995
3996                 /**
3997                 tx_range=sqrt(5.0*tx_site[0].alt);
3998                 rx_range=sqrt(5.0*altitudeLR);
3999                 **/
4000
4001                 /* deg_range determines the maximum
4002                    amount of topo data we read */
4003
4004                 deg_range=(tx_range+rx_range)/69.0;
4005
4006                 /* max_range sets the maximum size of the
4007                    analysis.  A small, non-zero amount can
4008                    be used to shrink the size of the analysis
4009                    and limit the amount of topo data read by
4010                    SPLAT!  A very large number will only increase
4011                    the width of the analysis, not the size of
4012                    the map. */
4013
4014                 if (max_range==0.0)
4015                         max_range=tx_range+rx_range;
4016
4017                 if (max_range<(tx_range+rx_range))
4018                         deg_range=max_range/69.0;
4019
4020                 /* Prevent the demand for a really wide coverage
4021                    from allocating more slots than are available
4022                    in memory. */
4023
4024                 switch (MAXSLOTS)
4025                 {
4026                         case 2: deg_limit=0.25;
4027                                 break;
4028
4029                         case 4: deg_limit=0.5;
4030                                 break;
4031
4032                         case 9: deg_limit=1.0;
4033                                 break;
4034
4035                         case 16: deg_limit=2.0;
4036                                 break;
4037
4038                         case 25: deg_limit=3.0;
4039                 }
4040
4041                 if (tx_site[0].lat<70.0)
4042                         deg_range_lon=deg_range/cos(deg2rad*tx_site[0].lat);
4043                 else
4044                         deg_range_lon=deg_range/cos(deg2rad*70.0);
4045
4046                 /* Correct for squares in degrees not being square in miles */  
4047
4048                 if (deg_range>deg_limit)
4049                         deg_range=deg_limit;
4050
4051                 if (deg_range_lon>deg_limit)
4052                         deg_range_lon=deg_limit;
4053
4054                 north_min=(unsigned char)floor(tx_site[0].lat-deg_range);
4055                 north_max=(unsigned char)floor(tx_site[0].lat+deg_range);
4056                 west_min=(unsigned char)floor(tx_site[0].lon-deg_range_lon);
4057                 west_max=(unsigned char)floor(tx_site[0].lon+deg_range_lon);
4058
4059                 if (min_lat==0)
4060                         min_lat=north_min;
4061
4062                 else if (north_min<min_lat)
4063                         min_lat=north_min;
4064
4065                 if (min_lon==0)
4066                         min_lon=west_min;
4067
4068                 else if (west_min<min_lon)
4069                         min_lon=west_min;
4070
4071                 if (north_max>max_lat)
4072                         max_lat=north_max;
4073
4074                 if (west_max>max_lon)
4075                         max_lon=west_max;
4076
4077                 if (min_lat!=0 && min_lon!=0 && max_lat!=0 && max_lon!=0)
4078                 {
4079                         for (y=min_lon; y<=max_lon; y++)
4080                                 for (x=min_lat; x<=max_lat; x++)
4081                                 {
4082                                         sprintf(string,"%u:%u:%u:%u",x, x+1, y, y+1);
4083                                         LoadSDF(string);
4084                                 }
4085                 }
4086         }
4087
4088         if (mapfile[0])
4089                 map=1;
4090
4091         if (coverage)
4092         {
4093                 for (x=0; x<txsites; x++)
4094                 {
4095                         PlotCoverage(tx_site[x],altitude);
4096                         PlaceMarker(tx_site[x]);
4097
4098                         if (report!='N')
4099                                 SiteReport(tx_site[x]);
4100                 }
4101
4102                 map=1;
4103         }
4104
4105         else if (LRmap)
4106         {
4107                 PlotLRMap(tx_site[0],altitudeLR);
4108                 PlaceMarker(tx_site[0]);
4109
4110                 if (report!='N')
4111                         SiteReport(tx_site[0]);
4112                 
4113                 map=1;
4114         }
4115
4116         else       
4117         {
4118                 PlaceMarker(rx_site);
4119
4120                 for (x=0; x<txsites; x++)
4121                 {
4122                         PlaceMarker(tx_site[x]);
4123
4124                         if (report=='y')
4125                         {
4126                                 switch (x)
4127                                 {
4128                                         case 0:
4129                                                 PlotPath(tx_site[x],rx_site,1);
4130                                                 break;
4131
4132                                         case 1:
4133                                                 PlotPath(tx_site[x],rx_site,8);
4134                                                 break;
4135
4136                                         case 2:
4137                                                 PlotPath(tx_site[x],rx_site,16);
4138                                                 break;
4139
4140                                         case 3:
4141                                                 PlotPath(tx_site[x],rx_site,32);
4142                                 }
4143                         }
4144
4145                         if (report!='N')
4146                                 ObstructionReport(tx_site[x],rx_site,report);
4147                 }
4148         }
4149
4150         if (map)
4151         {
4152                 if (bfs)
4153                 {
4154                         for (x=0; x<bfs; x++)
4155                                 LoadBoundaries(boundary_file[x]);
4156                 }
4157
4158                 if (cities)
4159                 {
4160                         for (x=0; x<cities; x++)
4161                                 LoadCities(city_file[x]);
4162                 }
4163                                 
4164                 if (!LRmap)
4165                         WritePPM(mapfile);
4166                 else
4167                         WritePPMLR(mapfile);
4168         }
4169
4170         if (terrain_plot)
4171         {
4172                 if (txsites>1)
4173                 {
4174                         for (x=0; terrain_file[x]!='.' && terrain_file[x]!=0 && x<80; x++);
4175
4176                         if (terrain_file[x]=='.')  /* extension */
4177                         {
4178                                 ext[0]='.';
4179                                 for (y=1, z=x, x++; terrain_file[x]!=0 && x<253 && y<14; x++, y++)
4180                                         ext[y]=terrain_file[x];
4181
4182                                 ext[y]=0;
4183                                 terrain_file[z]=0;
4184                         }
4185
4186                         else
4187                         {
4188                                 ext[0]=0;  /* No extension */
4189                                 terrain_file[x]=0;
4190                         }
4191
4192                         for (count=0; count<txsites; count++)
4193                         {
4194                                 sprintf(string,"%s-%c%s%c",terrain_file,'1'+count,ext,0);
4195                                 GraphTerrain(tx_site[count],rx_site,string);
4196                         }
4197                 }
4198
4199                 else
4200                         GraphTerrain(tx_site[0],rx_site,terrain_file);
4201         }
4202
4203         if (elevation_plot)
4204         {
4205                 if (txsites>1)
4206                 {
4207                         for (x=0; elevation_file[x]!='.' && elevation_file[x]!=0 && x<80; x++);
4208
4209                         if (elevation_file[x]=='.')  /* extension */
4210                         {
4211                                 ext[0]='.';
4212                                 for (y=1, z=x, x++; elevation_file[x]!=0 && x<253 && y<14; x++, y++)
4213                                         ext[y]=elevation_file[x];
4214
4215                                 ext[y]=0;
4216                                 elevation_file[z]=0;
4217                         }
4218
4219                         else
4220                         {
4221                                 ext[0]=0;  /* No extension */
4222                                 elevation_file[x]=0;
4223                         }
4224
4225                         for (count=0; count<txsites; count++)
4226                         {
4227                                 sprintf(string,"%s-%c%s%c",elevation_file,'1'+count,ext,0);
4228                                 GraphElevation(tx_site[count],rx_site,string);
4229                         }
4230                 }
4231
4232                 else
4233                         GraphElevation(tx_site[0],rx_site,elevation_file);
4234         }
4235
4236         if (height_plot)
4237         {
4238                 if (txsites>1)
4239                 {
4240                         for (x=0; height_file[x]!='.' && height_file[x]!=0 && x<80; x++);
4241
4242                         if (height_file[x]=='.')  /* extension */
4243                         {
4244                                 ext[0]='.';
4245                                 for (y=1, z=x, x++; height_file[x]!=0 && x<253 && y<14; x++, y++)
4246                                         ext[y]=height_file[x];
4247
4248                                 ext[y]=0;
4249                                 height_file[z]=0;
4250                         }
4251
4252                         else
4253                         {
4254                                 ext[0]=0;  /* No extension */
4255                                 height_file[x]=0;
4256                         }
4257
4258                         for (count=0; count<txsites; count++)
4259                         {
4260                                 sprintf(string,"%s-%c%s%c",height_file,'1'+count,ext,0);
4261                                 GraphHeight(tx_site[count],rx_site,string);
4262                         }
4263                 }
4264
4265                 else
4266                         GraphHeight(tx_site[0],rx_site,height_file);
4267         }
4268
4269         if (longley_plot)
4270         {
4271                 if (txsites>1)
4272                 {
4273                         for (x=0; longley_file[x]!='.' && longley_file[x]!=0 && x<80; x++);
4274
4275                         if (longley_file[x]=='.')  /* extension */
4276                         {
4277                                 ext[0]='.';
4278                                 for (y=1, z=x, x++; longley_file[x]!=0 && x<253 && y<14; x++, y++)
4279                                         ext[y]=longley_file[x];
4280
4281                                 ext[y]=0;
4282                                 longley_file[z]=0;
4283                         }
4284
4285                         else
4286                         {
4287                                 ext[0]=0;  /* No extension */
4288                                 longley_file[x]=0;
4289                         }
4290
4291                         for (count=0; count<txsites; count++)
4292                         {
4293                                 sprintf(string,"%s-%c%s%c",longley_file,'1'+count,ext,0);
4294                                 GraphLongley(tx_site[count],rx_site,string);
4295                         }
4296                 }
4297
4298                 else
4299                         GraphLongley(tx_site[0],rx_site,longley_file);
4300         }
4301
4302         return 0;
4303 }