1 /****************************************************************************
2 * SPLAT: An RF Signal Propagation Loss and Terrain Analysis Tool *
3 * Last update: 21-Dec-2006 *
4 *****************************************************************************
5 * Project started in 1997 by John A. Magliacane, KD2BD *
6 *****************************************************************************
8 * Extensively modified by J. D. McDonald in Jan. 2004 to include *
9 * the Longley-Rice propagation model using C++ code from NTIA/ITS. *
11 * See: http://flattop.its.bldrdoc.gov/itm.html *
13 *****************************************************************************
15 * This program is free software; you can redistribute it and/or modify it *
16 * under the terms of the GNU General Public License as published by the *
17 * Free Software Foundation; either version 2 of the License or any later *
20 * This program is distributed in the hope that it will useful, but WITHOUT *
21 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
22 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License *
25 *****************************************************************************
26 g++ -Wall -O3 -s -lm -lbz2 -fomit-frame-pointer itm.cpp splat.cpp -o splat
27 *****************************************************************************/
37 #include "smallfont.h"
41 #define BZBUFFER 65536
44 #define ARRAYSIZE 4950
48 #define ARRAYSIZE 10870
52 #define ARRAYSIZE 19240
56 #define ARRAYSIZE 30025
59 char string[255], sdf_path[255], opened=0, *splat_version={"1.2.0"};
61 double TWOPI=6.283185307179586, HALFPI=1.570796326794896,
62 PI=3.141592653589793, deg2rad=1.74532925199e-02,
63 EARTHRADIUS=20902230.97, METERS_PER_MILE=1609.344,
64 METERS_PER_FOOT=0.3048, KM_PER_MILE=1.609344, earthradius,
67 int min_north=90, max_north=-90, min_west=360, max_west=-1,
68 max_elevation=-32768, min_elevation=32768, bzerror, maxdB=230;
70 unsigned char got_elevation_pattern=0, got_azimuth_pattern=0, metric=0;
72 struct site { double lat;
78 struct path { double lat[ARRAYSIZE];
79 double lon[ARRAYSIZE];
80 double elevation[ARRAYSIZE];
81 double distance[ARRAYSIZE];
85 struct dem { int min_north;
91 short data[1200][1200];
92 unsigned char mask[1200][1200];
95 struct LR { double eps_dielect;
96 double sgm_conductivity;
97 double eno_ns_surfref;
103 float antenna_pattern[361][1001];
106 double elev_l[ARRAYSIZE+10];
108 void point_to_point(double elev[], double tht_m, double rht_m,
109 double eps_dielect, double sgm_conductivity, double eno_ns_surfref,
110 double frq_mhz, int radio_climate, int pol, double conf,
111 double rel, double &dbloss, char *strmode, int &errnum);
113 double arccos(double x, double y)
115 /* This function implements the arc cosine function,
116 returning a value between 0 and TWOPI. */
129 int ReduceAngle(double angle)
131 /* This function normalizes the argument to
132 an integer angle between 0 and 180 degrees */
136 temp=acos(cos(angle*deg2rad));
138 return (int)rint(temp/deg2rad);
141 char *dec2dms(double decimal)
143 /* Converts decimal degrees to degrees, minutes, seconds,
144 (DMS) and returns the result as a character string. */
147 int degrees, minutes, seconds;
175 sprintf(string,"%d%c %d\' %d\"", degrees*sign, 176, minutes, seconds);
179 int OrMask(double lat, double lon, int value)
181 /* Lines, text, markings, and coverage areas are stored in a
182 mask that is combined with topology data when topographic
183 maps are generated by SPLAT!. This function sets bits in
184 the mask based on the latitude and longitude of the area
187 int x, y, indx, minlat, minlon;
190 minlat=(int)floor(lat);
191 minlon=(int)floor(lon);
193 for (indx=0, found=0; indx<MAXSLOTS && found==0;)
194 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
201 x=(int)(1199.0*(lat-floor(lat)));
202 y=(int)(1199.0*(lon-floor(lon)));
204 dem[indx].mask[x][y]|=value;
206 return (dem[indx].mask[x][y]);
213 int GetMask(double lat, double lon)
215 /* This function returns the mask bits based on the latitude
216 and longitude given. */
218 return (OrMask(lat,lon,0));
221 double GetElevation(struct site location)
223 /* This function returns the elevation (in feet) of any location
224 represented by the digital elevation model data in memory.
225 Function returns -5000.0 for locations not found in memory. */
228 int x, y, indx, minlat, minlon;
233 minlat=(int)floor(location.lat);
234 minlon=(int)floor(location.lon);
236 x=(int)(1199.0*(location.lat-floor(location.lat)));
237 y=(int)(1199.0*(location.lon-floor(location.lon)));
239 for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
241 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
244 elevation=3.28084*dem[indx].data[x][y];
252 int AddElevation(double lat, double lon, double height)
254 /* This function adds a user-defined terrain feature
255 (in meters AGL) to the digital elevation model data
256 in memory. Does nothing and returns 0 for locations
257 not found in memory. */
260 int x, y, indx, minlat, minlon;
262 minlat=(int)floor(lat);
263 minlon=(int)floor(lon);
265 x=(int)(1199.0*(lat-floor(lat)));
266 y=(int)(1199.0*(lon-floor(lon)));
268 for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
270 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
273 dem[indx].data[x][y]+=(short)rint(height);
281 double Distance(struct site site1, struct site site2)
283 /* This function returns the great circle distance
284 in miles between any two site locations. */
286 double lat1, lon1, lat2, lon2, distance;
288 lat1=site1.lat*deg2rad;
289 lon1=site1.lon*deg2rad;
290 lat2=site2.lat*deg2rad;
291 lon2=site2.lon*deg2rad;
293 distance=3959.0*acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos((lon1)-(lon2)));
298 double Azimuth(struct site source, struct site destination)
300 /* This function returns the azimuth (in degrees) to the
301 destination as seen from the location of the source. */
303 double dest_lat, dest_lon, src_lat, src_lon,
304 beta, azimuth, diff, num, den, fraction;
306 dest_lat=destination.lat*deg2rad;
307 dest_lon=destination.lon*deg2rad;
309 src_lat=source.lat*deg2rad;
310 src_lon=source.lon*deg2rad;
312 /* Calculate Surface Distance */
314 beta=acos(sin(src_lat)*sin(dest_lat)+cos(src_lat)*cos(dest_lat)*cos(src_lon-dest_lon));
316 /* Calculate Azimuth */
318 num=sin(dest_lat)-(sin(src_lat)*cos(beta));
319 den=cos(src_lat)*sin(beta);
322 /* Trap potential problems in acos() due to rounding */
330 /* Calculate azimuth */
332 azimuth=acos(fraction);
334 /* Reference it to True North */
336 diff=dest_lon-src_lon;
345 azimuth=TWOPI-azimuth;
347 return (azimuth/deg2rad);
350 double ElevationAngle(struct site source, struct site destination)
352 /* This function returns the angle of elevation (in degrees)
353 of the destination as seen from the source location.
354 A positive result represents an angle of elevation (uptilt),
355 while a negative result represents an angle of depression
356 (downtilt), as referenced to a normal to the center of
359 register double a, b, dx;
361 a=GetElevation(destination)+destination.alt+earthradius;
362 b=GetElevation(source)+source.alt+earthradius;
364 dx=5280.0*Distance(source,destination);
366 /* Apply the Law of Cosines */
368 return ((180.0*(acos(((b*b)+(dx*dx)-(a*a))/(2.0*b*dx)))/PI)-90.0);
371 void ReadPath(struct site source, struct site destination)
373 /* This function generates a sequence of latitude and
374 longitude positions between source and destination
375 locations along a great circle path, and stores
376 elevation and distance information for points
377 along that path in the "path" structure. */
380 double azimuth, distance, lat1, lon1, beta, den, num,
381 lat2, lon2, total_distance, x, y, path_length,
383 struct site tempsite;
385 lat1=source.lat*deg2rad;
386 lon1=source.lon*deg2rad;
388 lat2=destination.lat*deg2rad;
389 lon2=destination.lon*deg2rad;
391 azimuth=Azimuth(source,destination)*deg2rad;
393 total_distance=Distance(source,destination);
395 x=68755.0*acos(cos(lon1-lon2)); /* 1200 samples per degree */
396 y=68755.0*acos(cos(lat1-lat2)); /* 68755 samples per radian */
398 path_length=sqrt((x*x)+(y*y)); /* Total number of samples */
400 increment=total_distance/path_length; /* Miles per sample */
402 for (distance=0, c=0; distance<=total_distance; distance+=increment)
404 beta=distance/3959.0;
405 lat2=asin(sin(lat1)*cos(beta)+cos(azimuth)*sin(beta)*cos(lat1));
406 num=cos(beta)-(sin(lat1)*sin(lat2));
407 den=cos(lat1)*cos(lat2);
409 if (azimuth==0.0 && (beta>HALFPI-lat1))
412 else if (azimuth==HALFPI && (beta>HALFPI+lat1))
415 else if (fabs(num/den)>1.0)
420 if ((PI-azimuth)>=0.0)
421 lon2=lon1-arccos(num,den);
423 lon2=lon1+arccos(num,den);
441 path.elevation[c]=GetElevation(tempsite);
442 path.distance[c]=distance;
447 /* Make sure exact destination point is recorded at path.length-1 */
451 path.lat[c]=destination.lat;
452 path.lon[c]=destination.lon;
453 path.elevation[c]=GetElevation(destination);
454 path.distance[c]=total_distance;
461 path.length=ARRAYSIZE-1;
464 double ElevationAngle2(struct site source, struct site destination, double er)
466 /* This function returns the angle of elevation (in degrees)
467 of the destination as seen from the source location, UNLESS
468 the path between the sites is obstructed, in which case, the
469 elevation angle to the first obstruction is returned instead.
470 "er" represents the earth radius. */
474 double source_alt, destination_alt, cos_xmtr_angle,
475 cos_test_angle, test_alt, elevation, distance,
476 source_alt2, first_obstruction_angle=0.0;
481 ReadPath(source,destination);
483 distance=5280.0*Distance(source,destination);
484 source_alt=er+source.alt+GetElevation(source);
485 destination_alt=er+destination.alt+GetElevation(destination);
486 source_alt2=source_alt*source_alt;
488 /* Calculate the cosine of the elevation angle of the
489 destination (receiver) as seen by the source (transmitter). */
491 cos_xmtr_angle=((source_alt2)+(distance*distance)-(destination_alt*destination_alt))/(2.0*source_alt*distance);
493 /* Test all points in between source and destination locations to
494 see if the angle to a topographic feature generates a higher
495 elevation angle than that produced by the destination. Begin
496 at the source since we're interested in identifying the FIRST
497 obstruction along the path between source and destination. */
499 for (x=2, block=0; x<path.length && block==0; x++)
501 distance=5280.0*path.distance[x];
503 test_alt=earthradius+path.elevation[x];
505 cos_test_angle=((source_alt2)+(distance*distance)-(test_alt*test_alt))/(2.0*source_alt*distance);
507 /* Compare these two angles to determine if
508 an obstruction exists. Since we're comparing
509 the cosines of these angles rather than
510 the angles themselves, the sense of the
511 following "if" statement is reversed from
512 what it would be if the angles themselves
515 if (cos_xmtr_angle>cos_test_angle)
518 first_obstruction_angle=((acos(cos_test_angle))/deg2rad)-90.0;
523 elevation=first_obstruction_angle;
526 elevation=((acos(cos_xmtr_angle))/deg2rad)-90.0;
533 double AverageTerrain(struct site source, double azimuthx, double start_distance, double end_distance)
535 /* This function returns the average terrain calculated in
536 the direction of "azimuth" (degrees) between "start_distance"
537 and "end_distance" (miles) from the source location. If
538 the terrain is all water (non-critical error), -5000.0 is
539 returned. If not enough SDF data has been loaded into
540 memory to complete the survey (critical error), then
541 -9999.0 is returned. */
543 int c, samples, endpoint;
544 double beta, lat1, lon1, lat2, lon2, num, den, azimuth, terrain=0.0;
545 struct site destination;
547 lat1=source.lat*deg2rad;
548 lon1=source.lon*deg2rad;
550 /* Generate a path of elevations between the source
551 location and the remote location provided. */
553 beta=end_distance/3959.0;
555 azimuth=deg2rad*azimuthx;
557 lat2=asin(sin(lat1)*cos(beta)+cos(azimuth)*sin(beta)*cos(lat1));
558 num=cos(beta)-(sin(lat1)*sin(lat2));
559 den=cos(lat1)*cos(lat2);
561 if (azimuth==0.0 && (beta>HALFPI-lat1))
564 else if (azimuth==HALFPI && (beta>HALFPI+lat1))
567 else if (fabs(num/den)>1.0)
572 if ((PI-azimuth)>=0.0)
573 lon2=lon1-arccos(num,den);
575 lon2=lon1+arccos(num,den);
587 destination.lat=lat2;
588 destination.lon=lon2;
590 /* If SDF data is missing for the endpoint of
591 the radial, then the average terrain cannot
592 be accurately calculated. Return -9999.0 */
594 if (GetElevation(destination)<-4999.0)
598 ReadPath(source,destination);
600 endpoint=path.length;
602 /* Shrink the length of the radial if the
603 outermost portion is not over U.S. land. */
605 for (c=endpoint-1; c>=0 && path.elevation[c]==0.0; c--);
609 for (c=0, samples=0; c<endpoint; c++)
611 if (path.distance[c]>=start_distance)
613 terrain+=path.elevation[c];
619 terrain=-5000.0; /* No land */
621 terrain=(terrain/(double)samples);
627 double haat(struct site antenna)
629 /* This function returns the antenna's Height Above Average
630 Terrain (HAAT) based on FCC Part 73.313(d). If a critical
631 error occurs, such as a lack of SDF data to complete the
632 survey, -5000.0 is returned. */
636 double terrain, avg_terrain, haat, sum=0.0;
638 /* Calculate the average terrain between 2 and 10 miles
639 from the antenna site at azimuths of 0, 45, 90, 135,
640 180, 225, 270, and 315 degrees. */
642 for (c=0, azi=0; azi<=315 && error==0; azi+=45)
644 terrain=AverageTerrain(antenna, (double)azi, 2.0, 10.0);
646 if (terrain<-9998.0) /* SDF data is missing */
649 if (terrain>-4999.0) /* It's land, not water */
651 sum+=terrain; /* Sum of averages */
660 avg_terrain=(sum/(double)c);
661 haat=(antenna.alt+GetElevation(antenna))-avg_terrain;
666 float LonDiff(float lon1, float lon2)
668 /* This function returns the short path longitudinal
669 difference between longitude1 and longitude2
670 as an angle between -180.0 and +180.0 degrees.
671 If lon1 is west of lon2, the result is positive.
672 If lon1 is east of lon2, the result is negative. */
687 void PlaceMarker(struct site location)
689 /* This function places text and marker data in the mask array
690 for illustration on topographic maps generated by SPLAT!.
691 By default, SPLAT! centers text information BELOW the marker,
692 but may move it above, to the left, or to the right of the
693 marker depending on how much room is available on the map,
694 or depending on whether the area is already occupied by
695 another marker or label. If no room or clear space is
696 available on the map to place the marker and its associated
697 text, then the marker and text are not written to the map. */
700 char ok2print, occupied;
701 double x, y, lat, lon, textx=0.0, texty=0.0, xmin, xmax,
702 ymin, ymax, p1, p3, p6, p8, p12, p16, p24, label_length;
711 if (lat<xmax && lat>xmin && (LonDiff(lon,ymax)<0.0) && (LonDiff(lon,ymin)>0.0))
723 /* Is Marker Position Clear Of Text Or Other Markers? */
725 for (x=lat-p3; (x<=xmax && x>=xmin && x<=lat+p3); x+=p1)
726 for (y=lon-p3; (LonDiff(y,ymax)<=0.0) && (LonDiff(y,ymin)>=0.0) && (LonDiff(y,lon+p3)<=0.0); y+=p1)
727 occupied|=(GetMask(x,y)&2);
731 /* Determine Where Text Can Be Positioned */
733 /* label_length=length in pixels.
734 Each character is 8 pixels wide. */
736 label_length=p1*(double)(strlen(location.name)<<3);
738 if ((LonDiff(lon+label_length,ymax)<=0.0) && (LonDiff(lon-label_length,ymin)>=0.0))
740 /* Default: Centered Text */
742 texty=lon+label_length/2.0;
746 /* Position Text Below The Marker */
753 /* Is This Position Clear Of
754 Text Or Other Markers? */
756 for (a=0, occupied=0; a<16; a++)
758 for (b=0; b<(int)strlen(location.name); b++)
759 for (c=0; c<8; c++, y-=p1)
760 occupied|=(GetMask(x,y)&2);
774 /* Position Text Above The Marker */
781 /* Is This Position Clear Of
782 Text Or Other Markers? */
784 for (a=0, occupied=0; a<16; a++)
786 for (b=0; b<(int)strlen(location.name); b++)
787 for (c=0; c<8; c++, y-=p1)
788 occupied|=(GetMask(x,y)&2);
803 if (LonDiff(lon-label_length,ymin)>=0.0)
805 /* Position Text To The
806 Right Of The Marker */
814 /* Is This Position Clear Of
815 Text Or Other Markers? */
817 for (a=0, occupied=0; a<16; a++)
819 for (b=0; b<(int)strlen(location.name); b++)
820 for (c=0; c<8; c++, y-=p1)
821 occupied|=(GetMask(x,y)&2);
835 /* Position Text To The
836 Left Of The Marker */
839 texty=lon+p8+(label_length);
844 /* Is This Position Clear Of
845 Text Or Other Markers? */
847 for (a=0, occupied=0; a<16; a++)
849 for (b=0; b<(int)strlen(location.name); b++)
850 for (c=0; c<8; c++, y-=p1)
851 occupied|=(GetMask(x,y)&2);
864 /* textx and texty contain the latitude and longitude
865 coordinates that describe the placement of the text
875 for (a=0; a<16 && ok2print; a++)
877 for (b=0; b<(int)strlen(location.name); b++)
879 byte=fontdata[16*(location.name[b])+a];
881 for (c=128; c>0; c=c>>1, y-=p1)
890 /* Draw Square Marker Centered
891 On Location Specified */
893 for (x=lat-p3; (x<=xmax && x>=xmin && x<=lat+p3); x+=p1)
894 for (y=lon-p3; (LonDiff(y,ymax)<=0.0) && (LonDiff(y,ymin)>=0.0) && (LonDiff(y,lon+p3)<=0.0); y+=p1)
901 double ReadBearing(char *input)
903 /* This function takes numeric input in the form of a character
904 string, and returns an equivalent bearing in degrees as a
905 decimal number (double). The input may either be expressed
906 in decimal format (40.139722) or degree, minute, second
907 format (40 08 23). This function also safely handles
908 extra spaces found either leading, trailing, or
909 embedded within the numbers expressed in the
910 input string. Decimal seconds are permitted. */
912 double seconds, bearing=0.0;
914 int a, b, length, degrees, minutes;
916 /* Copy "input" to "string", and ignore any extra
917 spaces that might be present in the process. */
920 length=strlen(input);
922 for (a=0, b=0; a<length && a<18; a++)
924 if ((input[a]!=32 && input[a]!='\n') || (input[a]==32 && input[a+1]!=32 && input[a+1]!='\n' && b!=0))
933 /* Count number of spaces in the clean string. */
935 length=strlen(string);
937 for (a=0, b=0; a<length; a++)
941 if (b==0) /* Decimal Format (40.139722) */
942 sscanf(string,"%lf",&bearing);
944 if (b==2) /* Degree, Minute, Second Format (40 08 23) */
946 sscanf(string,"%d %d %lf",°rees, &minutes, &seconds);
948 bearing=(double)abs(degrees)+((double)abs(minutes)/60)+(fabs(seconds)/3600);
950 if ((degrees<0) || (minutes<0) || (seconds<0.0))
954 /* Anything else returns a 0.0 */
956 if (bearing>360.0 || bearing<-90.0)
962 struct site LoadQTH(char *filename)
964 /* This function reads SPLAT! .qth (site location) files.
965 The latitude and longitude may be expressed either in
966 decimal degrees, or in degree, minute, second format.
967 Antenna height is assumed to be expressed in feet above
968 ground level (AGL), unless followed by the letter 'M',
969 or 'm', or by the word "meters" or "Meters", in which
970 case meters is assumed, and is handled accordingly. */
973 char string[50], qthfile[255];
974 struct site tempsite;
977 for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
978 qthfile[x]=filename[x];
991 fd=fopen(qthfile,"r");
998 /* Strip <CR> and/or <LF> from end of site name */
1000 for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0; tempsite.name[x]=string[x], x++);
1005 fgets(string,49,fd);
1006 tempsite.lat=ReadBearing(string);
1008 /* Site Longitude */
1009 fgets(string,49,fd);
1010 tempsite.lon=ReadBearing(string);
1012 /* Antenna Height */
1013 fgets(string,49,fd);
1016 /* Remove <CR> and/or <LF> from antenna height string */
1017 for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0; x++);
1021 /* Antenna height may either be in feet or meters.
1022 If the letter 'M' or 'm' is discovered in
1023 the string, then this is an indication that
1024 the value given is expressed in meters, and
1025 must be converted to feet before exiting. */
1027 for (x=0; string[x]!='M' && string[x]!='m' && string[x]!=0 && x<48; x++);
1028 if (string[x]=='M' || string[x]=='m')
1031 sscanf(string,"%f",&tempsite.alt);
1032 tempsite.alt*=3.28084;
1038 sscanf(string,"%f",&tempsite.alt);
1045 void LoadPAT(char *filename)
1047 /* This function reads and processes antenna pattern (.az
1048 and .el) files that correspond in name to previously
1049 loaded SPLAT! .lrp files. */
1051 int a, b, w, x, y, z, last_index, next_index, span;
1052 char string[255], azfile[255], elfile[255], *pointer=NULL;
1053 float az, xx, elevation, amplitude, rotation, valid1, valid2,
1054 delta, azimuth[361], azimuth_pattern[361], el_pattern[10001],
1055 elevation_pattern[361][1001], slant_angle[361], tilt,
1056 mechanical_tilt, tilt_azimuth, tilt_increment, sum;
1058 unsigned char read_count[10001];
1060 for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
1062 azfile[x]=filename[x];
1063 elfile[x]=filename[x];
1078 /* Load .az antenna pattern file */
1080 fd=fopen(azfile,"r");
1084 /* Clear azimuth pattern array */
1086 for (x=0; x<=360; x++)
1093 /* Read azimuth pattern rotation
1094 in degrees measured clockwise
1097 fgets(string,254,fd);
1098 pointer=strchr(string,';');
1103 sscanf(string,"%f",&rotation);
1106 /* Read azimuth (degrees) and corresponding
1107 normalized field radiation pattern amplitude
1108 (0.0 to 1.0) until EOF is reached. */
1110 fgets(string,254,fd);
1111 pointer=strchr(string,';');
1116 sscanf(string,"%f %f",&az, &litude);
1122 if (x>=0 && x<=360 && fd!=NULL)
1124 azimuth[x]+=amplitude;
1128 fgets(string,254,fd);
1129 pointer=strchr(string,';');
1134 sscanf(string,"%f %f",&az, &litude);
1136 } while (feof(fd)==0);
1141 /* Handle 0=360 degree ambiguity */
1143 if ((read_count[0]==0) && (read_count[360]!=0))
1145 read_count[0]=read_count[360];
1146 azimuth[0]=azimuth[360];
1149 if ((read_count[0]!=0) && (read_count[360]==0))
1151 read_count[360]=read_count[0];
1152 azimuth[360]=azimuth[0];
1155 /* Average pattern values in case more than
1156 one was read for each degree of azimuth. */
1158 for (x=0; x<=360; x++)
1160 if (read_count[x]>1)
1161 azimuth[x]/=(float)read_count[x];
1164 /* Interpolate missing azimuths
1165 to completely fill the array */
1170 for (x=0; x<=360; x++)
1172 if (read_count[x]!=0)
1180 if (last_index!=-1 && next_index!=-1)
1182 valid1=azimuth[last_index];
1183 valid2=azimuth[next_index];
1185 span=next_index-last_index;
1186 delta=(valid2-valid1)/(float)span;
1188 for (y=last_index+1; y<next_index; y++)
1189 azimuth[y]=azimuth[y-1]+delta;
1196 /* Perform azimuth pattern rotation
1197 and load azimuth_pattern[361] with
1198 azimuth pattern data in its final form. */
1200 for (x=0; x<360; x++)
1202 y=x+(int)rintf(rotation);
1207 azimuth_pattern[y]=azimuth[x];
1210 azimuth_pattern[360]=azimuth_pattern[0];
1212 got_azimuth_pattern=255;
1215 /* Read and process .el file */
1217 fd=fopen(elfile,"r");
1221 for (x=0; x<=10000; x++)
1227 /* Read mechanical tilt (degrees) and
1228 tilt azimuth in degrees measured
1229 clockwise from true North. */
1231 fgets(string,254,fd);
1232 pointer=strchr(string,';');
1237 sscanf(string,"%f %f",&mechanical_tilt, &tilt_azimuth);
1239 /* Read elevation (degrees) and corresponding
1240 normalized field radiation pattern amplitude
1241 (0.0 to 1.0) until EOF is reached. */
1243 fgets(string,254,fd);
1244 pointer=strchr(string,';');
1249 sscanf(string,"%f %f", &elevation, &litude);
1253 /* Read in normalized radiated field values
1254 for every 0.01 degrees of elevation between
1255 -10.0 and +90.0 degrees */
1257 x=(int)rintf(100.0*(elevation+10.0));
1259 if (x>=0 && x<=10000)
1261 el_pattern[x]+=amplitude;
1265 fgets(string,254,fd);
1266 pointer=strchr(string,';');
1271 sscanf(string,"%f %f", &elevation, &litude);
1276 /* Average the field values in case more than
1277 one was read for each 0.01 degrees of elevation. */
1279 for (x=0; x<=10000; x++)
1281 if (read_count[x]>1)
1282 el_pattern[x]/=(float)read_count[x];
1285 /* Interpolate between missing elevations (if
1286 any) to completely fill the array and provide
1287 radiated field values for every 0.01 degrees of
1293 for (x=0; x<=10000; x++)
1295 if (read_count[x]!=0)
1303 if (last_index!=-1 && next_index!=-1)
1305 valid1=el_pattern[last_index];
1306 valid2=el_pattern[next_index];
1308 span=next_index-last_index;
1309 delta=(valid2-valid1)/(float)span;
1311 for (y=last_index+1; y<next_index; y++)
1312 el_pattern[y]=el_pattern[y-1]+delta;
1319 /* Fill slant_angle[] array with offset angles based
1320 on the antenna's mechanical beam tilt (if any)
1321 and tilt direction (azimuth). */
1323 if (mechanical_tilt==0.0)
1325 for (x=0; x<=360; x++)
1331 tilt_increment=mechanical_tilt/90.0;
1333 for (x=0; x<=360; x++)
1336 y=(int)rintf(tilt_azimuth+xx);
1345 slant_angle[y]=-(tilt_increment*(90.0-xx));
1348 slant_angle[y]=-(tilt_increment*(xx-270.0));
1352 slant_angle[360]=slant_angle[0]; /* 360 degree wrap-around */
1354 for (w=0; w<=360; w++)
1356 tilt=slant_angle[w];
1358 /** Convert tilt angle to
1359 an array index offset **/
1361 y=(int)rintf(100.0*tilt);
1363 /* Copy shifted el_pattern[10001] field
1364 values into elevation_pattern[361][1001]
1365 at the corresponding azimuth, downsampling
1366 (averaging) along the way in chunks of 10. */
1368 for (x=y, z=0; z<=1000; x+=10, z++)
1370 for (sum=0.0, a=0; a<10; a++)
1374 if (b>=0 && b<=10000)
1379 sum+=el_pattern[10000];
1382 elevation_pattern[w][z]=sum/10.0;
1386 got_elevation_pattern=255;
1389 for (x=0; x<=360; x++)
1391 for (y=0; y<=1000; y++)
1393 if (got_elevation_pattern)
1394 elevation=elevation_pattern[x][y];
1398 if (got_azimuth_pattern)
1399 az=azimuth_pattern[x];
1403 LR.antenna_pattern[x][y]=az*elevation;
1408 int LoadSDF_SDF(char *name)
1410 /* This function reads uncompressed SPLAT Data Files (.sdf)
1411 containing digital elevation model data into memory.
1412 Elevation data, maximum and minimum elevations, and
1413 quadrangle limits are stored in the first available
1416 int x, y, data, indx, minlat, minlon, maxlat, maxlon;
1417 char found, free_slot=0, line[20], sdf_file[255],
1418 path_plus_name[255];
1421 for (x=0; name[x]!='.' && name[x]!=0 && x<250; x++)
1422 sdf_file[x]=name[x];
1426 /* Parse filename for minimum latitude and longitude values */
1428 sscanf(sdf_file,"%d:%d:%d:%d",&minlat,&maxlat,&minlon,&maxlon);
1436 /* Is it already in memory? */
1438 for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
1440 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west && maxlat==dem[indx].max_north && maxlon==dem[indx].max_west)
1444 /* Is room available to load it? */
1448 for (indx=0, free_slot=0; indx<MAXSLOTS && free_slot==0; indx++)
1449 if (dem[indx].max_north==-90)
1455 if (free_slot && found==0 && indx>=0 && indx<MAXSLOTS)
1457 /* Search for SDF file in current working directory first */
1459 strncpy(path_plus_name,sdf_file,255);
1461 fd=fopen(path_plus_name,"rb");
1465 /* Next, try loading SDF file from path specified
1466 in $HOME/.splat_path file or by -d argument */
1468 strncpy(path_plus_name,sdf_path,255);
1469 strncat(path_plus_name,sdf_file,255);
1471 fd=fopen(path_plus_name,"rb");
1476 fprintf(stdout,"Loading \"%s\" into slot %d...",path_plus_name,indx+1);
1480 sscanf(line,"%d",&dem[indx].max_west);
1483 sscanf(line,"%d",&dem[indx].min_north);
1486 sscanf(line,"%d",&dem[indx].min_west);
1489 sscanf(line,"%d",&dem[indx].max_north);
1491 for (x=0; x<1200; x++)
1492 for (y=0; y<1200; y++)
1495 sscanf(line,"%d",&data);
1497 dem[indx].data[x][y]=data;
1499 if (data>dem[indx].max_el)
1500 dem[indx].max_el=data;
1502 if (data<dem[indx].min_el)
1503 dem[indx].min_el=data;
1508 if (dem[indx].min_el<min_elevation)
1509 min_elevation=dem[indx].min_el;
1511 if (dem[indx].max_el>max_elevation)
1512 max_elevation=dem[indx].max_el;
1515 max_north=dem[indx].max_north;
1517 else if (dem[indx].max_north>max_north)
1518 max_north=dem[indx].max_north;
1521 min_north=dem[indx].min_north;
1523 else if (dem[indx].min_north<min_north)
1524 min_north=dem[indx].min_north;
1527 max_west=dem[indx].max_west;
1531 if (abs(dem[indx].max_west-max_west)<180)
1533 if (dem[indx].max_west>max_west)
1534 max_west=dem[indx].max_west;
1539 if (dem[indx].max_west<max_west)
1540 max_west=dem[indx].max_west;
1545 min_west=dem[indx].min_west;
1549 if (abs(dem[indx].min_west-min_west)<180)
1551 if (dem[indx].min_west<min_west)
1552 min_west=dem[indx].min_west;
1557 if (dem[indx].min_west>min_west)
1558 min_west=dem[indx].min_west;
1562 fprintf(stdout," Done!\n");
1575 char *BZfgets(BZFILE *bzfd, unsigned length)
1577 /* This function returns at most one less than 'length' number
1578 of characters from a bz2 compressed file whose file descriptor
1579 is pointed to by *bzfd. In operation, a buffer is filled with
1580 uncompressed data (size = BZBUFFER), which is then parsed
1581 and doled out as NULL terminated character strings every time
1582 this function is invoked. A NULL string indicates an EOF
1583 or error condition. */
1585 static int x, y, nBuf;
1586 static char buffer[BZBUFFER+1], output[BZBUFFER+1];
1589 if (opened!=1 && bzerror==BZ_OK)
1591 /* First time through. Initialize everything! */
1602 if (x==nBuf && bzerror!=BZ_STREAM_END && bzerror==BZ_OK && opened)
1604 /* Uncompress data into a static buffer */
1606 nBuf=BZ2_bzRead(&bzerror, bzfd, buffer, BZBUFFER);
1611 /* Build a string from buffer contents */
1613 output[y]=buffer[x];
1615 if (output[y]=='\n' || output[y]==0 || y==(int)length-1)
1634 int LoadSDF_BZ(char *name)
1636 /* This function reads .bz2 compressed SPLAT Data Files containing
1637 digital elevation model data into memory. Elevation data,
1638 maximum and minimum elevations, and quadrangle limits are
1639 stored in the first available dem[] structure. */
1641 int x, y, data, indx, minlat, minlon, maxlat, maxlon;
1642 char found, free_slot=0, sdf_file[255], path_plus_name[255];
1646 for (x=0; name[x]!='.' && name[x]!=0 && x<247; x++)
1647 sdf_file[x]=name[x];
1651 /* Parse sdf_file name for minimum latitude and longitude values */
1653 sscanf(sdf_file,"%d:%d:%d:%d",&minlat,&maxlat,&minlon,&maxlon);
1665 /* Is it already in memory? */
1667 for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
1669 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west && maxlat==dem[indx].max_north && maxlon==dem[indx].max_west)
1673 /* Is room available to load it? */
1677 for (indx=0, free_slot=0; indx<MAXSLOTS && free_slot==0; indx++)
1678 if (dem[indx].max_north==-90)
1684 if (free_slot && found==0 && indx>=0 && indx<MAXSLOTS)
1686 /* Search for SDF file in current working directory first */
1688 strncpy(path_plus_name,sdf_file,255);
1690 fd=fopen(path_plus_name,"rb");
1691 bzfd=BZ2_bzReadOpen(&bzerror,fd,0,0,NULL,0);
1693 if (fd==NULL || bzerror!=BZ_OK)
1695 /* Next, try loading SDF file from path specified
1696 in $HOME/.splat_path file or by -d argument */
1698 strncpy(path_plus_name,sdf_path,255);
1699 strncat(path_plus_name,sdf_file,255);
1701 fd=fopen(path_plus_name,"rb");
1702 bzfd=BZ2_bzReadOpen(&bzerror,fd,0,0,NULL,0);
1705 if (fd!=NULL && bzerror==BZ_OK)
1707 fprintf(stdout,"Loading \"%s\" into slot %d...",path_plus_name,indx+1);
1710 sscanf(BZfgets(bzfd,255),"%d",&dem[indx].max_west);
1711 sscanf(BZfgets(bzfd,255),"%d",&dem[indx].min_north);
1712 sscanf(BZfgets(bzfd,255),"%d",&dem[indx].min_west);
1713 sscanf(BZfgets(bzfd,255),"%d",&dem[indx].max_north);
1715 for (x=0; x<1200; x++)
1716 for (y=0; y<1200; y++)
1718 sscanf(BZfgets(bzfd,20),"%d",&data);
1720 dem[indx].data[x][y]=data;
1722 if (data>dem[indx].max_el)
1723 dem[indx].max_el=data;
1725 if (data<dem[indx].min_el)
1726 dem[indx].min_el=data;
1731 BZ2_bzReadClose(&bzerror,bzfd);
1733 if (dem[indx].min_el<min_elevation)
1734 min_elevation=dem[indx].min_el;
1736 if (dem[indx].max_el>max_elevation)
1737 max_elevation=dem[indx].max_el;
1740 max_north=dem[indx].max_north;
1742 else if (dem[indx].max_north>max_north)
1743 max_north=dem[indx].max_north;
1746 min_north=dem[indx].min_north;
1748 else if (dem[indx].min_north<min_north)
1749 min_north=dem[indx].min_north;
1752 max_west=dem[indx].max_west;
1756 if (abs(dem[indx].max_west-max_west)<180)
1758 if (dem[indx].max_west>max_west)
1759 max_west=dem[indx].max_west;
1764 if (dem[indx].max_west<max_west)
1765 max_west=dem[indx].max_west;
1770 min_west=dem[indx].min_west;
1774 if (abs(dem[indx].min_west-min_west)<180)
1776 if (dem[indx].min_west<min_west)
1777 min_west=dem[indx].min_west;
1782 if (dem[indx].min_west>min_west)
1783 min_west=dem[indx].min_west;
1787 fprintf(stdout," Done!\n");
1800 char LoadSDF(char *name)
1802 /* This function loads the requested SDF file from the filesystem.
1803 It first tries to invoke the LoadSDF_SDF() function to load an
1804 uncompressed SDF file (since uncompressed files load slightly
1805 faster). If that attempt fails, then it tries to load a
1806 compressed SDF file by invoking the LoadSDF_BZ() function.
1807 If that fails, then we can assume that no elevation data
1808 exists for the region requested, and that the region
1809 requested must be entirely over water. */
1811 int x, y, indx, minlat, minlon, maxlat, maxlon;
1812 char found, free_slot=0;
1813 int return_value=-1;
1815 /* Try to load an uncompressed SDF first. */
1817 return_value=LoadSDF_SDF(name);
1819 /* If that fails, try loading a compressed SDF. */
1821 if (return_value==0 || return_value==-1)
1822 return_value=LoadSDF_BZ(name);
1824 /* If neither format can be found, then assume the area is water. */
1826 if (return_value==0 || return_value==-1)
1828 /* Parse SDF name for minimum latitude and longitude values */
1830 sscanf(name,"%d:%d:%d:%d",&minlat,&maxlat,&minlon,&maxlon);
1832 /* Is it already in memory? */
1834 for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
1836 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west && maxlat==dem[indx].max_north && maxlon==dem[indx].max_west)
1840 /* Is room available to load it? */
1844 for (indx=0, free_slot=0; indx<MAXSLOTS && free_slot==0; indx++)
1845 if (dem[indx].max_north==-90)
1851 if (free_slot && found==0 && indx>=0 && indx<MAXSLOTS)
1853 fprintf(stdout,"Region \"%s\" assumed as sea-level into slot %d...",name,indx+1);
1856 dem[indx].max_west=maxlon;
1857 dem[indx].min_north=minlat;
1858 dem[indx].min_west=minlon;
1859 dem[indx].max_north=maxlat;
1861 /* Fill DEM with sea-level topography */
1863 for (x=0; x<1200; x++)
1864 for (y=0; y<1200; y++)
1866 dem[indx].data[x][y]=0;
1868 if (dem[indx].min_el>0)
1872 if (dem[indx].min_el<min_elevation)
1873 min_elevation=dem[indx].min_el;
1875 if (dem[indx].max_el>max_elevation)
1876 max_elevation=dem[indx].max_el;
1879 max_north=dem[indx].max_north;
1881 else if (dem[indx].max_north>max_north)
1882 max_north=dem[indx].max_north;
1885 min_north=dem[indx].min_north;
1887 else if (dem[indx].min_north<min_north)
1888 min_north=dem[indx].min_north;
1891 max_west=dem[indx].max_west;
1895 if (abs(dem[indx].max_west-max_west)<180)
1897 if (dem[indx].max_west>max_west)
1898 max_west=dem[indx].max_west;
1903 if (dem[indx].max_west<max_west)
1904 max_west=dem[indx].max_west;
1909 min_west=dem[indx].min_west;
1913 if (abs(dem[indx].min_west-min_west)<180)
1915 if (dem[indx].min_west<min_west)
1916 min_west=dem[indx].min_west;
1921 if (dem[indx].min_west>min_west)
1922 min_west=dem[indx].min_west;
1926 fprintf(stdout," Done!\n");
1933 return return_value;
1936 void LoadCities(char *filename)
1938 /* This function reads SPLAT! city/site files, and plots
1939 the locations and names of the cities and site locations
1940 read on topographic maps generated by SPLAT! */
1943 char input[80], str[3][80];
1944 struct site city_site;
1947 fd=fopen(filename,"r");
1953 fprintf(stdout,"Reading \"%s\"... ",filename);
1956 while (fd!=NULL && feof(fd)==0)
1958 /* Parse line for name, latitude, and longitude */
1960 for (x=0, y=0, z=0; x<78 && input[x]!=0 && z<3; x++)
1962 if (input[x]!=',' && y<78)
1976 strncpy(city_site.name,str[0],49);
1977 city_site.lat=ReadBearing(str[1]);
1978 city_site.lon=ReadBearing(str[2]);
1981 PlaceMarker(city_site);
1987 fprintf(stdout,"Done!\n");
1992 fprintf(stderr,"*** ERROR: \"%s\": not found!\n",filename);
1995 void LoadUDT(char *filename)
1997 /* This function reads a file containing User-Defined Terrain
1998 features for their addition to the digital elevation model
1999 data used by SPLAT!. Elevations in the UDT file are evaluated
2000 and then copied into a temporary file under /tmp. Then the
2001 contents of the temp file are scanned, and if found to be unique,
2002 are added to the ground elevations described by the digital
2003 elevation data already loaded into memory. */
2005 int i, x, y, z, fd=0;
2006 char input[80], str[3][80], tempname[15], *pointer=NULL;
2007 double latitude, longitude, height, templat, templon,
2008 tempheight, one_pixel;
2009 FILE *fd1=NULL, *fd2=NULL;
2011 strcpy(tempname,"/tmp/XXXXXX\0");
2012 one_pixel=1.0/1200.0;
2014 fd1=fopen(filename,"r");
2018 fd=mkstemp(tempname);
2019 fd2=fopen(tempname,"w");
2021 fgets(input,78,fd1);
2023 pointer=strchr(input,';');
2028 fprintf(stdout,"Reading \"%s\"... ",filename);
2031 while (feof(fd1)==0)
2033 /* Parse line for latitude, longitude, height */
2035 for (x=0, y=0, z=0; x<78 && input[x]!=0 && z<3; x++)
2037 if (input[x]!=',' && y<78)
2051 latitude=ReadBearing(str[0]);
2052 longitude=ReadBearing(str[1]);
2054 /* Remove <CR> and/or <LF> from antenna height string */
2056 for (i=0; str[2][i]!=13 && str[2][i]!=10 && str[2][i]!=0; i++);
2060 /* The terrain feature may be expressed in either
2061 feet or meters. If the letter 'M' or 'm' is
2062 discovered in the string, then this is an
2063 indication that the value given is expressed
2064 in meters. Otherwise the height is interpreted
2065 as being expressed in feet. */
2067 for (i=0; str[2][i]!='M' && str[2][i]!='m' && str[2][i]!=0 && i<48; i++);
2069 if (str[2][i]=='M' || str[2][i]=='m')
2072 height=rint(atof(str[2]));
2078 height=rint(3.28084*atof(str[2]));
2082 fprintf(fd2,"%f, %f, %f\n",latitude, longitude, height);
2084 fgets(input,78,fd1);
2086 pointer=strchr(input,';');
2096 fprintf(stdout,"Done!\n");
2099 fd1=fopen(tempname,"r");
2100 fd2=fopen(tempname,"r");
2102 fscanf(fd1,"%lf, %lf, %lf", &latitude, &longitude, &height);
2104 for (y=0; feof(fd1)==0; y++)
2108 fscanf(fd2,"%lf, %lf, %lf", &templat, &templon, &tempheight);
2110 for (x=0, z=0; feof(fd2)==0; x++)
2113 if (fabs(latitude-templat)<=one_pixel && fabs(longitude-templon)<=one_pixel)
2116 fscanf(fd2,"%lf, %lf, %lf", &templat, &templon, &tempheight);
2120 AddElevation(latitude, longitude, height);
2122 fscanf(fd1,"%lf, %lf, %lf", &latitude, &longitude, &height);
2131 fprintf(stderr,"*** ERROR: \"%s\": not found!\n",filename);
2134 void LoadBoundaries(char *filename)
2136 /* This function reads Cartographic Boundary Files available from
2137 the U.S. Census Bureau, and plots the data contained in those
2138 files on the PPM Map generated by SPLAT!. Such files contain
2139 the coordinates that describe the boundaries of cities,
2140 counties, and states. */
2143 double lat0, lon0, lat1, lon1;
2145 struct site source, destination;
2148 fd=fopen(filename,"r");
2152 fgets(string,78,fd);
2154 fprintf(stdout,"Reading \"%s\"... ",filename);
2159 fgets(string,78,fd);
2160 sscanf(string,"%lf %lf", &lon0, &lat0);
2161 fgets(string,78,fd);
2165 sscanf(string,"%lf %lf", &lon1, &lat1);
2172 destination.lat=lat1;
2173 destination.lon=lon1;
2175 ReadPath(source,destination);
2177 for (x=0; x<path.length; x++)
2178 OrMask(path.lat[x],path.lon[x],4);
2183 fgets(string,78,fd);
2185 } while (strncmp(string,"END",3)!=0 && feof(fd)==0);
2187 fgets(string,78,fd);
2189 } while (strncmp(string,"END",3)!=0 && feof(fd)==0);
2193 fprintf(stdout,"Done!\n");
2198 fprintf(stderr,"*** ERROR: \"%s\": not found!\n",filename);
2201 void ReadLRParm(char *txsite_filename)
2203 /* This function reads Longley-Rice parameter data for the
2204 transmitter site. The file name is the same as the txsite,
2205 except the filename extension is .lrp. If the needed file
2206 is not found, then the file "splat.lrp" is read from the
2207 current working directory. Failure to load this file will
2208 result in the default parameters hard coded into this
2209 function to be used and written to "splat.lrp". */
2212 char filename[255], string[80], *pointer=NULL;
2214 FILE *fd=NULL, *outfile=NULL;
2216 /* Default parameters in case things go bad */
2218 LR.eps_dielect=15.0;
2219 LR.sgm_conductivity=0.005;
2220 LR.eno_ns_surfref=301.0;
2227 /* Modify txsite filename to one with a .lrp extension. */
2229 strncpy(filename,txsite_filename,255);
2231 for (x=0; filename[x]!='.' && filename[x]!=0 && filename[x]!='\n' && x<249; x++);
2239 fd=fopen(filename,"r");
2243 /* Load default "splat.lrp" file */
2245 strncpy(filename,"splat.lrp\0",10);
2246 fd=fopen(filename,"r");
2251 fgets(string,80,fd);
2253 pointer=strchr(string,';');
2258 ok=sscanf(string,"%lf", &din);
2264 fgets(string,80,fd);
2266 pointer=strchr(string,';');
2271 ok=sscanf(string,"%lf", &din);
2276 LR.sgm_conductivity=din;
2278 fgets(string,80,fd);
2280 pointer=strchr(string,';');
2285 ok=sscanf(string,"%lf", &din);
2290 LR.eno_ns_surfref=din;
2292 fgets(string,80,fd);
2294 pointer=strchr(string,';');
2299 ok=sscanf(string,"%lf", &din);
2306 fgets(string,80,fd);
2308 pointer=strchr(string,';');
2313 ok=sscanf(string,"%d", &iin);
2318 LR.radio_climate=iin;
2320 fgets(string,80,fd);
2322 pointer=strchr(string,';');
2327 ok=sscanf(string,"%d", &iin);
2334 fgets(string,80,fd);
2336 pointer=strchr(string,';');
2341 ok=sscanf(string,"%lf", &din);
2348 fgets(string,80,fd);
2350 pointer=strchr(string,';');
2355 ok=sscanf(string,"%lf", &din);
2369 /* Create a "splat.lrp" file since one
2370 could not be successfully loaded. */
2372 outfile=fopen("splat.lrp","w");
2374 fprintf(outfile,"%.3f\t; Earth Dielectric Constant (Relative permittivity)\n",LR.eps_dielect);
2375 fprintf(outfile,"%.3f\t; Earth Conductivity (Siemens per meter)\n", LR.sgm_conductivity);
2376 fprintf(outfile,"%.3f\t; Atmospheric Bending Constant (N-Units)\n",LR.eno_ns_surfref);
2377 fprintf(outfile,"%.3f\t; Frequency in MHz (20 MHz to 20 GHz)\n", LR.frq_mhz);
2378 fprintf(outfile,"%d\t; Radio Climate\n",LR.radio_climate);
2379 fprintf(outfile,"%d\t; Polarization (0 = Horizontal, 1 = Vertical)\n", LR.pol);
2380 fprintf(outfile,"%.2f\t; Fraction of situations\n",LR.conf);
2381 fprintf(outfile, "%.2f\t; Fraction of time\n",LR.rel);
2382 fprintf(outfile,"\nPlease consult SPLAT! documentation for the meaning and use of this data.\n");
2386 fprintf(stderr,"\n%c*** There were problems reading your \"%s\" file! ***\nA \"splat.lrp\" file was written to your directory with default data.\n",7,filename);
2389 if (fd==NULL || ok==0)
2390 fprintf(stderr,"Longley-Rice default parameters have been assumed for this analysis.\n");
2393 void PlotPath(struct site source, struct site destination, char mask_value)
2395 /* This function analyzes the path between the source and
2396 destination locations. It determines which points along
2397 the path have line-of-sight visibility to the source.
2398 Points along with path having line-of-sight visibility
2399 to the source at an AGL altitude equal to that of the
2400 destination location are stored by setting bit 1 in the
2401 mask[][] array, which are displayed in green when PPM
2402 maps are later generated by SPLAT!. */
2406 register double cos_xmtr_angle, cos_test_angle, test_alt;
2407 double distance, rx_alt, tx_alt;
2409 ReadPath(source,destination);
2411 for (y=0; y<path.length; y++)
2413 /* Test this point only if it hasn't been already
2414 tested and found to be free of obstructions. */
2416 if ((GetMask(path.lat[y],path.lon[y])&mask_value)==0)
2418 distance=5280.0*path.distance[y];
2419 tx_alt=earthradius+source.alt+path.elevation[0];
2420 rx_alt=earthradius+destination.alt+path.elevation[y];
2422 /* Calculate the cosine of the elevation of the
2423 transmitter as seen at the temp rx point. */
2425 cos_xmtr_angle=((rx_alt*rx_alt)+(distance*distance)-(tx_alt*tx_alt))/(2.0*rx_alt*distance);
2427 for (x=y, block=0; x>=0 && block==0; x--)
2429 distance=5280.0*(path.distance[y]-path.distance[x]);
2430 test_alt=earthradius+path.elevation[x];
2432 cos_test_angle=((rx_alt*rx_alt)+(distance*distance)-(test_alt*test_alt))/(2.0*rx_alt*distance);
2434 /* Compare these two angles to determine if
2435 an obstruction exists. Since we're comparing
2436 the cosines of these angles rather than
2437 the angles themselves, the following "if"
2438 statement is reversed from what it would
2439 be if the actual angles were compared. */
2441 if (cos_xmtr_angle>cos_test_angle)
2446 OrMask(path.lat[y],path.lon[y],mask_value);
2451 void PlotLRPath(struct site source, struct site destination, FILE *fd)
2453 /* This function plots the RF path loss between source and
2454 destination points based on the Longley-Rice propagation
2455 model, taking into account antenna pattern data, if available. */
2457 char block=0, strmode[100];
2459 double loss, azimuth, pattern=0.0,
2460 source_alt, dest_alt, source_alt2, dest_alt2,
2461 cos_xmtr_angle, cos_test_angle=0.0, test_alt,
2462 elevation, distance=0.0,
2466 ReadPath(source,destination);
2468 four_thirds_earth=EARTHRADIUS*(4.0/3.0);
2470 /* Copy elevations along path into the elev_l[] array. */
2472 for (x=0; x<path.length; x++)
2473 elev_l[x+2]=path.elevation[x]*METERS_PER_FOOT;
2475 /* Since the only energy the Longley-Rice model considers
2476 reaching the destination is based on what is scattered
2477 or deflected from the first obstruction along the path,
2478 we first need to find the location and elevation angle
2479 of that first obstruction (if it exists). This is done
2480 using a 4/3rds Earth radius to match the model used by
2481 Longley-Rice. This information is required for properly
2482 integrating the antenna's elevation pattern into the
2483 calculation for overall path loss. (Using path.length-1
2484 below avoids a Longley-Rice model error from occuring at
2485 the destination point.) */
2487 for (y=2; (y<(path.length-1) && path.distance[y]<=max_range); y++)
2489 /* Process this point only if it
2490 has not already been processed. */
2492 if (GetMask(path.lat[y],path.lon[y])==0)
2494 distance=5280.0*path.distance[y];
2495 source_alt=four_thirds_earth+source.alt+path.elevation[0];
2496 dest_alt=four_thirds_earth+destination.alt+path.elevation[y];
2497 dest_alt2=dest_alt*dest_alt;
2498 source_alt2=source_alt*source_alt;
2500 /* Calculate the cosine of the elevation of
2501 the receiver as seen by the transmitter. */
2503 cos_xmtr_angle=((source_alt2)+(distance*distance)-(dest_alt2))/(2.0*source_alt*distance);
2505 if (got_elevation_pattern || fd!=NULL)
2507 /* If no antenna elevation pattern is available, and
2508 no output file is designated, the following code
2509 that determines the elevation angle to the first
2510 obstruction along the path is bypassed. */
2512 for (x=2, block=0; x<y && block==0; x++)
2514 distance=5280.0*(path.distance[y]-path.distance[x]);
2515 test_alt=four_thirds_earth+path.elevation[x];
2517 /* Calculate the cosine of the elevation
2518 angle of the terrain (test point)
2519 as seen by the transmitter. */
2521 cos_test_angle=((source_alt2)+(distance*distance)-(test_alt*test_alt))/(2.0*source_alt*distance);
2523 /* Compare these two angles to determine if
2524 an obstruction exists. Since we're comparing
2525 the cosines of these angles rather than
2526 the angles themselves, the sense of the
2527 following "if" statement is reversed from
2528 what it would be if the angles themselves
2531 if (cos_xmtr_angle>cos_test_angle)
2535 /* At this point, we have the elevation angle
2536 to the first obstruction (if it exists). */
2539 /* Determine attenuation for each point along the
2540 path using Longley-Rice's point_to_point mode
2541 starting at y=2 (number_of_points = 1), the
2542 shortest distance terrain can play a role in
2545 elev_l[0]=y-1; /* (number of points - 1) */
2547 /* Distance between elevation samples */
2548 elev_l[1]=METERS_PER_MILE*(path.distance[y]-path.distance[y-1]);
2550 point_to_point(elev_l,source.alt*METERS_PER_FOOT,
2551 destination.alt*METERS_PER_FOOT, LR.eps_dielect,
2552 LR.sgm_conductivity, LR.eno_ns_surfref, LR.frq_mhz,
2553 LR.radio_climate, LR.pol, LR.conf, LR.rel, loss,
2557 elevation=((acos(cos_test_angle))/deg2rad)-90.0;
2560 elevation=((acos(cos_xmtr_angle))/deg2rad)-90.0;
2562 temp.lat=path.lat[y];
2563 temp.lon=path.lon[y];
2565 azimuth=(Azimuth(source,temp));
2569 /* Write path loss data to output file */
2571 fprintf(fd,"%.7f, %.7f, %.3f, %.3f, %.2f\n",path.lat[y], path.lon[y], azimuth, elevation, loss);
2574 /* Integrate the antenna's radiation
2575 pattern into the overall path loss. */
2577 x=(int)rint(10.0*(10.0-elevation));
2579 if (x>=0 && x<=1000)
2581 azimuth=rint(azimuth);
2583 pattern=(double)LR.antenna_pattern[(int)azimuth][x];
2587 pattern=20.0*log10(pattern);
2602 OrMask(path.lat[y],path.lon[y],((unsigned char)(loss))<<3);
2605 else if (GetMask(path.lat[y],path.lon[y])==0 && path.distance[y]>max_range)
2606 OrMask(path.lat[y],path.lon[y],1);
2610 void PlotCoverage(struct site source, double altitude)
2612 /* This function performs a 360 degree sweep around the
2613 transmitter site (source location), and plots the
2614 line-of-sight coverage of the transmitter on the SPLAT!
2615 generated topographic map based on a receiver located
2616 at the specified altitude (in feet AGL). Results
2617 are stored in memory, and written out in the form
2618 of a topographic map when the WritePPM() function
2619 is later invoked. */
2621 float lat, lon, one_pixel;
2622 static unsigned char mask_value;
2625 unsigned char symbol[4], x;
2627 /* Initialize mask_value */
2629 if (mask_value!=8 && mask_value!=16 && mask_value!=32)
2632 one_pixel=1.0/1200.0;
2641 fprintf(stdout,"\nComputing line-of-sight coverage of %s with an RX antenna\nat %.2f %s AGL:\n\n 0%c to 25%c ",source.name,metric?altitude*METERS_PER_FOOT:altitude,metric?"meters":"feet",37,37);
2644 /* 18.75=1200 pixels/degree divided by 64 loops
2645 per progress indicator symbol (.oOo) printed. */
2647 z=(int)(18.75*ReduceAngle(max_west-min_west));
2649 for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2658 PlotPath(source,edge,mask_value);
2663 fprintf(stdout,"%c",symbol[x]);
2675 fprintf(stdout,"\n25%c to 50%c ",37,37);
2678 z=(int)(18.75*(max_north-min_north));
2680 for (lat=max_north, x=0; lat>=min_north; lat-=one_pixel)
2686 PlotPath(source,edge,mask_value);
2691 fprintf(stdout,"%c",symbol[x]);
2703 fprintf(stdout,"\n50%c to 75%c ",37,37);
2706 z=(int)(18.75*ReduceAngle(max_west-min_west));
2708 for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2717 PlotPath(source,edge,mask_value);
2722 fprintf(stdout,"%c",symbol[x]);
2734 fprintf(stdout,"\n75%c to 100%c ",37,37);
2737 z=(int)(18.75*(max_north-min_north));
2739 for (lat=min_north, x=0; lat<=max_north; lat+=one_pixel)
2745 PlotPath(source,edge,mask_value);
2750 fprintf(stdout,"%c",symbol[x]);
2761 fprintf(stdout,"\nDone!\n");
2764 /* Assign next mask value */
2781 void PlotLRMap(struct site source, double altitude, char *plo_filename)
2783 /* This function performs a 360 degree sweep around the
2784 transmitter site (source location), and plots the
2785 Longley-Rice attenuation on the SPLAT! generated
2786 topographic map based on a receiver located at
2787 the specified altitude (in feet AGL). Results
2788 are stored in memory, and written out in the form
2789 of a topographic map when the WritePPMLR() function
2790 is later invoked. */
2794 float lat, lon, one_pixel;
2795 unsigned char symbol[4], x;
2798 one_pixel=1.0/1200.0;
2807 fprintf(stdout,"\nComputing Longley-Rice coverage of %s ", source.name);
2809 fprintf(stdout,"out to a radius\nof %.2f %s with an RX antenna at %.2f %s AGL:\n\n 0%c to 25%c ",metric?max_range*KM_PER_MILE:max_range,metric?"kilometers":"miles",metric?altitude*METERS_PER_FOOT:altitude,metric?"meters":"feet",37,37);
2812 if (plo_filename[0]!=0)
2813 fd=fopen(plo_filename,"wb");
2817 /* Write header information to output file */
2819 fprintf(fd,"%d, %d\t; max_west, min_west\n%d, %d\t; max_north, min_north\n",max_west, min_west, max_north, min_north);
2822 /* 18.75=1200 pixels/degree divided by 64 loops
2823 per progress indicator symbol (.oOo) printed. */
2825 z=(int)(18.75*ReduceAngle(max_west-min_west));
2827 for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2836 PlotLRPath(source,edge,fd);
2841 fprintf(stdout,"%c",symbol[x]);
2853 fprintf(stdout,"\n25%c to 50%c ",37,37);
2856 z=(int)(18.75*(max_north-min_north));
2858 for (lat=max_north, x=0; lat>=min_north; lat-=one_pixel)
2864 PlotLRPath(source,edge,fd);
2869 fprintf(stdout,"%c",symbol[x]);
2881 fprintf(stdout,"\n50%c to 75%c ",37,37);
2884 z=(int)(18.75*ReduceAngle(max_west-min_west));
2886 for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2895 PlotLRPath(source,edge,fd);
2900 fprintf(stdout,"%c",symbol[x]);
2912 fprintf(stdout,"\n75%c to 100%c ",37,37);
2915 z=(int)(18.75*(max_north-min_north));
2917 for (lat=min_north, x=0; lat<=max_north; lat+=one_pixel)
2923 PlotLRPath(source,edge,fd);
2928 fprintf(stdout,"%c",symbol[x]);
2942 fprintf(stdout,"\nDone!\n");
2946 void WritePPM(char *filename, unsigned char geo)
2948 /* This function generates a topographic map in Portable Pix Map
2949 (PPM) format based on logarithmically scaled topology data,
2950 as well as the content of flags held in the mask[][] array.
2951 The image created is rotated counter-clockwise 90 degrees
2952 from its representation in dem[][] so that north points
2953 up and east points right in the image generated. */
2955 char mapfile[255], geofile[255];
2956 unsigned char found, mask;
2957 unsigned width, height, output;
2958 int indx, x, y, x0=0, y0=0, minlat, minlon;
2959 float lat, lon, one_pixel, conversion, one_over_gamma;
2962 one_pixel=1.0/1200.0;
2963 one_over_gamma=1.0/GAMMA;
2964 conversion=255.0/pow((double)(max_elevation-min_elevation),one_over_gamma);
2966 width=(unsigned)(1200*ReduceAngle(max_west-min_west));
2967 height=(unsigned)(1200*ReduceAngle(max_north-min_north));
2970 strncpy(mapfile, "map.ppm\0",8);
2973 for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
2975 mapfile[x]=filename[x];
2976 geofile[x]=filename[x];
2993 fd=fopen(geofile,"wb");
2995 fprintf(fd,"FILENAME\t%s\n",mapfile);
2996 fprintf(fd,"#\t\tX\tY\tLong\t\tLat\n");
2997 fprintf(fd,"TIEPOINT\t0\t0\t%d.000\t\t%d.000\n",(max_west<180?-max_west:360-max_west),max_north);
2998 fprintf(fd,"TIEPOINT\t%u\t%u\t%d.000\t\t%d.000\n",width-1,height-1,(min_west<180?-min_west:360-min_west),min_north);
2999 fprintf(fd,"IMAGESIZE\t%u\t%u\n",width,height);
3000 fprintf(fd,"#\n# Auto Generated by SPLAT! v%s\n#\n",splat_version);
3005 fd=fopen(mapfile,"wb");
3007 fprintf(fd,"P6\n%u %u\n255\n",width,height);
3009 fprintf(stdout,"\nWriting \"%s\" (%ux%u pixmap image)... ",mapfile,width,height);
3012 for (y=0, lat=((double)max_north)-one_pixel; y<(int)height; y++, lat-=one_pixel)
3014 minlat=(int)floor(lat);
3016 for (x=0, lon=((double)max_west)-one_pixel; x<(int)width; x++, lon-=one_pixel)
3021 minlon=(int)floor(lon);
3023 for (indx=0, found=0; indx<MAXSLOTS && found==0;)
3024 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
3031 x0=(int)(1199.0*(lat-floor(lat)));
3032 y0=(int)(1199.0*(lon-floor(lon)));
3034 mask=dem[indx].mask[x0][y0];
3037 /* Text Labels: Red */
3038 fprintf(fd,"%c%c%c",255,0,0);
3041 /* County Boundaries: Light Cyan */
3042 fprintf(fd,"%c%c%c",128,128,255);
3044 else switch (mask&57)
3048 fprintf(fd,"%c%c%c",0,255,0);
3053 fprintf(fd,"%c%c%c",0,255,255);
3057 /* TX1 + TX2: Yellow */
3058 fprintf(fd,"%c%c%c",255,255,0);
3062 /* TX3: Medium Violet */
3063 fprintf(fd,"%c%c%c",147,112,219);
3067 /* TX1 + TX3: Pink */
3068 fprintf(fd,"%c%c%c",255,192,203);
3072 /* TX2 + TX3: Orange */
3073 fprintf(fd,"%c%c%c",255,165,0);
3077 /* TX1 + TX2 + TX3: Dark Green */
3078 fprintf(fd,"%c%c%c",0,100,0);
3083 fprintf(fd,"%c%c%c",255,130,71);
3087 /* TX1 + TX4: Green Yellow */
3088 fprintf(fd,"%c%c%c",173,255,47);
3092 /* TX2 + TX4: Dark Sea Green 1 */
3093 fprintf(fd,"%c%c%c",193,255,193);
3097 /* TX1 + TX2 + TX4: Blanched Almond */
3098 fprintf(fd,"%c%c%c",255,235,205);
3102 /* TX3 + TX4: Dark Turquoise */
3103 fprintf(fd,"%c%c%c",0,206,209);
3107 /* TX1 + TX3 + TX4: Medium Spring Green */
3108 fprintf(fd,"%c%c%c",0,250,154);
3112 /* TX2 + TX3 + TX4: Tan */
3113 fprintf(fd,"%c%c%c",210,180,140);
3117 /* TX1 + TX2 + TX3 + TX4: Gold2 */
3118 fprintf(fd,"%c%c%c",238,201,0);
3122 /* Water: Medium Blue */
3123 if (dem[indx].data[x0][y0]==0)
3124 fprintf(fd,"%c%c%c",0,0,170);
3127 /* Elevation: Greyscale */
3128 output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
3129 fprintf(fd,"%c%c%c",output,output,output);
3136 /* We should never get here, but if */
3137 /* we do, display the region as black */
3139 fprintf(fd,"%c%c%c",0,0,0);
3145 fprintf(stdout,"Done!\n");
3149 void WritePPMLR(char *filename, unsigned char geo)
3151 /* This function generates a topographic map in Portable Pix Map
3152 (PPM) format based on the content of flags held in the mask[][]
3153 array (only). The image created is rotated counter-clockwise
3154 90 degrees from its representation in dem[][] so that north
3155 points up and east points right in the image generated. */
3157 char mapfile[255], geofile[255];
3158 unsigned width, height, output;
3159 unsigned char found, mask, cityorcounty;
3160 int indx, x, y, t, t2, x0, y0, minlat, minlon, loss;
3161 float lat, lon, one_pixel, conversion, one_over_gamma;
3164 one_pixel=1.0/1200.0;
3165 one_over_gamma=1.0/GAMMA;
3166 conversion=255.0/pow((double)(max_elevation-min_elevation),one_over_gamma);
3168 width=(unsigned)(1200*ReduceAngle(max_west-min_west));
3169 height=(unsigned)(1200*ReduceAngle(max_north-min_north));
3172 strncpy(mapfile, "map.ppm\0",8);
3175 for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
3177 mapfile[x]=filename[x];
3178 geofile[x]=filename[x];
3195 fd=fopen(geofile,"wb");
3197 fprintf(fd,"FILENAME\t%s\n",mapfile);
3198 fprintf(fd,"#\t\tX\tY\tLong\t\tLat\n");
3199 fprintf(fd,"TIEPOINT\t0\t0\t%d.000\t\t%d.000\n",(max_west<180?-max_west:360-max_west),max_north);
3200 fprintf(fd,"TIEPOINT\t%u\t%u\t%d.000\t\t%.3f\n",width-1,height+29,(min_west<180?-min_west:360-min_west),(double)(min_north-0.025));
3201 fprintf(fd,"IMAGESIZE\t%u\t%u\n",width,height+30);
3202 fprintf(fd,"#\n# Auto Generated by SPLAT! v%s\n#\n",splat_version);
3207 fd=fopen(mapfile,"wb");
3209 fprintf(fd,"P6\n%u %u\n255\n",width,height+30);
3211 fprintf(stdout,"\nWriting \"%s\" (%ux%u pixmap image)... ",mapfile,width,height+30);
3214 for (y=0, lat=((double)max_north)-one_pixel; y<(int)height; y++, lat-=one_pixel)
3216 minlat=(int)floor(lat);
3218 for (x=0, lon=((double)max_west)-one_pixel; x<(int)width; x++, lon-=one_pixel)
3223 minlon=(int)floor(lon);
3225 for (indx=0, found=0; indx<MAXSLOTS && found==0;)
3226 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
3232 x0=(int)(1199.0*(lat-floor(lat)));
3233 y0=(int)(1199.0*(lon-floor(lon)));
3235 mask=dem[indx].mask[x0][y0];
3236 loss=70+(10*(int)((mask&248)>>3));
3241 /* Text Labels - Black or Red */
3243 if ((mask&120) && (loss<=90))
3244 fprintf(fd,"%c%c%c",0,0,0);
3246 fprintf(fd,"%c%c%c",255,0,0);
3253 /* County Boundaries: Black */
3255 fprintf(fd,"%c%c%c",0,0,0);
3260 if (cityorcounty==0)
3264 { /* Display land or sea elevation */
3266 if (dem[indx].data[x0][y0]==0)
3267 fprintf(fd,"%c%c%c",0,0,170);
3270 /* Elevation: Greyscale */
3271 output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
3272 fprintf(fd,"%c%c%c",output,output,output);
3278 /* Plot signal loss in color */
3281 fprintf(fd,"%c%c%c",255,0,0);
3285 fprintf(fd,"%c%c%c",255,128,0);
3289 fprintf(fd,"%c%c%c",255,165,0);
3293 fprintf(fd,"%c%c%c",255,206,0);
3297 fprintf(fd,"%c%c%c",255,255,0);
3301 fprintf(fd,"%c%c%c",184,255,0);
3305 fprintf(fd,"%c%c%c",0,255,0);
3309 fprintf(fd,"%c%c%c",0,208,0);
3313 fprintf(fd,"%c%c%c",0,196,196);
3317 fprintf(fd,"%c%c%c",0,148,255);
3321 fprintf(fd,"%c%c%c",80,80,255);
3325 fprintf(fd,"%c%c%c",0,38,255);
3329 fprintf(fd,"%c%c%c",142,63,255);
3333 fprintf(fd,"%c%c%c",196,54,255);
3337 fprintf(fd,"%c%c%c",255,0,255);
3341 fprintf(fd,"%c%c%c",255,194,204);
3346 if (dem[indx].data[x0][y0]==0)
3347 fprintf(fd,"%c%c%c",0,0,170);
3350 /* Elevation: Greyscale */
3351 output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
3352 fprintf(fd,"%c%c%c",output,output,output);
3360 /* We should never get here, but if */
3361 /* we do, display the region as black */
3363 fprintf(fd,"%c%c%c",0,0,0);
3368 /* Display legend along bottom of image */
3372 for (y0=0; y0<30; y0++)
3374 for (indx=0; indx<16; indx++)
3376 for (x=0; x<x0; x++)
3381 if (y0>=10 && y0<=18)
3386 if (smallfont[t2/10][y0-10][x-11])
3391 if (smallfont[t2%10][y0-10][x-19])
3395 if (smallfont[0][y0-10][x-27])
3402 fprintf(fd,"%c%c%c",255,0,0);
3406 fprintf(fd,"%c%c%c",255,128,0);
3410 fprintf(fd,"%c%c%c",255,165,0);
3414 fprintf(fd,"%c%c%c",255,206,0);
3418 fprintf(fd,"%c%c%c",255,255,0);
3422 fprintf(fd,"%c%c%c",184,255,0);
3426 fprintf(fd,"%c%c%c",0,255,0);
3430 fprintf(fd,"%c%c%c",0,208,0);
3434 fprintf(fd,"%c%c%c",0,196,196);
3438 fprintf(fd,"%c%c%c",0,148,255);
3442 fprintf(fd,"%c%c%c",80,80,255);
3446 fprintf(fd,"%c%c%c",0,38,255);
3450 fprintf(fd,"%c%c%c",142,63,255);
3454 fprintf(fd,"%c%c%c",196,54,255);
3458 fprintf(fd,"%c%c%c",255,0,255);
3463 fprintf(fd,"%c%c%c",0,0,0);
3467 fprintf(fd,"%c%c%c",255,194,204);
3474 fprintf(stdout,"Done!\n");
3478 void GraphTerrain(struct site source, struct site destination, char *name)
3480 /* This function invokes gnuplot to generate an appropriate
3481 output file indicating the terrain profile between the source
3482 and destination locations. "filename" is the name assigned
3483 to the output file generated by gnuplot. The filename extension
3484 is used to set gnuplot's terminal setting and output file type.
3485 If no extension is found, .png is assumed. */
3488 char filename[255], term[30], ext[15];
3491 ReadPath(destination,source);
3493 fd=fopen("profile.gp","wb");
3495 for (x=0; x<path.length; x++)
3498 fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*path.elevation[x]);
3501 fprintf(fd,"%f\t%f\n",path.distance[x],path.elevation[x]);
3508 /* Default filename and output file type */
3510 strncpy(filename,"profile\0",8);
3511 strncpy(term,"png\0",4);
3512 strncpy(ext,"png\0",4);
3517 /* Grab extension and terminal type from "name" */
3519 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
3520 filename[x]=name[x];
3524 for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
3526 term[y]=tolower(name[x]);
3536 { /* No extension -- Default is png */
3539 strncpy(term,"png\0",4);
3540 strncpy(ext,"png\0",4);
3544 /* Either .ps or .postscript may be used
3545 as an extension for postscript output. */
3547 if (strncmp(term,"postscript",10)==0)
3548 strncpy(ext,"ps\0",3);
3550 else if (strncmp(ext,"ps",2)==0)
3551 strncpy(term,"postscript enhanced color\0",26);
3553 fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
3556 fd=fopen("splat.gp","w");
3557 fprintf(fd,"set grid\n");
3558 fprintf(fd,"set autoscale\n");
3559 fprintf(fd,"set encoding iso_8859_1\n");
3560 fprintf(fd,"set term %s\n",term);
3561 fprintf(fd,"set title \"SPLAT! Terrain Profile Between %s and %s (%.2f%c Azimuth)\"\n",destination.name, source.name, Azimuth(destination,source),176);
3565 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(source,destination));
3566 fprintf(fd,"set ylabel \"Ground Elevation Above Sea Level (meters)\"\n");
3573 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(source,destination));
3574 fprintf(fd,"set ylabel \"Ground Elevation Above Sea Level (feet)\"\n");
3577 fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
3578 fprintf(fd,"plot \"profile.gp\" title \"\" with lines\n");
3581 x=system("gnuplot splat.gp");
3586 unlink("profile.gp");
3587 fprintf(stdout," Done!\n");
3592 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
3595 void GraphElevation(struct site source, struct site destination, char *name)
3597 /* This function invokes gnuplot to generate an appropriate
3598 output file indicating the terrain profile between the source
3599 and destination locations. "filename" is the name assigned
3600 to the output file generated by gnuplot. The filename extension
3601 is used to set gnuplot's terminal setting and output file type.
3602 If no extension is found, .png is assumed. */
3605 char filename[255], term[30], ext[15];
3606 double angle, refangle, maxangle=-90.0;
3608 FILE *fd=NULL, *fd2=NULL;
3610 ReadPath(destination,source); /* destination=RX, source=TX */
3611 refangle=ElevationAngle(destination,source);
3613 fd=fopen("profile.gp","wb");
3614 fd2=fopen("reference.gp","wb");
3616 for (x=1; x<path.length-1; x++)
3618 remote.lat=path.lat[x];
3619 remote.lon=path.lon[x];
3621 angle=ElevationAngle(destination,remote);
3625 fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[x],angle);
3626 fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[x],refangle);
3631 fprintf(fd,"%f\t%f\n",path.distance[x],angle);
3632 fprintf(fd2,"%f\t%f\n",path.distance[x],refangle);
3641 fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],refangle);
3642 fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],refangle);
3647 fprintf(fd,"%f\t%f\n",path.distance[path.length-1],refangle);
3648 fprintf(fd2,"%f\t%f\n",path.distance[path.length-1],refangle);
3656 /* Default filename and output file type */
3658 strncpy(filename,"profile\0",8);
3659 strncpy(term,"png\0",4);
3660 strncpy(ext,"png\0",4);
3665 /* Grab extension and terminal type from "name" */
3667 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
3668 filename[x]=name[x];
3672 for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
3674 term[y]=tolower(name[x]);
3684 { /* No extension -- Default is png */
3687 strncpy(term,"png\0",4);
3688 strncpy(ext,"png\0",4);
3692 /* Either .ps or .postscript may be used
3693 as an extension for postscript output. */
3695 if (strncmp(term,"postscript",10)==0)
3696 strncpy(ext,"ps\0",3);
3698 else if (strncmp(ext,"ps",2)==0)
3699 strncpy(term,"postscript enhanced color\0",26);
3701 fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
3704 fd=fopen("splat.gp","w");
3706 fprintf(fd,"set grid\n");
3707 fprintf(fd,"set yrange [%2.3f to %2.3f]\n", (-fabs(refangle)-0.25), maxangle+0.25);
3708 fprintf(fd,"set encoding iso_8859_1\n");
3709 fprintf(fd,"set term %s\n",term);
3710 fprintf(fd,"set title \"SPLAT! Elevation Profile Between %s and %s (%.2f%c azimuth)\"\n",destination.name,source.name,Azimuth(destination,source),176);
3713 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(source,destination));
3715 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(source,destination));
3718 fprintf(fd,"set ylabel \"Elevation Angle Along LOS Path Between %s and %s (degrees)\"\n",destination.name,source.name);
3719 fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
3720 fprintf(fd,"plot \"profile.gp\" title \"Real Earth Profile\" with lines, \"reference.gp\" title \"Line of Sight Path (%.2f%c elevation)\" with lines\n",refangle,176);
3724 x=system("gnuplot splat.gp");
3729 unlink("profile.gp");
3730 unlink("reference.gp");
3732 fprintf(stdout," Done!\n");
3737 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
3740 void GraphHeight(struct site source, struct site destination, char *name, double f, unsigned char n)
3742 /* This function invokes gnuplot to generate an appropriate
3743 output file indicating the terrain profile between the source
3744 and destination locations referenced to the line-of-sight path
3745 between the receive and transmit sites. "filename" is the name
3746 assigned to the output file generated by gnuplot. The filename
3747 extension is used to set gnuplot's terminal setting and output
3748 file type. If no extension is found, .png is assumed. */
3751 char filename[255], term[30], ext[15];
3752 double a, b, c, height=0.0, refangle, cangle, maxheight=-100000.0,
3753 minheight=100000.0, lambda=0.0, f_zone=0.0, fpt6_zone=0.0,
3754 nm=0.0, nb=0.0, ed=0.0, es=0.0, r=0.0, d=0.0, d1=0.0,
3755 terrain, azimuth, distance, dheight=0.0, minterrain=100000.0,
3756 minearth=100000.0, miny, maxy, min2y, max2y;
3758 FILE *fd=NULL, *fd2=NULL, *fd3=NULL, *fd4=NULL, *fd5=NULL;
3760 ReadPath(destination,source); /* destination=RX, source=TX */
3761 azimuth=Azimuth(destination,source);
3762 distance=Distance(destination,source);
3763 refangle=ElevationAngle(destination,source);
3764 b=GetElevation(destination)+destination.alt+earthradius;
3766 /* Wavelength and path distance (great circle) in feet. */
3770 lambda=9.8425e8/(f*1e6);
3771 d=5280.0*path.distance[path.length-1];
3776 ed=GetElevation(destination);
3777 es=GetElevation(source);
3778 nb=-destination.alt-ed;
3779 nm=(-source.alt-es-nb)/(path.distance[path.length-1]);
3782 fd=fopen("profile.gp","wb");
3783 fd2=fopen("reference.gp","wb");
3784 fd5=fopen("curvature.gp", "wb");
3788 fd3=fopen("fresnel.gp", "wb");
3789 fd4=fopen("fresnel_pt_6.gp", "wb");
3792 for (x=0; x<path.length; x++)
3794 remote.lat=path.lat[x];
3795 remote.lon=path.lon[x];
3798 terrain=GetElevation(remote);
3801 terrain+=destination.alt; /* RX antenna spike */
3803 a=terrain+earthradius;
3804 cangle=5280.0*Distance(destination,remote)/earthradius;
3805 c=b*sin(refangle*deg2rad+HALFPI)/sin(HALFPI-refangle*deg2rad-cangle);
3809 /* Per Fink and Christiansen, Electronics
3810 * Engineers' Handbook, 1989:
3812 * H = sqrt(lamba * d1 * (d - d1)/d)
3814 * where H is the distance from the LOS
3815 * path to the first Fresnel zone boundary.
3820 d1=5280.0*path.distance[x];
3821 f_zone=-1*sqrt(lambda*d1*(d-d1)/d);
3822 fpt6_zone=0.6*f_zone;
3827 r=-(nm*path.distance[x])-nb;
3842 fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*height);
3843 fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*r);
3844 fprintf(fd5,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*(height-terrain));
3849 fprintf(fd,"%f\t%f\n",path.distance[x],height);
3850 fprintf(fd2,"%f\t%f\n",path.distance[x],r);
3851 fprintf(fd5,"%f\t%f\n",path.distance[x],height-terrain);
3858 fprintf(fd3,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*f_zone);
3859 fprintf(fd4,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*fpt6_zone);
3864 fprintf(fd3,"%f\t%f\n",path.distance[x],f_zone);
3865 fprintf(fd4,"%f\t%f\n",path.distance[x],fpt6_zone);
3868 if (f_zone<minheight)
3872 if (height>maxheight)
3875 if (height<minheight)
3881 if (terrain<minterrain)
3884 if ((height-terrain)<minearth)
3885 minearth=height-terrain;
3889 r=-(nm*path.distance[path.length-1])-nb;
3895 fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3896 fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3901 fprintf(fd,"%f\t%f\n",path.distance[path.length-1],r);
3902 fprintf(fd2,"%f\t%f\n",path.distance[path.length-1],r);
3909 fprintf(fd3,"%f\t%f\n",path.distance[path.length-1],r);
3910 fprintf(fd4,"%f\t%f\n",path.distance[path.length-1],r);
3915 fprintf(fd3,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3916 fprintf(fd4,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3938 /* Default filename and output file type */
3940 strncpy(filename,"height\0",8);
3941 strncpy(term,"png\0",4);
3942 strncpy(ext,"png\0",4);
3947 /* Grab extension and terminal type from "name" */
3949 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
3950 filename[x]=name[x];
3954 for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
3956 term[y]=tolower(name[x]);
3966 { /* No extension -- Default is png */
3969 strncpy(term,"png\0",4);
3970 strncpy(ext,"png\0",4);
3974 /* Either .ps or .postscript may be used
3975 as an extension for postscript output. */
3977 if (strncmp(term,"postscript",10)==0)
3978 strncpy(ext,"ps\0",3);
3980 else if (strncmp(ext,"ps",2)==0)
3981 strncpy(term,"postscript enhanced color\0",26);
3983 fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
3986 fd=fopen("splat.gp","w");
3988 dheight=maxheight-minheight;
3989 miny=minheight-0.15*dheight;
3990 maxy=maxheight+0.05*dheight;
3995 dheight=maxheight-minheight;
3996 min2y=miny-minterrain+0.05*dheight;
4000 miny-=min2y-minearth+0.05*dheight;
4001 min2y=minearth-0.05*dheight;
4004 max2y=min2y+maxy-miny;
4006 fprintf(fd,"set grid\n");
4007 fprintf(fd,"set yrange [%2.3f to %2.3f]\n", metric?miny*METERS_PER_FOOT:miny, metric?maxy*METERS_PER_FOOT:maxy);
4008 fprintf(fd,"set y2range [%2.3f to %2.3f]\n", metric?min2y*METERS_PER_FOOT:min2y, metric?max2y*METERS_PER_FOOT:max2y);
4009 fprintf(fd,"set xrange [-0.5 to %2.3f]\n",metric?KM_PER_MILE*rint(distance+0.5):rint(distance+0.5));
4010 fprintf(fd,"set encoding iso_8859_1\n");
4011 fprintf(fd,"set term %s\n",term);
4014 fprintf(fd,"set title \"SPLAT! Path Profile Between %s and %s (%.2f%c azimuth)\\nWith First Fresnel Zone\"\n",destination.name, source.name, azimuth,176);
4017 fprintf(fd,"set title \"SPLAT! Height Profile Between %s and %s (%.2f%c azimuth)\"\n",destination.name, source.name, azimuth,176);
4020 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(source,destination));
4022 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(source,destination));
4027 fprintf(fd,"set ylabel \"Normalized Height Referenced To LOS Path Between\\n%s and %s (meters)\"\n",destination.name,source.name);
4030 fprintf(fd,"set ylabel \"Normalized Height Referenced To LOS Path Between\\n%s and %s (feet)\"\n",destination.name,source.name);
4037 fprintf(fd,"set ylabel \"Height Referenced To LOS Path Between %s and %s (meters)\"\n",destination.name,source.name);
4040 fprintf(fd,"set ylabel \"Height Referenced To LOS Path Between %s and %s (feet)\"\n",destination.name,source.name);
4043 fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
4046 fprintf(fd,"plot \"profile.gp\" title \"Point-to-Point Profile\" with lines, \"reference.gp\" title \"Line of Sight Path\" with lines, \"curvature.gp\" axes x1y2 title \"Earth's Curvature Contour\" with lines, \"fresnel.gp\" axes x1y1 title \"First Fresnel Zone (%.3f MHz)\" with lines, \"fresnel_pt_6.gp\" title \"60%% of First Fresnel Zone\" with lines\n",f);
4049 fprintf(fd,"plot \"profile.gp\" title \"Point-to-Point Profile\" with lines, \"reference.gp\" title \"Line Of Sight Path\" with lines, \"curvature.gp\" axes x1y2 title \"Earth's Curvature Contour\" with lines\n");
4053 x=system("gnuplot splat.gp");
4058 unlink("profile.gp");
4059 unlink("reference.gp");
4060 unlink("curvature.gp");
4064 unlink("fresnel.gp");
4065 unlink("fresnel_pt_6.gp");
4068 fprintf(stdout," Done!\n");
4073 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
4076 void GraphLongley(struct site source, struct site destination, char *name)
4078 /* This function invokes gnuplot to generate an appropriate
4079 output file indicating the Longley-Rice model loss between
4080 the source and destination locations. "filename" is the
4081 name assigned to the output file generated by gnuplot.
4082 The filename extension is used to set gnuplot's terminal
4083 setting and output file type. If no extension is found,
4086 int x, y, z, errnum, errflag=0;
4087 char filename[255], term[30], ext[15], strmode[100],
4088 report_name[80], block=0;
4089 double maxloss=-100000.0, minloss=100000.0, loss, haavt,
4090 angle1, angle2, azimuth, pattern=1.0, patterndB=0.0,
4091 total_loss=0.0, cos_xmtr_angle, cos_test_angle=0.0,
4092 source_alt, test_alt, dest_alt, source_alt2, dest_alt2,
4093 distance, elevation, four_thirds_earth;
4094 FILE *fd=NULL, *fd2=NULL;
4096 sprintf(report_name,"%s-to-%s.lro",source.name,destination.name);
4098 four_thirds_earth=EARTHRADIUS*(4.0/3.0);
4100 for (x=0; report_name[x]!=0; x++)
4101 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
4104 fd2=fopen(report_name,"w");
4106 fprintf(fd2,"\n\t--==[ SPLAT! v%s Longley-Rice Model Path Loss Report ]==--\n\n",splat_version);
4107 fprintf(fd2,"Analysis of RF path conditions between %s and %s:\n",source.name, destination.name);
4108 fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4109 fprintf(fd2,"Transmitter site: %s\n",source.name);
4111 if (source.lat>=0.0)
4113 fprintf(fd2,"Site location: %.4f North / %.4f West",source.lat, source.lon);
4114 fprintf(fd2, " (%s N / ", dec2dms(source.lat));
4120 fprintf(fd2,"Site location: %.4f South / %.4f West",-source.lat, source.lon);
4121 fprintf(fd2, " (%s S / ", dec2dms(source.lat));
4124 fprintf(fd2, "%s W)\n", dec2dms(source.lon));
4128 fprintf(fd2,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(source));
4129 fprintf(fd2,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*source.alt,METERS_PER_FOOT*(source.alt+GetElevation(source)));
4134 fprintf(fd2,"Ground elevation: %.2f feet AMSL\n",GetElevation(source));
4135 fprintf(fd2,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",source.alt, source.alt+GetElevation(source));
4143 fprintf(fd2,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4145 fprintf(fd2,"Antenna height above average terrain: %.2f feet\n",haavt);
4148 azimuth=Azimuth(source,destination);
4149 angle1=ElevationAngle(source,destination);
4150 angle2=ElevationAngle2(source,destination,earthradius);
4152 if (got_azimuth_pattern || got_elevation_pattern)
4154 x=(int)rint(10.0*(10.0-angle2));
4156 if (x>=0 && x<=1000)
4157 pattern=(double)LR.antenna_pattern[(int)rint(azimuth)][x];
4159 patterndB=20.0*log10(pattern);
4161 fprintf(fd2,"Antenna pattern between %s and %s: %.3f (%.2f dB)\n",source.name, destination.name, pattern, patterndB);
4166 fprintf(fd2,"Distance to %s: %.2f kilometers\n",destination.name,METERS_PER_FOOT*Distance(source,destination));
4169 fprintf(fd2,"Distance to %s: %.2f miles.\n",destination.name,Distance(source,destination));
4171 fprintf(fd2,"Azimuth to %s: %.2f degrees\n",destination.name,azimuth);
4174 fprintf(fd2,"Elevation angle to %s: %+.4f degrees\n",destination.name,angle1);
4177 fprintf(fd2,"Depression angle to %s: %+.4f degrees\n",destination.name,angle1);
4182 fprintf(fd2,"Depression");
4184 fprintf(fd2,"Elevation");
4186 fprintf(fd2," angle to the first obstruction: %+.4f degrees\n",angle2);
4189 fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4193 fprintf(fd2,"Receiver site: %s\n",destination.name);
4195 if (destination.lat>=0.0)
4197 fprintf(fd2,"Site location: %.4f North / %.4f West",destination.lat, destination.lon);
4198 fprintf(fd2, " (%s N / ", dec2dms(destination.lat));
4203 fprintf(fd2,"Site location: %.4f South / %.4f West",-destination.lat, destination.lon);
4204 fprintf(fd2, " (%s S / ", dec2dms(destination.lat));
4207 fprintf(fd2, "%s W)\n", dec2dms(destination.lon));
4211 fprintf(fd2,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(destination));
4212 fprintf(fd2,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*destination.alt, METERS_PER_FOOT*(destination.alt+GetElevation(destination)));
4217 fprintf(fd2,"Ground elevation: %.2f feet AMSL\n",GetElevation(destination));
4218 fprintf(fd2,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",destination.alt, destination.alt+GetElevation(destination));
4221 haavt=haat(destination);
4226 fprintf(fd2,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4228 fprintf(fd2,"Antenna height above average terrain: %.2f feet\n",haavt);
4232 fprintf(fd2,"Distance to %s: %.2f kilometers\n",source.name,KM_PER_MILE*Distance(source,destination));
4235 fprintf(fd2,"Distance to %s: %.2f miles\n",source.name,Distance(source,destination));
4237 azimuth=Azimuth(destination,source);
4239 angle1=ElevationAngle(destination,source);
4240 angle2=ElevationAngle2(destination,source,earthradius);
4242 fprintf(fd2,"Azimuth to %s: %.2f degrees.\n",source.name,azimuth);
4245 fprintf(fd2,"Elevation angle to %s: %+.4f degrees\n",source.name,angle1);
4248 fprintf(fd2,"Depression angle to %s: %+.4f degrees\n",source.name,angle1);
4253 fprintf(fd2,"Depression");
4255 fprintf(fd2,"Elevation");
4257 fprintf(fd2," angle to the first obstruction: %+.4f degrees\n",angle2);
4260 fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4262 fprintf(fd2,"Longley-Rice path calculation parameters used in this analysis:\n\n");
4263 fprintf(fd2,"Earth's Dielectric Constant: %.3lf\n",LR.eps_dielect);
4264 fprintf(fd2,"Earth's Conductivity: %.3lf\n",LR.sgm_conductivity);
4265 fprintf(fd2,"Atmospheric Bending Constant (N): %.3lf\n",LR.eno_ns_surfref);
4266 fprintf(fd2,"Frequency: %.3lf (MHz)\n",LR.frq_mhz);
4267 fprintf(fd2,"Radio Climate: %d (",LR.radio_climate);
4269 switch (LR.radio_climate)
4272 fprintf(fd2,"Equatorial");
4276 fprintf(fd2,"Continental Subtropical");
4280 fprintf(fd2,"Maritime Subtropical");
4284 fprintf(fd2,"Desert");
4288 fprintf(fd2,"Continental Temperate");
4292 fprintf(fd2,"Martitime Temperate, Over Land");
4296 fprintf(fd2,"Maritime Temperate, Over Sea");
4300 fprintf(fd2,"Unknown");
4303 fprintf(fd2,")\nPolarization: %d (",LR.pol);
4306 fprintf(fd2,"Horizontal");
4309 fprintf(fd2,"Vertical");
4311 fprintf(fd2,")\nFraction of Situations: %.1lf%c\n",LR.conf*100.0,37);
4312 fprintf(fd2,"Fraction of Time: %.1lf%c\n",LR.rel*100.0,37);
4313 fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4315 fprintf(fd2,"Analysis Results:\n\n");
4317 ReadPath(source, destination); /* source=TX, destination=RX */
4319 /* Copy elevations along path into the elev_l[] array. */
4321 for (x=0; x<path.length; x++)
4322 elev_l[x+2]=path.elevation[x]*METERS_PER_FOOT;
4325 fprintf(fd2,"Distance (km)");
4327 fprintf(fd2,"Distance (mi)");
4329 fprintf(fd2,"\tLoss (dB)\tErrnum\tComment\n\n");
4331 fd=fopen("profile.gp","w");
4333 azimuth=rint(Azimuth(source,destination));
4335 for (y=2; y<(path.length-1); y++) /* path.length-1 avoids LR error */
4337 distance=5280.0*path.distance[y];
4338 source_alt=four_thirds_earth+source.alt+path.elevation[0];
4339 dest_alt=four_thirds_earth+destination.alt+path.elevation[y];
4340 dest_alt2=dest_alt*dest_alt;
4341 source_alt2=source_alt*source_alt;
4343 /* Calculate the cosine of the elevation of
4344 the receiver as seen by the transmitter. */
4346 cos_xmtr_angle=((source_alt2)+(distance*distance)-(dest_alt2))/(2.0*source_alt*distance);
4348 if (got_elevation_pattern)
4350 /* If an antenna elevation pattern is available, the
4351 following code determines the elevation angle to
4352 the first obstruction along the path. */
4354 for (x=2, block=0; x<y && block==0; x++)
4356 distance=5280.0*(path.distance[y]-path.distance[x]);
4357 test_alt=four_thirds_earth+path.elevation[x];
4359 /* Calculate the cosine of the elevation
4360 angle of the terrain (test point)
4361 as seen by the transmitter. */
4363 cos_test_angle=((source_alt2)+(distance*distance)-(test_alt*test_alt))/(2.0*source_alt*distance);
4365 /* Compare these two angles to determine if
4366 an obstruction exists. Since we're comparing
4367 the cosines of these angles rather than
4368 the angles themselves, the sense of the
4369 following "if" statement is reversed from
4370 what it would be if the angles themselves
4373 if (cos_xmtr_angle>cos_test_angle)
4377 /* At this point, we have the elevation angle
4378 to the first obstruction (if it exists). */
4381 /* Determine path loss for each point along the
4382 path using Longley-Rice's point_to_point mode
4383 starting at x=2 (number_of_points = 1), the
4384 shortest distance terrain can play a role in
4387 elev_l[0]=y-1; /* (number of points - 1) */
4389 /* Distance between elevation samples */
4390 elev_l[1]=METERS_PER_MILE*(path.distance[y]-path.distance[y-1]);
4392 point_to_point(elev_l, source.alt*METERS_PER_FOOT,
4393 destination.alt*METERS_PER_FOOT, LR.eps_dielect,
4394 LR.sgm_conductivity, LR.eno_ns_surfref, LR.frq_mhz,
4395 LR.radio_climate, LR.pol, LR.conf, LR.rel, loss,
4399 elevation=((acos(cos_test_angle))/deg2rad)-90.0;
4401 elevation=((acos(cos_xmtr_angle))/deg2rad)-90.0;
4403 /* Integrate the antenna's radiation
4404 pattern into the overall path loss. */
4406 x=(int)rint(10.0*(10.0-elevation));
4408 if (x>=0 && x<=1000)
4410 pattern=(double)LR.antenna_pattern[(int)azimuth][x];
4413 patterndB=20.0*log10(pattern);
4419 total_loss=loss-patterndB;
4423 fprintf(fd,"%f\t%f\n",KM_PER_MILE*(path.distance[path.length-1]-path.distance[y]),total_loss);
4424 fprintf(fd2,"%7.2f\t\t%7.2f\t\t %d\t%s\n",KM_PER_MILE*path.distance[y],total_loss, errnum, strmode);
4429 fprintf(fd,"%f\t%f\n",path.distance[path.length-1]-path.distance[y],total_loss);
4430 fprintf(fd2,"%7.2f\t\t%7.2f\t\t %d\t%s\n",path.distance[y],total_loss, errnum, strmode);
4435 if (total_loss>maxloss)
4438 if (total_loss<minloss)
4444 fprintf(fd2,"\nLongley-Rice path loss between %s and %s is %.2f dB.\n",source.name, destination.name, loss);
4447 fprintf(fd2,"Total path loss including TX antenna pattern is %.2f dB.\n",total_loss);
4451 fprintf(fd2,"\nNotes on \"errnum\":\n\n");
4452 fprintf(fd2," 0: No error. :-)\n");
4453 fprintf(fd2," 1: Warning! Some parameters are nearly out of range.\n");
4454 fprintf(fd2," Results should be used with caution.\n");
4455 fprintf(fd2," 2: Note: Default parameters have been substituted for impossible ones.\n");
4456 fprintf(fd2," 3: Warning! A combination of parameters is out of range.\n");
4457 fprintf(fd2," Results are probably invalid.\n");
4458 fprintf(fd2," Other: Warning! Some parameters are out of range.\n");
4459 fprintf(fd2," Results are probably invalid.\n\nEnd of Report\n");
4462 fprintf(stdout,"Longley-Rice path loss between %s and %s is %.2f dB.\n",source.name, destination.name, loss);
4465 fprintf(stdout,"Total path loss including TX antenna pattern is %.2f dB.\n",total_loss);
4467 fprintf(stdout,"Path Loss Report written to: \"%s\"\n",report_name);
4474 /* Default filename and output file type */
4476 strncpy(filename,"loss\0",5);
4477 strncpy(term,"png\0",4);
4478 strncpy(ext,"png\0",4);
4483 /* Grab extension and terminal type from "name" */
4485 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
4486 filename[x]=name[x];
4490 for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
4492 term[y]=tolower(name[x]);
4502 { /* No extension -- Default is png */
4505 strncpy(term,"png\0",4);
4506 strncpy(ext,"png\0",4);
4510 /* Either .ps or .postscript may be used
4511 as an extension for postscript output. */
4513 if (strncmp(term,"postscript",10)==0)
4514 strncpy(ext,"ps\0",3);
4516 else if (strncmp(ext,"ps",2)==0)
4517 strncpy(term,"postscript enhanced color\0",26);
4519 fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
4522 fd=fopen("splat.gp","w");
4524 fprintf(fd,"set grid\n");
4525 fprintf(fd,"set yrange [%2.3f to %2.3f]\n", minloss, maxloss);
4526 fprintf(fd,"set encoding iso_8859_1\n");
4527 fprintf(fd,"set term %s\n",term);
4528 fprintf(fd,"set title \"SPLAT! Loss Profile Along Path Between %s and %s (%.2f%c azimuth)\"\n",destination.name, source.name, Azimuth(destination,source),176);
4531 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(destination,source));
4533 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(destination,source));
4535 if (got_azimuth_pattern || got_elevation_pattern)
4536 fprintf(fd,"set ylabel \"Total Path Loss (including TX antenna pattern) (dB)");
4538 fprintf(fd,"set ylabel \"Longley-Rice Path Loss (dB)");
4540 fprintf(fd,"\"\nset output \"%s.%s\"\n",filename,ext);
4541 fprintf(fd,"plot \"profile.gp\" title \"Path Loss\" with lines\n");
4545 x=system("gnuplot splat.gp");
4550 unlink("profile.gp");
4551 unlink("reference.gp");
4553 fprintf(stdout," Done!\n");
4558 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
4561 void ObstructionReport(struct site xmtr, struct site rcvr, char report, double f)
4565 double h_r, h_t, h_x, h_r_orig, cos_tx_angle, cos_test_angle,
4566 cos_tx_angle_f1, cos_tx_angle_fpt6, haavt, d_tx, d_x,
4567 h_r_f1, h_r_fpt6, h_f, h_los, lambda=0.0, azimuth,
4568 pattern, patterndB, distance, angle1, angle2;
4569 char report_name[80], string[255], string_fpt6[255],
4573 sprintf(report_name,"%s-to-%s.txt",xmtr.name,rcvr.name);
4575 for (x=0; report_name[x]!=0; x++)
4576 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
4579 fd=fopen(report_name,"w");
4581 fprintf(fd,"\n\t\t--==[ SPLAT! v%s Obstruction Report ]==--\n\n",splat_version);
4582 fprintf(fd,"Analysis of great circle path between %s and %s:\n",xmtr.name, rcvr.name);
4583 fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
4584 fprintf(fd,"Transmitter site: %s\n",xmtr.name);
4589 fprintf(fd,"Site location: %.4f North / %.4f West",xmtr.lat, xmtr.lon);
4590 fprintf(fd, " (%s N / ", dec2dms(xmtr.lat));
4595 fprintf(fd,"Site location: %.4f South / %.4f West",-xmtr.lat, xmtr.lon);
4596 fprintf(fd, " (%s S / ", dec2dms(xmtr.lat));
4599 fprintf(fd, "%s W)\n", dec2dms(xmtr.lon));
4603 fprintf(fd,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(xmtr));
4604 fprintf(fd,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*xmtr.alt, METERS_PER_FOOT*(xmtr.alt+GetElevation(xmtr)));
4609 fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(xmtr));
4610 fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",xmtr.alt, xmtr.alt+GetElevation(xmtr));
4618 fprintf(fd,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4620 fprintf(fd,"Antenna height above average terrain: %.2f feet\n",haavt);
4625 distance=Distance(xmtr,rcvr);
4626 azimuth=Azimuth(xmtr,rcvr);
4627 angle1=ElevationAngle(xmtr,rcvr);
4628 angle2=ElevationAngle2(xmtr,rcvr,earthradius);
4630 if (got_azimuth_pattern || got_elevation_pattern)
4632 x=(int)rint(10.0*(10.0-angle2));
4634 if (x>=0 && x<=1000)
4635 pattern=(double)LR.antenna_pattern[(int)rint(azimuth)][x];
4639 fprintf(fd,"Antenna pattern toward %s: %.3f",rcvr.name,pattern);
4640 patterndB=20.0*log10(pattern);
4641 fprintf(fd," (%.2f dB)\n",patterndB);
4646 fprintf(fd,"Distance to %s: %.2f kilometers\n",rcvr.name,KM_PER_MILE*distance);
4649 fprintf(fd,"Distance to %s: %.2f miles\n",rcvr.name,distance);
4651 fprintf(fd,"Azimuth to %s: %.2f degrees\n",rcvr.name,azimuth);
4654 fprintf(fd,"Elevation angle to %s: %+.4f degrees\n",rcvr.name,angle1);
4657 fprintf(fd,"Depression angle to %s: %+.4f degrees\n",rcvr.name,angle1);
4662 fprintf(fd,"Depression");
4664 fprintf(fd,"Elevation");
4666 fprintf(fd," angle to the first obstruction: %+.4f degrees\n",angle2);
4669 fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
4673 fprintf(fd,"Receiver site: %s\n",rcvr.name);
4677 fprintf(fd,"Site location: %.4f North / %.4f West",rcvr.lat, rcvr.lon);
4678 fprintf(fd, " (%s N / ", dec2dms(rcvr.lat));
4683 fprintf(fd,"Site location: %.4f South / %.4f West",-rcvr.lat, rcvr.lon);
4684 fprintf(fd, " (%s S / ", dec2dms(rcvr.lat));
4687 fprintf(fd, "%s W)\n", dec2dms(rcvr.lon));
4691 fprintf(fd,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(rcvr));
4692 fprintf(fd,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*rcvr.alt, METERS_PER_FOOT*(rcvr.alt+GetElevation(rcvr)));
4697 fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(rcvr));
4698 fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",rcvr.alt, rcvr.alt+GetElevation(rcvr));
4706 fprintf(fd,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4708 fprintf(fd,"Antenna height above average terrain: %.2f feet\n",haavt);
4711 azimuth=Azimuth(rcvr,xmtr);
4712 angle1=ElevationAngle(rcvr,xmtr);
4713 angle2=ElevationAngle2(rcvr,xmtr,earthradius);
4716 fprintf(fd,"Distance to %s: %.2f kilometers\n",xmtr.name,KM_PER_MILE*distance);
4718 fprintf(fd,"Distance to %s: %.2f miles\n",xmtr.name,distance);
4720 fprintf(fd,"Azimuth to %s: %.2f degrees\n",xmtr.name,azimuth);
4723 fprintf(fd,"Elevation to %s: %+.4f degrees\n",xmtr.name,angle1);
4726 fprintf(fd,"Depression angle to %s: %+.4f degrees\n",xmtr.name,angle1);
4731 fprintf(fd,"Depression");
4733 fprintf(fd,"Elevation");
4735 fprintf(fd," angle to the first obstruction: %+.4f degrees\n",angle2);
4739 fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
4743 /* Generate profile of the terrain. Create the path
4744 from transmitter to receiver because that's the
4745 way the original los() function did it, and going
4746 the other way can yield slightly different results. */
4748 ReadPath(xmtr,rcvr);
4749 h_r=GetElevation(rcvr)+rcvr.alt+earthradius;
4753 h_t=GetElevation(xmtr)+xmtr.alt+earthradius;
4754 d_tx=5280.0*Distance(rcvr,xmtr);
4755 cos_tx_angle=((h_r*h_r)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r*d_tx);
4756 cos_tx_angle_f1=cos_tx_angle;
4757 cos_tx_angle_fpt6=cos_tx_angle;
4760 lambda=9.8425e8/(f*1e6);
4762 /* At each point along the path calculate the cosine
4763 of a sort of "inverse elevation angle" at the receiver.
4764 From the antenna, 0 deg. looks at the ground, and 90 deg.
4765 is parallel to the ground.
4767 Start at the receiver. If this is the lowest antenna,
4768 then terrain obstructions will be nearest to it. (Plus,
4769 that's the way the original los() did it.)
4771 Calculate cosines only. That's sufficient to compare
4772 angles and it saves the extra computational burden of
4773 acos(). However, note the inverted comparison: if
4774 acos(A) > acos(B), then B > A. */
4776 for (x=path.length-1; x>0; x--)
4778 site_x.lat=path.lat[x];
4779 site_x.lon=path.lon[x];
4782 h_x=GetElevation(site_x)+earthradius;
4783 d_x=5280.0*Distance(rcvr,site_x);
4785 /* Deal with the LOS path first. */
4787 cos_test_angle=((h_r*h_r)+(d_x*d_x)-(h_x*h_x))/(2.0*h_r*d_x);
4789 if (cos_tx_angle>cos_test_angle)
4792 fprintf(fd,"SPLAT! detected obstructions at:\n\n");
4794 if (site_x.lat>=0.0)
4797 fprintf(fd,"\t%.4f N, %.4f W, %5.2f kilometers, %6.2f meters AMSL\n",site_x.lat, site_x.lon, KM_PER_MILE*(d_x/5280.0), METERS_PER_FOOT*(h_x-earthradius));
4799 fprintf(fd,"\t%.4f N, %.4f W, %5.2f miles, %6.2f feet AMSL\n",site_x.lat, site_x.lon, d_x/5280.0, h_x-earthradius);
4805 fprintf(fd,"\t%.4f S, %.4f W, %5.2f kilometers, %6.2f meters AMSL\n",-site_x.lat, site_x.lon, KM_PER_MILE*(d_x/5280.0), METERS_PER_FOOT*(h_x-earthradius));
4808 fprintf(fd,"\t%.4f S, %.4f W, %5.2f miles, %6.2f feet AMSL\n",-site_x.lat, site_x.lon, d_x/5280.0, h_x-earthradius);
4812 while (cos_tx_angle>cos_test_angle)
4815 cos_test_angle=((h_r*h_r)+(d_x*d_x)-(h_x*h_x))/(2.0*h_r*d_x);
4816 cos_tx_angle=((h_r*h_r)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r*d_tx);
4821 /* Now clear the first Fresnel zone, but don't
4822 clutter the obstruction report. */
4824 cos_tx_angle_f1=((h_r_f1*h_r_f1)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r_f1*d_tx);
4825 h_los=sqrt(h_r_f1*h_r_f1+d_x*d_x-2*h_r_f1*d_x*cos_tx_angle_f1);
4826 h_f=h_los-sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4831 cos_tx_angle_f1=((h_r_f1*h_r_f1)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r_f1*d_tx);
4832 h_los=sqrt(h_r_f1*h_r_f1+d_x*d_x-2*h_r_f1*d_x*cos_tx_angle_f1);
4833 h_f=h_los-sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4836 /* And clear the 60% F1 zone. */
4838 cos_tx_angle_fpt6=((h_r_fpt6*h_r_fpt6)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r_fpt6*d_tx);
4839 h_los=sqrt(h_r_fpt6*h_r_fpt6+d_x*d_x-2*h_r_fpt6*d_x*cos_tx_angle_fpt6);
4840 h_f=h_los-0.6*sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4845 cos_tx_angle_fpt6=((h_r_fpt6*h_r_fpt6)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r_fpt6*d_tx);
4846 h_los=sqrt(h_r_fpt6*h_r_fpt6+d_x*d_x-2*h_r_fpt6*d_x*cos_tx_angle_fpt6);
4847 h_f=h_los-0.6*sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4855 sprintf(string,"\nAntenna at %s must be raised to at least %.2f meters AGL\nto clear all obstructions detected by SPLAT!\n",rcvr.name, METERS_PER_FOOT*(h_r-GetElevation(rcvr)-earthradius));
4857 sprintf(string,"\nAntenna at %s must be raised to at least %.2f feet AGL\nto clear all obstructions detected by SPLAT!\n",rcvr.name, h_r-GetElevation(rcvr)-earthradius);
4861 sprintf(string,"\nNo obstructions to LOS path due to terrain were detected by SPLAT!\n");
4865 if (h_r_fpt6>h_r_orig)
4868 sprintf(string_fpt6,"\nAntenna at %s must be raised to at least %.2f meters AGL\nto clear 60%c of the first Fresnel zone.\n",rcvr.name, METERS_PER_FOOT*(h_r_fpt6-GetElevation(rcvr)-earthradius),37);
4871 sprintf(string_fpt6,"\nAntenna at %s must be raised to at least %.2f feet AGL\nto clear 60%c of the first Fresnel zone.\n",rcvr.name, h_r_fpt6-GetElevation(rcvr)-earthradius,37);
4875 sprintf(string_fpt6,"\n60%c of the first Fresnel zone is clear.\n",37);
4877 if (h_r_f1>h_r_orig)
4880 sprintf(string_f1,"\nAntenna at %s must be raised to at least %.2f meters AGL\nto clear the first Fresnel zone.\n",rcvr.name, METERS_PER_FOOT*(h_r_f1-GetElevation(rcvr)-earthradius));
4883 sprintf(string_f1,"\nAntenna at %s must be raised to at least %.2f feet AGL\nto clear the first Fresnel zone.\n",rcvr.name, h_r_f1-GetElevation(rcvr)-earthradius);
4888 sprintf(string_f1,"\nThe first Fresnel zone is clear.\n\n");
4892 fprintf(fd,"%s",string);
4896 fprintf(fd,"%s",string_f1);
4897 fprintf(fd,"%s",string_fpt6);
4902 /* Display report summary on terminal */
4904 /* Line-of-sight status */
4906 fprintf(stdout,"%s",string);
4910 /* Fresnel zone status */
4912 fprintf(stdout,"%s",string_f1);
4913 fprintf(stdout,"%s",string_fpt6);
4916 fprintf(stdout, "\nObstruction report written to: \"%s\"\n",report_name);
4921 void SiteReport(struct site xmtr)
4923 char report_name[80];
4928 sprintf(report_name,"%s-site_report.txt",xmtr.name);
4930 for (x=0; report_name[x]!=0; x++)
4931 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
4934 fd=fopen(report_name,"w");
4936 fprintf(fd,"\n\t--==[ SPLAT! v%s Site Analysis Report For: %s ]==--\n\n",splat_version,xmtr.name);
4938 fprintf(fd,"---------------------------------------------------------------------------\n\n");
4942 fprintf(fd,"Site location: %.4f North / %.4f West",xmtr.lat, xmtr.lon);
4943 fprintf(fd, " (%s N / ",dec2dms(xmtr.lat));
4948 fprintf(fd,"Site location: %.4f South / %.4f West",-xmtr.lat, xmtr.lon);
4949 fprintf(fd, " (%s S / ",dec2dms(xmtr.lat));
4952 fprintf(fd, "%s W)\n",dec2dms(xmtr.lon));
4956 fprintf(fd,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(xmtr));
4957 fprintf(fd,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*xmtr.alt, METERS_PER_FOOT*(xmtr.alt+GetElevation(xmtr)));
4962 fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(xmtr));
4963 fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",xmtr.alt, xmtr.alt+GetElevation(xmtr));
4968 if (terrain>-4999.0)
4971 fprintf(fd,"Antenna height above average terrain: %.2f meters\n\n",METERS_PER_FOOT*terrain);
4973 fprintf(fd,"Antenna height above average terrain: %.2f feet\n\n",terrain);
4975 /* Display the average terrain between 2 and 10 miles
4976 from the transmitter site at azimuths of 0, 45, 90,
4977 135, 180, 225, 270, and 315 degrees. */
4979 for (azi=0; azi<=315; azi+=45)
4981 fprintf(fd,"Average terrain at %3d degrees azimuth: ",azi);
4982 terrain=AverageTerrain(xmtr,(double)azi,2.0,10.0);
4984 if (terrain>-4999.0)
4987 fprintf(fd,"%.2f meters AMSL\n",METERS_PER_FOOT*terrain);
4989 fprintf(fd,"%.2f feet AMSL\n",terrain);
4993 fprintf(fd,"No terrain\n");
4997 fprintf(fd,"\n---------------------------------------------------------------------------\n\n");
4999 fprintf(stdout,"\nSite analysis report written to: \"%s\"\n",report_name);
5002 void LoadTopoData(int max_lon, int min_lon, int max_lat, int min_lat)
5004 /* This function loads the SDF files required
5005 to cover the limits of the region specified. */
5007 int x, y, width, ymin, ymax;
5009 width=ReduceAngle(max_lon-min_lon);
5011 if ((max_lon-min_lon)<=180.0)
5013 for (y=0; y<=width; y++)
5014 for (x=min_lat; x<=max_lat; x++)
5016 ymin=(int)(min_lon+(double)y);
5032 sprintf(string,"%d:%d:%d:%d",x, x+1, ymin, ymax);
5039 for (y=0; y<=width; y++)
5040 for (x=min_lat; x<=max_lat; x++)
5058 sprintf(string,"%d:%d:%d:%d",x, x+1, ymin, ymax);
5064 int LoadPLI(char *filename)
5066 int error=0, max_west, min_west, max_north, min_north;
5067 char string[80], *pointer=NULL;
5068 double latitude=0.0, longitude=0.0, azimuth=0.0, elevation=0.0,
5072 fd=fopen(filename,"r");
5076 fgets(string,78,fd);
5077 pointer=strchr(string,';');
5082 sscanf(string,"%d, %d",&max_west, &min_west);
5084 fgets(string,78,fd);
5085 pointer=strchr(string,';');
5090 sscanf(string,"%d, %d",&max_north, &min_north);
5092 fgets(string,78,fd);
5093 pointer=strchr(string,';');
5098 LoadTopoData(max_west-1, min_west, max_north-1, min_north);
5100 fprintf(stdout,"\nReading \"%s\"... ",filename);
5103 fscanf(fd,"%lf, %lf, %lf, %lf, %lf",&latitude, &longitude, &azimuth, &elevation, &loss);
5117 if (loss<=(double)maxdB)
5118 OrMask(latitude,longitude,((unsigned char)(loss))<<3);
5120 fscanf(fd,"%lf, %lf, %lf, %lf, %lf",&latitude, &longitude, &azimuth, &elevation, &loss);
5125 fprintf(stdout," Done!\n");
5135 void WriteKML(struct site source, struct site destination)
5138 char block, report_name[80];
5139 double distance, rx_alt, tx_alt, cos_xmtr_angle,
5140 azimuth, cos_test_angle, test_alt;
5143 ReadPath(source,destination);
5145 sprintf(report_name,"%s-to-%s.kml",source.name,destination.name);
5147 for (x=0; report_name[x]!=0; x++)
5148 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
5151 fd=fopen(report_name,"w");
5153 fprintf(fd,"<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n");
5154 fprintf(fd,"<kml xmlns=\"http://earth.google.com/kml/2.0\">\n");
5155 fprintf(fd,"<!-- Generated by SPLAT! Version %s -->\n",splat_version);
5156 fprintf(fd,"<Folder>\n");
5157 fprintf(fd,"<name>SPLAT! Path</name>\n");
5158 fprintf(fd,"<open>1</open>\n");
5159 fprintf(fd,"<description>Path Between %s and %s</description>\n",source.name,destination.name);
5161 fprintf(fd,"<Placemark>\n");
5162 fprintf(fd," <name>%s</name>\n",source.name);
5163 fprintf(fd," <description>\n");
5164 fprintf(fd," Transmit Site\n");
5166 if (source.lat>=0.0)
5167 fprintf(fd," <BR>%s North</BR>\n",dec2dms(source.lat));
5169 fprintf(fd," <BR>%s South</BR>\n",dec2dms(source.lat));
5171 fprintf(fd," <BR>%s West</BR>\n",dec2dms(source.lon));
5173 azimuth=Azimuth(source,destination);
5174 distance=Distance(source,destination);
5177 fprintf(fd," <BR>%.2f km",distance*KM_PER_MILE);
5179 fprintf(fd," <BR>%.2f miles",distance);
5181 fprintf(fd," to %s</BR>\n <BR>toward an azimuth of %.2f%c</BR>\n",destination.name,azimuth,176);
5183 fprintf(fd," </description>\n");
5184 fprintf(fd," <visibility>1</visibility>\n");
5185 fprintf(fd," <Style>\n");
5186 fprintf(fd," <IconStyle>\n");
5187 fprintf(fd," <Icon>\n");
5188 fprintf(fd," <href>root://icons/palette-5.png</href>\n");
5189 fprintf(fd," <x>224</x>\n");
5190 fprintf(fd," <y>224</y>\n");
5191 fprintf(fd," <w>32</w>\n");
5192 fprintf(fd," <h>32</h>\n");
5193 fprintf(fd," </Icon>\n");
5194 fprintf(fd," </IconStyle>\n");
5195 fprintf(fd," </Style>\n");
5196 fprintf(fd," <Point>\n");
5197 fprintf(fd," <extrude>1</extrude>\n");
5198 fprintf(fd," <altitudeMode>relativeToGround</altitudeMode>\n");
5199 fprintf(fd," <coordinates>%f,%f,30</coordinates>\n",(source.lon<180.0?-source.lon:360.0-source.lon),source.lat);
5200 fprintf(fd," </Point>\n");
5201 fprintf(fd,"</Placemark>\n");
5203 fprintf(fd,"<Placemark>\n");
5204 fprintf(fd," <name>%s</name>\n",destination.name);
5205 fprintf(fd," <description>\n");
5206 fprintf(fd," Receive Site\n");
5208 if (destination.lat>=0.0)
5209 fprintf(fd," <BR>%s North</BR>\n",dec2dms(destination.lat));
5211 fprintf(fd," <BR>%s South</BR>\n",dec2dms(destination.lat));
5213 fprintf(fd," <BR>%s West</BR>\n",dec2dms(destination.lon));
5217 fprintf(fd," <BR>%.2f km",distance*KM_PER_MILE);
5219 fprintf(fd," <BR>%.2f miles",distance);
5221 fprintf(fd," to %s</BR>\n <BR>toward an azimuth of %.2f%c</BR>\n",source.name,Azimuth(destination,source),176);
5223 fprintf(fd," </description>\n");
5224 fprintf(fd," <visibility>1</visibility>\n");
5225 fprintf(fd," <Style>\n");
5226 fprintf(fd," <IconStyle>\n");
5227 fprintf(fd," <Icon>\n");
5228 fprintf(fd," <href>root://icons/palette-5.png</href>\n");
5229 fprintf(fd," <x>224</x>\n");
5230 fprintf(fd," <y>224</y>\n");
5231 fprintf(fd," <w>32</w>\n");
5232 fprintf(fd," <h>32</h>\n");
5233 fprintf(fd," </Icon>\n");
5234 fprintf(fd," </IconStyle>\n");
5235 fprintf(fd," </Style>\n");
5236 fprintf(fd," <Point>\n");
5237 fprintf(fd," <extrude>1</extrude>\n");
5238 fprintf(fd," <altitudeMode>relativeToGround</altitudeMode>\n");
5239 fprintf(fd," <coordinates>%f,%f,30</coordinates>\n",(destination.lon<180.0?-destination.lon:360.0-destination.lon),destination.lat);
5240 fprintf(fd," </Point>\n");
5241 fprintf(fd,"</Placemark>\n");
5243 fprintf(fd,"<Placemark>\n");
5244 fprintf(fd,"<name>Point-to-Point Path</name>\n");
5245 fprintf(fd," <visibility>1</visibility>\n");
5246 fprintf(fd," <open>0</open>\n");
5247 fprintf(fd," <Style>\n");
5248 fprintf(fd," <LineStyle>\n");
5249 fprintf(fd," <color>7fffffff</color>\n");
5250 fprintf(fd," </LineStyle>\n");
5251 fprintf(fd," <PolyStyle>\n");
5252 fprintf(fd," <color>7fffffff</color>\n");
5253 fprintf(fd," </PolyStyle>\n");
5254 fprintf(fd," </Style>\n");
5255 fprintf(fd," <LineString>\n");
5256 fprintf(fd," <extrude>1</extrude>\n");
5257 fprintf(fd," <tessellate>1</tessellate>\n");
5258 fprintf(fd," <altitudeMode>relativeToGround</altitudeMode>\n");
5259 fprintf(fd," <coordinates>\n");
5261 for (x=0; x<path.length; x++)
5262 fprintf(fd," %f,%f,5\n",(path.lon[x]<180.0?-path.lon[x]:360.0-path.lon[x]),path.lat[x]);
5264 fprintf(fd," </coordinates>\n");
5265 fprintf(fd," </LineString>\n");
5266 fprintf(fd,"</Placemark>\n");
5268 fprintf(fd,"<Placemark>\n");
5269 fprintf(fd,"<name>Line-of-Sight Path</name>\n");
5270 fprintf(fd," <visibility>1</visibility>\n");
5271 fprintf(fd," <open>0</open>\n");
5272 fprintf(fd," <Style>\n");
5273 fprintf(fd," <LineStyle>\n");
5274 fprintf(fd," <color>ff00ff00</color>\n");
5275 fprintf(fd," </LineStyle>\n");
5276 fprintf(fd," <PolyStyle>\n");
5277 fprintf(fd," <color>7f00ff00</color>\n");
5278 fprintf(fd," </PolyStyle>\n");
5279 fprintf(fd," </Style>\n");
5280 fprintf(fd," <LineString>\n");
5281 fprintf(fd," <extrude>1</extrude>\n");
5282 fprintf(fd," <tessellate>1</tessellate>\n");
5283 fprintf(fd," <altitudeMode>relativeToGround</altitudeMode>\n");
5284 fprintf(fd," <coordinates>\n");
5286 /* Walk across the "path", indentifying obstructions along the way */
5288 for (y=0; y<path.length; y++)
5290 distance=5280.0*path.distance[y];
5291 tx_alt=earthradius+source.alt+path.elevation[0];
5292 rx_alt=earthradius+destination.alt+path.elevation[y];
5294 /* Calculate the cosine of the elevation of the
5295 transmitter as seen at the temp rx point. */
5297 cos_xmtr_angle=((rx_alt*rx_alt)+(distance*distance)-(tx_alt*tx_alt))/(2.0*rx_alt*distance);
5299 for (x=y, block=0; x>=0 && block==0; x--)
5301 distance=5280.0*(path.distance[y]-path.distance[x]);
5302 test_alt=earthradius+path.elevation[x];
5304 cos_test_angle=((rx_alt*rx_alt)+(distance*distance)-(test_alt*test_alt))/(2.0*rx_alt*distance);
5306 /* Compare these two angles to determine if
5307 an obstruction exists. Since we're comparing
5308 the cosines of these angles rather than
5309 the angles themselves, the following "if"
5310 statement is reversed from what it would
5311 be if the actual angles were compared. */
5313 if (cos_xmtr_angle>cos_test_angle)
5318 fprintf(fd," %f,%f,-30\n",(path.lon[y]<180.0?-path.lon[y]:360.0-path.lon[y]),path.lat[y]);
5320 fprintf(fd," %f,%f,5\n",(path.lon[y]<180.0?-path.lon[y]:360.0-path.lon[y]),path.lat[y]);
5323 fprintf(fd," </coordinates>\n");
5324 fprintf(fd," </LineString>\n");
5325 fprintf(fd,"</Placemark>\n");
5327 fprintf(fd," <LookAt>\n");
5328 fprintf(fd," <longitude>%f</longitude>\n",(source.lon<180.0?-source.lon:360.0-source.lon));
5329 fprintf(fd," <latitude>%f</latitude>\n",source.lat);
5330 fprintf(fd," <range>300.0</range>\n");
5331 fprintf(fd," <tilt>45.0</tilt>\n");
5332 fprintf(fd," <heading>%f</heading>\n",azimuth);
5333 fprintf(fd," </LookAt>\n");
5335 fprintf(fd,"</Folder>\n");
5336 fprintf(fd,"</kml>\n");
5340 fprintf(stdout, "KML file written to: \"%s\"\n",report_name);
5345 int main(char argc, char *argv[])
5347 int x, y, z=0, min_lat, min_lon, max_lat, max_lon,
5348 rxlat, rxlon, txlat, txlon, west_min, west_max,
5349 north_min, north_max;
5351 unsigned char coverage=0, LRmap=0, ext[20], terrain_plot=0,
5352 elevation_plot=0, height_plot=0,
5353 longley_plot=0, cities=0, bfs=0, txsites=0,
5354 count, report='y', norm=0, topomap=0, geo=0,
5357 char mapfile[255], header[80], city_file[5][255],
5358 elevation_file[255], height_file[255],
5359 longley_file[255], terrain_file[255],
5360 string[255], rxfile[255], *env=NULL,
5361 txfile[255], map=0, boundary_file[5][255],
5362 udt_file[255], rxsite=0, plo_filename[255],
5363 pli_filename[255], nf=0;
5365 double altitude=0.0, altitudeLR=0.0, tx_range=0.0,
5366 rx_range=0.0, deg_range=0.0, deg_limit,
5367 deg_range_lon, er_mult, freq=0.0;
5369 struct site tx_site[4], rx_site;
5376 fprintf(stdout,"\n\t\t --==[ SPLAT! v%s Available Options... ]==--\n\n",splat_version);
5377 fprintf(stdout," -t txsite(s).qth (max of 4)\n");
5378 fprintf(stdout," -r rxsite.qth\n");
5379 fprintf(stdout," -c plot coverage of TX(s) with an RX antenna at X feet/meters AGL\n");
5380 fprintf(stdout," -L plot path loss map of TX based on an RX at X feet/meters AGL\n");
5381 fprintf(stdout," -s filename(s) of city/site file(s) to import (5 max)\n");
5382 fprintf(stdout," -b filename(s) of cartographic boundary file(s) to import (max of 5)\n");
5383 fprintf(stdout," -p filename of terrain profile graph to plot\n");
5384 fprintf(stdout," -e filename of terrain elevation graph to plot\n");
5385 fprintf(stdout," -h filename of terrain height graph to plot\n");
5386 fprintf(stdout," -H filename of normalized terrain height graph to plot\n");
5387 fprintf(stdout," -l filename of Longley-Rice graph to plot\n");
5388 fprintf(stdout," -o filename of topographic map to generate (.ppm)\n");
5389 fprintf(stdout," -u filename of user-defined terrain file to import\n");
5390 fprintf(stdout," -d sdf file directory path (overrides path in ~/.splat_path file)\n");
5391 fprintf(stdout," -n no analysis, brief report\n");
5392 fprintf(stdout," -N no analysis, no report\n");
5393 fprintf(stdout," -m earth radius multiplier\n");
5394 fprintf(stdout," -f frequency for Fresnel zone calculation (MHz)\n");
5395 fprintf(stdout," -R modify default range for -c or -L (miles/kilometers)\n");
5396 fprintf(stdout," -db maximum loss contour to display on path loss maps (80-230 dB)\n");
5397 fprintf(stdout," -nf do not plot Fresnel zones in height plots\n");
5398 fprintf(stdout," -plo filename of path-loss output file\n");
5399 fprintf(stdout," -pli filename of path-loss input file\n");
5400 fprintf(stdout," -udt filename of user defined terrain input file\n");
5401 fprintf(stdout," -geo generate an Xastir .geo georeference file (with .ppm output)\n");
5402 fprintf(stdout," -kml generate a Google Earth .kml file (for point-to-point links)\n");
5403 fprintf(stdout," -metric employ metric rather than imperial units for all user I/O\n\n");
5405 fprintf(stdout,"If that flew by too fast, consider piping the output through 'less':\n");
5406 fprintf(stdout,"\n\tsplat | less\n\n");
5407 fprintf(stdout,"Type 'man splat', or see the documentation for more details.\n\n");
5419 elevation_file[0]=0;
5429 earthradius=EARTHRADIUS;
5431 sprintf(header,"\n\t\t--==[ Welcome To SPLAT! v%s ]==--\n\n", splat_version);
5435 tx_site[x].lat=91.0;
5436 tx_site[x].lon=361.0;
5439 for (x=0; x<MAXSLOTS; x++)
5441 dem[x].min_el=32768;
5442 dem[x].max_el=-32768;
5443 dem[x].min_north=90;
5444 dem[x].max_north=-90;
5445 dem[x].min_west=360;
5449 /* Scan for command line arguments */
5451 for (x=1; x<=y; x++)
5453 if (strcmp(argv[x],"-R")==0)
5457 if (z<=y && argv[z][0] && argv[z][0]!='-')
5459 sscanf(argv[z],"%lf",&max_range);
5464 if (max_range>1000.0)
5469 if (strcmp(argv[x],"-m")==0)
5473 if (z<=y && argv[z][0] && argv[z][0]!='-')
5475 sscanf(argv[z],"%lf",&er_mult);
5483 earthradius*=er_mult;
5487 if (strcmp(argv[x],"-o")==0)
5491 if (z<=y && argv[z][0] && argv[z][0]!='-')
5492 strncpy(mapfile,argv[z],253);
5496 if (strcmp(argv[x],"-u")==0)
5500 if (z<=y && argv[z][0] && argv[z][0]!='-')
5501 strncpy(udt_file,argv[z],253);
5504 if (strcmp(argv[x],"-c")==0)
5508 if (z<=y && argv[z][0] && argv[z][0]!='-')
5510 sscanf(argv[z],"%lf",&altitude);
5515 if (strcmp(argv[x],"-db")==0 || strcmp(argv[x],"-dB")==0)
5519 if (z<=y && argv[z][0] && argv[z][0]!='-')
5521 sscanf(argv[z],"%d",&maxdB);
5533 if (strcmp(argv[x],"-p")==0)
5537 if (z<=y && argv[z][0] && argv[z][0]!='-')
5539 strncpy(terrain_file,argv[z],253);
5544 if (strcmp(argv[x],"-e")==0)
5548 if (z<=y && argv[z][0] && argv[z][0]!='-')
5550 strncpy(elevation_file,argv[z],253);
5555 if (strcmp(argv[x],"-h")==0 || strcmp(argv[x],"-H")==0)
5559 if (z<=y && argv[z][0] && argv[z][0]!='-')
5561 strncpy(height_file,argv[z],253);
5565 if (strcmp(argv[x],"-H")==0)
5571 if (strcmp(argv[x],"-n")==0)
5577 if (strcmp(argv[x],"-N")==0)
5583 if (strcmp(argv[x],"-metric")==0)
5586 if (strcmp(argv[x],"-geo")==0)
5589 if (strcmp(argv[x],"-kml")==0)
5592 if (strcmp(argv[x],"-nf")==0)
5595 if (strcmp(argv[x],"-d")==0)
5599 if (z<=y && argv[z][0] && argv[z][0]!='-')
5600 strncpy(sdf_path,argv[z],253);
5603 if (strcmp(argv[x],"-t")==0)
5605 /* Read Transmitter Location */
5609 while (z<=y && argv[z][0] && argv[z][0]!='-' && txsites<4)
5611 strncpy(txfile,argv[z],253);
5612 tx_site[txsites]=LoadQTH(txfile);
5620 if (strcmp(argv[x],"-L")==0)
5624 if (z<=y && argv[z][0] && argv[z][0]!='-')
5626 sscanf(argv[z],"%lf",&altitudeLR);
5629 fprintf(stdout,"c and L are exclusive options, ignoring L.\n");
5638 if (strcmp(argv[x],"-l")==0)
5642 if (z<=y && argv[z][0] && argv[z][0]!='-')
5644 strncpy(longley_file,argv[z],253);
5646 /* Doing this twice is harmless */
5651 if (strcmp(argv[x],"-r")==0)
5653 /* Read Receiver Location */
5657 if (z<=y && argv[z][0] && argv[z][0]!='-')
5659 strncpy(rxfile,argv[z],253);
5660 rx_site=LoadQTH(rxfile);
5665 if (strcmp(argv[x],"-s")==0)
5667 /* Read city file(s) */
5671 while (z<=y && argv[z][0] && argv[z][0]!='-' && cities<5)
5673 strncpy(city_file[cities],argv[z],253);
5681 if (strcmp(argv[x],"-b")==0)
5683 /* Read Boundary File(s) */
5687 while (z<=y && argv[z][0] && argv[z][0]!='-' && bfs<5)
5689 strncpy(boundary_file[bfs],argv[z],253);
5697 if (strcmp(argv[x],"-f")==0)
5701 if (z<=y && argv[z][0] && argv[z][0]!='-')
5703 sscanf(argv[z],"%lf",&freq);
5713 if (strcmp(argv[x],"-plo")==0)
5717 if (z<=y && argv[z][0] && argv[z][0]!='-')
5718 strncpy(plo_filename,argv[z],253);
5721 if (strcmp(argv[x],"-pli")==0)
5725 if (z<=y && argv[z][0] && argv[z][0]!='-')
5726 strncpy(pli_filename,argv[z],253);
5730 /* Perform some error checking on the arguments
5731 and switches parsed from the command-line.
5732 If an error is encountered, print a message
5733 and exit gracefully. */
5737 fprintf(stderr,"\n%c*** ERROR: No transmitter site(s) specified!\n\n",7);
5741 for (x=0, y=0; x<txsites; x++)
5743 if (tx_site[x].lat==91.0 && tx_site[x].lon==361.0)
5745 fprintf(stderr,"\n*** ERROR: Transmitter site #%d not found!",x+1);
5752 fprintf(stderr,"%c\n\n",7);
5756 if ((coverage+LRmap+pli_filename[0])==0 && rx_site.lat==91.0 && rx_site.lon==361.0)
5758 if (max_range!=0.0 && txsites!=0)
5760 /* Plot topographic map of radius "max_range" */
5769 fprintf(stderr,"\n%c*** ERROR: No receiver site found or specified!\n\n",7);
5774 /* No major errors were detected. Whew! :-) */
5776 /* Adjust input parameters if -metric option is used */
5780 altitudeLR/=METERS_PER_FOOT; /* meters --> feet */
5781 max_range/=KM_PER_MILE; /* kilometers --> miles */
5782 altitude/=METERS_PER_FOOT; /* meters --> feet */
5785 /* If no SDF path was specified on the command line (-d), check
5786 for a path specified in the $HOME/.splat_path file. If the
5787 file is not found, then sdf_path[] remains NULL, and the
5788 current working directory is assumed to contain the SDF
5794 sprintf(string,"%s/.splat_path",env);
5795 fd=fopen(string,"r");
5799 fgets(string,253,fd);
5801 /* Remove <CR> and/or <LF> from string */
5803 for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0 && x<253; x++);
5806 strncpy(sdf_path,string,253);
5812 /* Ensure a trailing '/' is present in sdf_path */
5818 if (sdf_path[x-1]!='/' && x!=0)
5825 fprintf(stdout,"%s",header);
5828 if (pli_filename[0])
5830 y=LoadPLI(pli_filename);
5832 for (x=0; x<txsites; x++)
5833 PlaceMarker(tx_site[x]);
5836 PlaceMarker(rx_site);
5840 for (x=0; x<bfs; x++)
5841 LoadBoundaries(boundary_file[x]);
5846 for (x=0; x<cities; x++)
5847 LoadCities(city_file[x]);
5850 WritePPMLR(mapfile,geo);
5861 min_lon=(int)floor(tx_site[0].lon);
5862 max_lon=(int)floor(tx_site[0].lon);
5864 for (y=0, z=0; z<txsites; z++)
5866 txlat=(int)floor(tx_site[z].lat);
5867 txlon=(int)floor(tx_site[z].lon);
5875 if (LonDiff(txlon,min_lon)<0.0)
5878 if (LonDiff(txlon,max_lon)>0.0)
5884 rxlat=(int)floor(rx_site.lat);
5885 rxlon=(int)floor(rx_site.lon);
5893 if (LonDiff(rxlon,min_lon)<0.0)
5896 if (LonDiff(rxlon,max_lon)>0.0)
5900 /* Load the required SDF files */
5902 LoadTopoData(max_lon, min_lon, max_lat, min_lat);
5904 if (coverage | LRmap | topomap)
5909 for (z=0; z<txsites; z++)
5911 /* "Ball park" estimates used to load any additional
5912 SDF files required to conduct this analysis. */
5914 tx_range=sqrt(1.5*(tx_site[z].alt+GetElevation(tx_site[z])));
5917 rx_range=sqrt(1.5*altitudeLR);
5919 rx_range=sqrt(1.5*altitude);
5921 /* deg_range determines the maximum
5922 amount of topo data we read */
5924 deg_range=(tx_range+rx_range)/69.0;
5926 /* max_range sets the maximum size of the
5927 analysis. A small, non-zero amount can
5928 be used to shrink the size of the analysis
5929 and limit the amount of topo data read by
5930 SPLAT! A very large number will only increase
5931 the width of the analysis, not the size of
5935 max_range=tx_range+rx_range;
5937 if (max_range<(tx_range+rx_range))
5938 deg_range=max_range/69.0;
5940 /* Prevent the demand for a really wide coverage
5941 from allocating more slots than are available
5946 case 2: deg_limit=0.25;
5949 case 4: deg_limit=0.5;
5952 case 9: deg_limit=1.0;
5955 case 16: deg_limit=2.0;
5958 case 25: deg_limit=3.0;
5961 if (tx_site[z].lat<70.0)
5962 deg_range_lon=deg_range/cos(deg2rad*tx_site[z].lat);
5964 deg_range_lon=deg_range/cos(deg2rad*70.0);
5966 /* Correct for squares in degrees not being square in miles */
5968 if (deg_range>deg_limit)
5969 deg_range=deg_limit;
5971 if (deg_range_lon>deg_limit)
5972 deg_range_lon=deg_limit;
5975 north_min=(int)floor(tx_site[z].lat-deg_range);
5976 north_max=(int)floor(tx_site[z].lat+deg_range);
5978 west_min=(int)floor(tx_site[z].lon-deg_range_lon);
5983 while (west_min>=360)
5986 west_max=(int)floor(tx_site[z].lon+deg_range_lon);
5991 while (west_max>=360)
5994 if (north_min<min_lat)
5997 if (north_max>max_lat)
6000 if (LonDiff(west_min,min_lon)<0.0)
6003 if (LonDiff(west_max,max_lon)>0.0)
6007 /* Load any additional SDF files, if required */
6009 LoadTopoData(max_lon, min_lon, max_lat, min_lat);
6015 if (mapfile[0] && topomap==0)
6018 if (freq==0.0 && nf==0)
6023 if (coverage | LRmap)
6025 for (x=0; x<txsites; x++)
6028 PlotCoverage(tx_site[x],altitude);
6031 PlotLRMap(tx_site[x],altitudeLR,plo_filename);
6033 PlaceMarker(tx_site[x]);
6036 SiteReport(tx_site[x]);
6042 if (coverage==0 && LRmap==0)
6044 PlaceMarker(rx_site);
6046 for (x=0; x<txsites; x++)
6048 PlaceMarker(tx_site[x]);
6055 PlotPath(tx_site[x],rx_site,1);
6059 PlotPath(tx_site[x],rx_site,8);
6063 PlotPath(tx_site[x],rx_site,16);
6067 PlotPath(tx_site[x],rx_site,32);
6072 ObstructionReport(tx_site[x],rx_site,report,freq);
6075 WriteKML(tx_site[x],rx_site);
6083 for (x=0; x<bfs; x++)
6084 LoadBoundaries(boundary_file[x]);
6089 for (x=0; x<cities; x++)
6090 LoadCities(city_file[x]);
6094 WritePPMLR(mapfile,geo);
6096 WritePPM(mapfile,geo);
6103 for (x=0; terrain_file[x]!='.' && terrain_file[x]!=0 && x<80; x++);
6105 if (terrain_file[x]=='.') /* extension */
6108 for (y=1, z=x, x++; terrain_file[x]!=0 && x<253 && y<14; x++, y++)
6109 ext[y]=terrain_file[x];
6117 ext[0]=0; /* No extension */
6121 for (count=0; count<txsites; count++)
6123 sprintf(string,"%s-%c%s%c",terrain_file,'1'+count,ext,0);
6124 GraphTerrain(tx_site[count],rx_site,string);
6129 GraphTerrain(tx_site[0],rx_site,terrain_file);
6136 for (x=0; elevation_file[x]!='.' && elevation_file[x]!=0 && x<80; x++);
6138 if (elevation_file[x]=='.') /* extension */
6141 for (y=1, z=x, x++; elevation_file[x]!=0 && x<253 && y<14; x++, y++)
6142 ext[y]=elevation_file[x];
6145 elevation_file[z]=0;
6150 ext[0]=0; /* No extension */
6151 elevation_file[x]=0;
6154 for (count=0; count<txsites; count++)
6156 sprintf(string,"%s-%c%s%c",elevation_file,'1'+count,ext,0);
6157 GraphElevation(tx_site[count],rx_site,string);
6162 GraphElevation(tx_site[0],rx_site,elevation_file);
6169 for (x=0; height_file[x]!='.' && height_file[x]!=0 && x<80; x++);
6171 if (height_file[x]=='.') /* extension */
6174 for (y=1, z=x, x++; height_file[x]!=0 && x<253 && y<14; x++, y++)
6175 ext[y]=height_file[x];
6183 ext[0]=0; /* No extension */
6187 for (count=0; count<txsites; count++)
6189 sprintf(string,"%s-%c%s%c",height_file,'1'+count,ext,0);
6190 GraphHeight(tx_site[count],rx_site,string,freq,norm);
6195 GraphHeight(tx_site[0],rx_site,height_file,freq,norm);
6202 for (x=0; longley_file[x]!='.' && longley_file[x]!=0 && x<80; x++);
6204 if (longley_file[x]=='.') /* extension */
6207 for (y=1, z=x, x++; longley_file[x]!=0 && x<253 && y<14; x++, y++)
6208 ext[y]=longley_file[x];
6216 ext[0]=0; /* No extension */
6220 for (count=0; count<txsites; count++)
6222 sprintf(string,"%s-%c%s%c",longley_file,'1'+count,ext,0);
6223 GraphLongley(tx_site[count],rx_site,string);
6228 GraphLongley(tx_site[0],rx_site,longley_file);