1 /****************************************************************************
2 * SPLAT: An RF Signal Propagation Loss and Terrain Analysis Tool *
3 * Last update: 15-Mar-2007 *
4 *****************************************************************************
5 * Project started in 1997 by John A. Magliacane, KD2BD *
6 *****************************************************************************
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.0b"};
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, four_thirds_earth;
2465 ReadPath(source,destination);
2467 four_thirds_earth=EARTHRADIUS*(4.0/3.0);
2469 /* Copy elevations along path into the elev_l[] array. */
2471 for (x=0; x<path.length; x++)
2472 elev_l[x+2]=path.elevation[x]*METERS_PER_FOOT;
2474 /* Since the only energy the Longley-Rice model considers
2475 reaching the destination is based on what is scattered
2476 or deflected from the first obstruction along the path,
2477 we first need to find the location and elevation angle
2478 of that first obstruction (if it exists). This is done
2479 using a 4/3rds Earth radius to match the model used by
2480 Longley-Rice. This information is required for properly
2481 integrating the antenna's elevation pattern into the
2482 calculation for overall path loss. (Using path.length-1
2483 below avoids a Longley-Rice model error from occuring at
2484 the destination point.) */
2486 for (y=2; (y<(path.length-1) && path.distance[y]<=max_range); y++)
2488 /* Process this point only if it
2489 has not already been processed. */
2491 if (GetMask(path.lat[y],path.lon[y])==0)
2493 distance=5280.0*path.distance[y];
2494 source_alt=four_thirds_earth+source.alt+path.elevation[0];
2495 dest_alt=four_thirds_earth+destination.alt+path.elevation[y];
2496 dest_alt2=dest_alt*dest_alt;
2497 source_alt2=source_alt*source_alt;
2499 /* Calculate the cosine of the elevation of
2500 the receiver as seen by the transmitter. */
2502 cos_xmtr_angle=((source_alt2)+(distance*distance)-(dest_alt2))/(2.0*source_alt*distance);
2504 if (got_elevation_pattern || fd!=NULL)
2506 /* If no antenna elevation pattern is available, and
2507 no output file is designated, the following code
2508 that determines the elevation angle to the first
2509 obstruction along the path is bypassed. */
2511 for (x=2, block=0; (x<y && block==0); x++)
2513 distance=5280.0*path.distance[x];
2514 test_alt=four_thirds_earth+path.elevation[x];
2516 /* Calculate the cosine of the elevation
2517 angle of the terrain (test point)
2518 as seen by the transmitter. */
2520 cos_test_angle=((source_alt2)+(distance*distance)-(test_alt*test_alt))/(2.0*source_alt*distance);
2522 /* Compare these two angles to determine if
2523 an obstruction exists. Since we're comparing
2524 the cosines of these angles rather than
2525 the angles themselves, the sense of the
2526 following "if" statement is reversed from
2527 what it would be if the angles themselves
2530 if (cos_xmtr_angle>cos_test_angle)
2534 /* At this point, we have the elevation angle
2535 to the first obstruction (if it exists). */
2538 /* Determine attenuation for each point along the
2539 path using Longley-Rice's point_to_point mode
2540 starting at y=2 (number_of_points = 1), the
2541 shortest distance terrain can play a role in
2544 elev_l[0]=y-1; /* (number of points - 1) */
2546 /* Distance between elevation samples */
2547 elev_l[1]=METERS_PER_MILE*(path.distance[y]-path.distance[y-1]);
2549 point_to_point(elev_l,source.alt*METERS_PER_FOOT,
2550 destination.alt*METERS_PER_FOOT, LR.eps_dielect,
2551 LR.sgm_conductivity, LR.eno_ns_surfref, LR.frq_mhz,
2552 LR.radio_climate, LR.pol, LR.conf, LR.rel, loss,
2556 elevation=((acos(cos_test_angle))/deg2rad)-90.0;
2559 elevation=((acos(cos_xmtr_angle))/deg2rad)-90.0;
2561 temp.lat=path.lat[y];
2562 temp.lon=path.lon[y];
2564 azimuth=(Azimuth(source,temp));
2568 /* Write path loss data to output file */
2570 fprintf(fd,"%.7f, %.7f, %.3f, %.3f, %.2f\n",path.lat[y], path.lon[y], azimuth, elevation, loss);
2573 /* Integrate the antenna's radiation
2574 pattern into the overall path loss. */
2576 x=(int)rint(10.0*(10.0-elevation));
2578 if (x>=0 && x<=1000)
2580 azimuth=rint(azimuth);
2582 pattern=(double)LR.antenna_pattern[(int)azimuth][x];
2586 pattern=20.0*log10(pattern);
2601 OrMask(path.lat[y],path.lon[y],((unsigned char)(loss))<<3);
2604 else if (GetMask(path.lat[y],path.lon[y])==0 && path.distance[y]>max_range)
2605 OrMask(path.lat[y],path.lon[y],1);
2609 void PlotCoverage(struct site source, double altitude)
2611 /* This function performs a 360 degree sweep around the
2612 transmitter site (source location), and plots the
2613 line-of-sight coverage of the transmitter on the SPLAT!
2614 generated topographic map based on a receiver located
2615 at the specified altitude (in feet AGL). Results
2616 are stored in memory, and written out in the form
2617 of a topographic map when the WritePPM() function
2618 is later invoked. */
2620 float lat, lon, one_pixel;
2621 static unsigned char mask_value;
2624 unsigned char symbol[4], x;
2626 /* Initialize mask_value */
2628 if (mask_value!=8 && mask_value!=16 && mask_value!=32)
2631 one_pixel=1.0/1200.0;
2640 fprintf(stdout,"\nComputing line-of-sight coverage of %s with an RX antenna\nat %.2f %s AGL:\n\n 0%c to 25%c ",source.name,metric?altitude*METERS_PER_FOOT:altitude,metric?"meters":"feet",37,37);
2643 /* 18.75=1200 pixels/degree divided by 64 loops
2644 per progress indicator symbol (.oOo) printed. */
2646 z=(int)(18.75*ReduceAngle(max_west-min_west));
2648 for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2657 PlotPath(source,edge,mask_value);
2662 fprintf(stdout,"%c",symbol[x]);
2674 fprintf(stdout,"\n25%c to 50%c ",37,37);
2677 z=(int)(18.75*(max_north-min_north));
2679 for (lat=max_north, x=0; lat>=min_north; lat-=one_pixel)
2685 PlotPath(source,edge,mask_value);
2690 fprintf(stdout,"%c",symbol[x]);
2702 fprintf(stdout,"\n50%c to 75%c ",37,37);
2705 z=(int)(18.75*ReduceAngle(max_west-min_west));
2707 for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2716 PlotPath(source,edge,mask_value);
2721 fprintf(stdout,"%c",symbol[x]);
2733 fprintf(stdout,"\n75%c to 100%c ",37,37);
2736 z=(int)(18.75*(max_north-min_north));
2738 for (lat=min_north, x=0; lat<=max_north; lat+=one_pixel)
2744 PlotPath(source,edge,mask_value);
2749 fprintf(stdout,"%c",symbol[x]);
2760 fprintf(stdout,"\nDone!\n");
2763 /* Assign next mask value */
2780 void PlotLRMap(struct site source, double altitude, char *plo_filename)
2782 /* This function performs a 360 degree sweep around the
2783 transmitter site (source location), and plots the
2784 Longley-Rice attenuation on the SPLAT! generated
2785 topographic map based on a receiver located at
2786 the specified altitude (in feet AGL). Results
2787 are stored in memory, and written out in the form
2788 of a topographic map when the WritePPMLR() function
2789 is later invoked. */
2793 float lat, lon, one_pixel;
2794 unsigned char symbol[4], x;
2797 one_pixel=1.0/1200.0;
2806 fprintf(stdout,"\nComputing Longley-Rice coverage of %s ", source.name);
2808 fprintf(stdout,"out to a radius\nof %.2f %s with an RX antenna at %.2f %s AGL:\n\n 0%c to 25%c ",metric?max_range*KM_PER_MILE:max_range,metric?"kilometers":"miles",metric?altitude*METERS_PER_FOOT:altitude,metric?"meters":"feet",37,37);
2811 if (plo_filename[0]!=0)
2812 fd=fopen(plo_filename,"wb");
2816 /* Write header information to output file */
2818 fprintf(fd,"%d, %d\t; max_west, min_west\n%d, %d\t; max_north, min_north\n",max_west, min_west, max_north, min_north);
2821 /* 18.75=1200 pixels/degree divided by 64 loops
2822 per progress indicator symbol (.oOo) printed. */
2824 z=(int)(18.75*ReduceAngle(max_west-min_west));
2826 for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2835 PlotLRPath(source,edge,fd);
2840 fprintf(stdout,"%c",symbol[x]);
2852 fprintf(stdout,"\n25%c to 50%c ",37,37);
2855 z=(int)(18.75*(max_north-min_north));
2857 for (lat=max_north, x=0; lat>=min_north; lat-=one_pixel)
2863 PlotLRPath(source,edge,fd);
2868 fprintf(stdout,"%c",symbol[x]);
2880 fprintf(stdout,"\n50%c to 75%c ",37,37);
2883 z=(int)(18.75*ReduceAngle(max_west-min_west));
2885 for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
2894 PlotLRPath(source,edge,fd);
2899 fprintf(stdout,"%c",symbol[x]);
2911 fprintf(stdout,"\n75%c to 100%c ",37,37);
2914 z=(int)(18.75*(max_north-min_north));
2916 for (lat=min_north, x=0; lat<=max_north; lat+=one_pixel)
2922 PlotLRPath(source,edge,fd);
2927 fprintf(stdout,"%c",symbol[x]);
2941 fprintf(stdout,"\nDone!\n");
2945 void WritePPM(char *filename, unsigned char geo)
2947 /* This function generates a topographic map in Portable Pix Map
2948 (PPM) format based on logarithmically scaled topology data,
2949 as well as the content of flags held in the mask[][] array.
2950 The image created is rotated counter-clockwise 90 degrees
2951 from its representation in dem[][] so that north points
2952 up and east points right in the image generated. */
2954 char mapfile[255], geofile[255];
2955 unsigned char found, mask;
2956 unsigned width, height, output;
2957 int indx, x, y, x0=0, y0=0, minlat, minlon;
2958 float lat, lon, one_pixel, conversion, one_over_gamma;
2961 one_pixel=1.0/1200.0;
2962 one_over_gamma=1.0/GAMMA;
2963 conversion=255.0/pow((double)(max_elevation-min_elevation),one_over_gamma);
2965 width=(unsigned)(1200*ReduceAngle(max_west-min_west));
2966 height=(unsigned)(1200*ReduceAngle(max_north-min_north));
2969 strncpy(mapfile, "map.ppm\0",8);
2972 for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
2974 mapfile[x]=filename[x];
2975 geofile[x]=filename[x];
2992 fd=fopen(geofile,"wb");
2994 fprintf(fd,"FILENAME\t%s\n",mapfile);
2995 fprintf(fd,"#\t\tX\tY\tLong\t\tLat\n");
2996 fprintf(fd,"TIEPOINT\t0\t0\t%d.000\t\t%d.000\n",(max_west<180?-max_west:360-max_west),max_north);
2997 fprintf(fd,"TIEPOINT\t%u\t%u\t%d.000\t\t%d.000\n",width-1,height-1,(min_west<180?-min_west:360-min_west),min_north);
2998 fprintf(fd,"IMAGESIZE\t%u\t%u\n",width,height);
2999 fprintf(fd,"#\n# Auto Generated by SPLAT! v%s\n#\n",splat_version);
3004 fd=fopen(mapfile,"wb");
3006 fprintf(fd,"P6\n%u %u\n255\n",width,height);
3008 fprintf(stdout,"\nWriting \"%s\" (%ux%u pixmap image)... ",mapfile,width,height);
3011 for (y=0, lat=((double)max_north)-one_pixel; y<(int)height; y++, lat-=one_pixel)
3013 minlat=(int)floor(lat);
3015 for (x=0, lon=((double)max_west)-one_pixel; x<(int)width; x++, lon-=one_pixel)
3020 minlon=(int)floor(lon);
3022 for (indx=0, found=0; indx<MAXSLOTS && found==0;)
3023 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
3030 x0=(int)(1199.0*(lat-floor(lat)));
3031 y0=(int)(1199.0*(lon-floor(lon)));
3033 mask=dem[indx].mask[x0][y0];
3036 /* Text Labels: Red */
3037 fprintf(fd,"%c%c%c",255,0,0);
3040 /* County Boundaries: Light Cyan */
3041 fprintf(fd,"%c%c%c",128,128,255);
3043 else switch (mask&57)
3047 fprintf(fd,"%c%c%c",0,255,0);
3052 fprintf(fd,"%c%c%c",0,255,255);
3056 /* TX1 + TX2: Yellow */
3057 fprintf(fd,"%c%c%c",255,255,0);
3061 /* TX3: Medium Violet */
3062 fprintf(fd,"%c%c%c",147,112,219);
3066 /* TX1 + TX3: Pink */
3067 fprintf(fd,"%c%c%c",255,192,203);
3071 /* TX2 + TX3: Orange */
3072 fprintf(fd,"%c%c%c",255,165,0);
3076 /* TX1 + TX2 + TX3: Dark Green */
3077 fprintf(fd,"%c%c%c",0,100,0);
3082 fprintf(fd,"%c%c%c",255,130,71);
3086 /* TX1 + TX4: Green Yellow */
3087 fprintf(fd,"%c%c%c",173,255,47);
3091 /* TX2 + TX4: Dark Sea Green 1 */
3092 fprintf(fd,"%c%c%c",193,255,193);
3096 /* TX1 + TX2 + TX4: Blanched Almond */
3097 fprintf(fd,"%c%c%c",255,235,205);
3101 /* TX3 + TX4: Dark Turquoise */
3102 fprintf(fd,"%c%c%c",0,206,209);
3106 /* TX1 + TX3 + TX4: Medium Spring Green */
3107 fprintf(fd,"%c%c%c",0,250,154);
3111 /* TX2 + TX3 + TX4: Tan */
3112 fprintf(fd,"%c%c%c",210,180,140);
3116 /* TX1 + TX2 + TX3 + TX4: Gold2 */
3117 fprintf(fd,"%c%c%c",238,201,0);
3121 /* Water: Medium Blue */
3122 if (dem[indx].data[x0][y0]==0)
3123 fprintf(fd,"%c%c%c",0,0,170);
3126 /* Elevation: Greyscale */
3127 output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
3128 fprintf(fd,"%c%c%c",output,output,output);
3135 /* We should never get here, but if */
3136 /* we do, display the region as black */
3138 fprintf(fd,"%c%c%c",0,0,0);
3144 fprintf(stdout,"Done!\n");
3148 void WritePPMLR(char *filename, unsigned char geo)
3150 /* This function generates a topographic map in Portable Pix Map
3151 (PPM) format based on the content of flags held in the mask[][]
3152 array (only). The image created is rotated counter-clockwise
3153 90 degrees from its representation in dem[][] so that north
3154 points up and east points right in the image generated. */
3156 char mapfile[255], geofile[255];
3157 unsigned width, height, output;
3158 unsigned char found, mask, cityorcounty;
3159 int indx, x, y, t, t2, x0, y0, minlat, minlon, loss;
3160 float lat, lon, one_pixel, conversion, one_over_gamma;
3163 one_pixel=1.0/1200.0;
3164 one_over_gamma=1.0/GAMMA;
3165 conversion=255.0/pow((double)(max_elevation-min_elevation),one_over_gamma);
3167 width=(unsigned)(1200*ReduceAngle(max_west-min_west));
3168 height=(unsigned)(1200*ReduceAngle(max_north-min_north));
3171 strncpy(mapfile, "map.ppm\0",8);
3174 for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
3176 mapfile[x]=filename[x];
3177 geofile[x]=filename[x];
3194 fd=fopen(geofile,"wb");
3196 fprintf(fd,"FILENAME\t%s\n",mapfile);
3197 fprintf(fd,"#\t\tX\tY\tLong\t\tLat\n");
3198 fprintf(fd,"TIEPOINT\t0\t0\t%d.000\t\t%d.000\n",(max_west<180?-max_west:360-max_west),max_north);
3199 fprintf(fd,"TIEPOINT\t%u\t%u\t%d.000\t\t%.3f\n",width-1,height+29,(min_west<180?-min_west:360-min_west),(double)(min_north-0.025));
3200 fprintf(fd,"IMAGESIZE\t%u\t%u\n",width,height+30);
3201 fprintf(fd,"#\n# Auto Generated by SPLAT! v%s\n#\n",splat_version);
3206 fd=fopen(mapfile,"wb");
3208 fprintf(fd,"P6\n%u %u\n255\n",width,height+30);
3210 fprintf(stdout,"\nWriting \"%s\" (%ux%u pixmap image)... ",mapfile,width,height+30);
3213 for (y=0, lat=((double)max_north)-one_pixel; y<(int)height; y++, lat-=one_pixel)
3215 minlat=(int)floor(lat);
3217 for (x=0, lon=((double)max_west)-one_pixel; x<(int)width; x++, lon-=one_pixel)
3222 minlon=(int)floor(lon);
3224 for (indx=0, found=0; indx<MAXSLOTS && found==0;)
3225 if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
3231 x0=(int)(1199.0*(lat-floor(lat)));
3232 y0=(int)(1199.0*(lon-floor(lon)));
3234 mask=dem[indx].mask[x0][y0];
3235 loss=70+(10*(int)((mask&248)>>3));
3240 /* Text Labels - Black or Red */
3242 if ((mask&120) && (loss<=90))
3243 fprintf(fd,"%c%c%c",0,0,0);
3245 fprintf(fd,"%c%c%c",255,0,0);
3252 /* County Boundaries: Black */
3254 fprintf(fd,"%c%c%c",0,0,0);
3259 if (cityorcounty==0)
3263 { /* Display land or sea elevation */
3265 if (dem[indx].data[x0][y0]==0)
3266 fprintf(fd,"%c%c%c",0,0,170);
3269 /* Elevation: Greyscale */
3270 output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
3271 fprintf(fd,"%c%c%c",output,output,output);
3277 /* Plot signal loss in color */
3280 fprintf(fd,"%c%c%c",255,0,0);
3284 fprintf(fd,"%c%c%c",255,128,0);
3288 fprintf(fd,"%c%c%c",255,165,0);
3292 fprintf(fd,"%c%c%c",255,206,0);
3296 fprintf(fd,"%c%c%c",255,255,0);
3300 fprintf(fd,"%c%c%c",184,255,0);
3304 fprintf(fd,"%c%c%c",0,255,0);
3308 fprintf(fd,"%c%c%c",0,208,0);
3312 fprintf(fd,"%c%c%c",0,196,196);
3316 fprintf(fd,"%c%c%c",0,148,255);
3320 fprintf(fd,"%c%c%c",80,80,255);
3324 fprintf(fd,"%c%c%c",0,38,255);
3328 fprintf(fd,"%c%c%c",142,63,255);
3332 fprintf(fd,"%c%c%c",196,54,255);
3336 fprintf(fd,"%c%c%c",255,0,255);
3340 fprintf(fd,"%c%c%c",255,194,204);
3345 if (dem[indx].data[x0][y0]==0)
3346 fprintf(fd,"%c%c%c",0,0,170);
3349 /* Elevation: Greyscale */
3350 output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
3351 fprintf(fd,"%c%c%c",output,output,output);
3359 /* We should never get here, but if */
3360 /* we do, display the region as black */
3362 fprintf(fd,"%c%c%c",0,0,0);
3367 /* Display legend along bottom of image */
3371 for (y0=0; y0<30; y0++)
3373 for (indx=0; indx<16; indx++)
3375 for (x=0; x<x0; x++)
3380 if (y0>=10 && y0<=18)
3385 if (smallfont[t2/10][y0-10][x-11])
3390 if (smallfont[t2%10][y0-10][x-19])
3394 if (smallfont[0][y0-10][x-27])
3401 fprintf(fd,"%c%c%c",255,0,0);
3405 fprintf(fd,"%c%c%c",255,128,0);
3409 fprintf(fd,"%c%c%c",255,165,0);
3413 fprintf(fd,"%c%c%c",255,206,0);
3417 fprintf(fd,"%c%c%c",255,255,0);
3421 fprintf(fd,"%c%c%c",184,255,0);
3425 fprintf(fd,"%c%c%c",0,255,0);
3429 fprintf(fd,"%c%c%c",0,208,0);
3433 fprintf(fd,"%c%c%c",0,196,196);
3437 fprintf(fd,"%c%c%c",0,148,255);
3441 fprintf(fd,"%c%c%c",80,80,255);
3445 fprintf(fd,"%c%c%c",0,38,255);
3449 fprintf(fd,"%c%c%c",142,63,255);
3453 fprintf(fd,"%c%c%c",196,54,255);
3457 fprintf(fd,"%c%c%c",255,0,255);
3462 fprintf(fd,"%c%c%c",0,0,0);
3466 fprintf(fd,"%c%c%c",255,194,204);
3473 fprintf(stdout,"Done!\n");
3477 void GraphTerrain(struct site source, struct site destination, char *name)
3479 /* This function invokes gnuplot to generate an appropriate
3480 output file indicating the terrain profile between the source
3481 and destination locations. "filename" is the name assigned
3482 to the output file generated by gnuplot. The filename extension
3483 is used to set gnuplot's terminal setting and output file type.
3484 If no extension is found, .png is assumed. */
3487 char filename[255], term[30], ext[15];
3490 ReadPath(destination,source);
3492 fd=fopen("profile.gp","wb");
3494 for (x=0; x<path.length; x++)
3497 fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*path.elevation[x]);
3500 fprintf(fd,"%f\t%f\n",path.distance[x],path.elevation[x]);
3507 /* Default filename and output file type */
3509 strncpy(filename,"profile\0",8);
3510 strncpy(term,"png\0",4);
3511 strncpy(ext,"png\0",4);
3516 /* Grab extension and terminal type from "name" */
3518 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
3519 filename[x]=name[x];
3523 for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
3525 term[y]=tolower(name[x]);
3535 { /* No extension -- Default is png */
3538 strncpy(term,"png\0",4);
3539 strncpy(ext,"png\0",4);
3543 /* Either .ps or .postscript may be used
3544 as an extension for postscript output. */
3546 if (strncmp(term,"postscript",10)==0)
3547 strncpy(ext,"ps\0",3);
3549 else if (strncmp(ext,"ps",2)==0)
3550 strncpy(term,"postscript enhanced color\0",26);
3552 fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
3555 fd=fopen("splat.gp","w");
3556 fprintf(fd,"set grid\n");
3557 fprintf(fd,"set autoscale\n");
3558 fprintf(fd,"set encoding iso_8859_1\n");
3559 fprintf(fd,"set term %s\n",term);
3560 fprintf(fd,"set title \"SPLAT! Terrain Profile Between %s and %s (%.2f%c Azimuth)\"\n",destination.name, source.name, Azimuth(destination,source),176);
3564 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(source,destination));
3565 fprintf(fd,"set ylabel \"Ground Elevation Above Sea Level (meters)\"\n");
3572 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(source,destination));
3573 fprintf(fd,"set ylabel \"Ground Elevation Above Sea Level (feet)\"\n");
3576 fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
3577 fprintf(fd,"plot \"profile.gp\" title \"\" with lines\n");
3580 x=system("gnuplot splat.gp");
3585 unlink("profile.gp");
3586 fprintf(stdout," Done!\n");
3591 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
3594 void GraphElevation(struct site source, struct site destination, char *name)
3596 /* This function invokes gnuplot to generate an appropriate
3597 output file indicating the terrain profile between the source
3598 and destination locations. "filename" is the name assigned
3599 to the output file generated by gnuplot. The filename extension
3600 is used to set gnuplot's terminal setting and output file type.
3601 If no extension is found, .png is assumed. */
3604 char filename[255], term[30], ext[15];
3605 double angle, refangle, maxangle=-90.0;
3607 FILE *fd=NULL, *fd2=NULL;
3609 ReadPath(destination,source); /* destination=RX, source=TX */
3610 refangle=ElevationAngle(destination,source);
3612 fd=fopen("profile.gp","wb");
3613 fd2=fopen("reference.gp","wb");
3615 for (x=1; x<path.length-1; x++)
3617 remote.lat=path.lat[x];
3618 remote.lon=path.lon[x];
3620 angle=ElevationAngle(destination,remote);
3624 fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[x],angle);
3625 fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[x],refangle);
3630 fprintf(fd,"%f\t%f\n",path.distance[x],angle);
3631 fprintf(fd2,"%f\t%f\n",path.distance[x],refangle);
3640 fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],refangle);
3641 fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],refangle);
3646 fprintf(fd,"%f\t%f\n",path.distance[path.length-1],refangle);
3647 fprintf(fd2,"%f\t%f\n",path.distance[path.length-1],refangle);
3655 /* Default filename and output file type */
3657 strncpy(filename,"profile\0",8);
3658 strncpy(term,"png\0",4);
3659 strncpy(ext,"png\0",4);
3664 /* Grab extension and terminal type from "name" */
3666 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
3667 filename[x]=name[x];
3671 for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
3673 term[y]=tolower(name[x]);
3683 { /* No extension -- Default is png */
3686 strncpy(term,"png\0",4);
3687 strncpy(ext,"png\0",4);
3691 /* Either .ps or .postscript may be used
3692 as an extension for postscript output. */
3694 if (strncmp(term,"postscript",10)==0)
3695 strncpy(ext,"ps\0",3);
3697 else if (strncmp(ext,"ps",2)==0)
3698 strncpy(term,"postscript enhanced color\0",26);
3700 fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
3703 fd=fopen("splat.gp","w");
3705 fprintf(fd,"set grid\n");
3706 fprintf(fd,"set yrange [%2.3f to %2.3f]\n", (-fabs(refangle)-0.25), maxangle+0.25);
3707 fprintf(fd,"set encoding iso_8859_1\n");
3708 fprintf(fd,"set term %s\n",term);
3709 fprintf(fd,"set title \"SPLAT! Elevation Profile Between %s and %s (%.2f%c azimuth)\"\n",destination.name,source.name,Azimuth(destination,source),176);
3712 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(source,destination));
3714 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(source,destination));
3717 fprintf(fd,"set ylabel \"Elevation Angle Along LOS Path Between %s and %s (degrees)\"\n",destination.name,source.name);
3718 fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
3719 fprintf(fd,"plot \"profile.gp\" title \"Real Earth Profile\" with lines, \"reference.gp\" title \"Line of Sight Path (%.2f%c elevation)\" with lines\n",refangle,176);
3723 x=system("gnuplot splat.gp");
3728 unlink("profile.gp");
3729 unlink("reference.gp");
3731 fprintf(stdout," Done!\n");
3736 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
3739 void GraphHeight(struct site source, struct site destination, char *name, double f, unsigned char n)
3741 /* This function invokes gnuplot to generate an appropriate
3742 output file indicating the terrain profile between the source
3743 and destination locations referenced to the line-of-sight path
3744 between the receive and transmit sites. "filename" is the name
3745 assigned to the output file generated by gnuplot. The filename
3746 extension is used to set gnuplot's terminal setting and output
3747 file type. If no extension is found, .png is assumed. */
3750 char filename[255], term[30], ext[15];
3751 double a, b, c, height=0.0, refangle, cangle, maxheight=-100000.0,
3752 minheight=100000.0, lambda=0.0, f_zone=0.0, fpt6_zone=0.0,
3753 nm=0.0, nb=0.0, ed=0.0, es=0.0, r=0.0, d=0.0, d1=0.0,
3754 terrain, azimuth, distance, dheight=0.0, minterrain=100000.0,
3755 minearth=100000.0, miny, maxy, min2y, max2y;
3757 FILE *fd=NULL, *fd2=NULL, *fd3=NULL, *fd4=NULL, *fd5=NULL;
3759 ReadPath(destination,source); /* destination=RX, source=TX */
3760 azimuth=Azimuth(destination,source);
3761 distance=Distance(destination,source);
3762 refangle=ElevationAngle(destination,source);
3763 b=GetElevation(destination)+destination.alt+earthradius;
3765 /* Wavelength and path distance (great circle) in feet. */
3769 lambda=9.8425e8/(f*1e6);
3770 d=5280.0*path.distance[path.length-1];
3775 ed=GetElevation(destination);
3776 es=GetElevation(source);
3777 nb=-destination.alt-ed;
3778 nm=(-source.alt-es-nb)/(path.distance[path.length-1]);
3781 fd=fopen("profile.gp","wb");
3782 fd2=fopen("reference.gp","wb");
3783 fd5=fopen("curvature.gp", "wb");
3787 fd3=fopen("fresnel.gp", "wb");
3788 fd4=fopen("fresnel_pt_6.gp", "wb");
3791 for (x=0; x<path.length; x++)
3793 remote.lat=path.lat[x];
3794 remote.lon=path.lon[x];
3797 terrain=GetElevation(remote);
3800 terrain+=destination.alt; /* RX antenna spike */
3802 a=terrain+earthradius;
3803 cangle=5280.0*Distance(destination,remote)/earthradius;
3804 c=b*sin(refangle*deg2rad+HALFPI)/sin(HALFPI-refangle*deg2rad-cangle);
3808 /* Per Fink and Christiansen, Electronics
3809 * Engineers' Handbook, 1989:
3811 * H = sqrt(lamba * d1 * (d - d1)/d)
3813 * where H is the distance from the LOS
3814 * path to the first Fresnel zone boundary.
3819 d1=5280.0*path.distance[x];
3820 f_zone=-1*sqrt(lambda*d1*(d-d1)/d);
3821 fpt6_zone=0.6*f_zone;
3826 r=-(nm*path.distance[x])-nb;
3841 fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*height);
3842 fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*r);
3843 fprintf(fd5,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*(height-terrain));
3848 fprintf(fd,"%f\t%f\n",path.distance[x],height);
3849 fprintf(fd2,"%f\t%f\n",path.distance[x],r);
3850 fprintf(fd5,"%f\t%f\n",path.distance[x],height-terrain);
3857 fprintf(fd3,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*f_zone);
3858 fprintf(fd4,"%f\t%f\n",KM_PER_MILE*path.distance[x],METERS_PER_FOOT*fpt6_zone);
3863 fprintf(fd3,"%f\t%f\n",path.distance[x],f_zone);
3864 fprintf(fd4,"%f\t%f\n",path.distance[x],fpt6_zone);
3867 if (f_zone<minheight)
3871 if (height>maxheight)
3874 if (height<minheight)
3880 if (terrain<minterrain)
3883 if ((height-terrain)<minearth)
3884 minearth=height-terrain;
3888 r=-(nm*path.distance[path.length-1])-nb;
3894 fprintf(fd,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3895 fprintf(fd2,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3900 fprintf(fd,"%f\t%f\n",path.distance[path.length-1],r);
3901 fprintf(fd2,"%f\t%f\n",path.distance[path.length-1],r);
3908 fprintf(fd3,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3909 fprintf(fd4,"%f\t%f\n",KM_PER_MILE*path.distance[path.length-1],METERS_PER_FOOT*r);
3914 fprintf(fd3,"%f\t%f\n",path.distance[path.length-1],r);
3915 fprintf(fd4,"%f\t%f\n",path.distance[path.length-1],r);
3937 /* Default filename and output file type */
3939 strncpy(filename,"height\0",8);
3940 strncpy(term,"png\0",4);
3941 strncpy(ext,"png\0",4);
3946 /* Grab extension and terminal type from "name" */
3948 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
3949 filename[x]=name[x];
3953 for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
3955 term[y]=tolower(name[x]);
3965 { /* No extension -- Default is png */
3968 strncpy(term,"png\0",4);
3969 strncpy(ext,"png\0",4);
3973 /* Either .ps or .postscript may be used
3974 as an extension for postscript output. */
3976 if (strncmp(term,"postscript",10)==0)
3977 strncpy(ext,"ps\0",3);
3979 else if (strncmp(ext,"ps",2)==0)
3980 strncpy(term,"postscript enhanced color\0",26);
3982 fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
3985 fd=fopen("splat.gp","w");
3987 dheight=maxheight-minheight;
3988 miny=minheight-0.15*dheight;
3989 maxy=maxheight+0.05*dheight;
3994 dheight=maxheight-minheight;
3995 min2y=miny-minterrain+0.05*dheight;
3999 miny-=min2y-minearth+0.05*dheight;
4000 min2y=minearth-0.05*dheight;
4003 max2y=min2y+maxy-miny;
4005 fprintf(fd,"set grid\n");
4006 fprintf(fd,"set yrange [%2.3f to %2.3f]\n", metric?miny*METERS_PER_FOOT:miny, metric?maxy*METERS_PER_FOOT:maxy);
4007 fprintf(fd,"set y2range [%2.3f to %2.3f]\n", metric?min2y*METERS_PER_FOOT:min2y, metric?max2y*METERS_PER_FOOT:max2y);
4008 fprintf(fd,"set xrange [-0.5 to %2.3f]\n",metric?KM_PER_MILE*rint(distance+0.5):rint(distance+0.5));
4009 fprintf(fd,"set encoding iso_8859_1\n");
4010 fprintf(fd,"set term %s\n",term);
4013 fprintf(fd,"set title \"SPLAT! Path Profile Between %s and %s (%.2f%c azimuth)\\nWith First Fresnel Zone\"\n",destination.name, source.name, azimuth,176);
4016 fprintf(fd,"set title \"SPLAT! Height Profile Between %s and %s (%.2f%c azimuth)\"\n",destination.name, source.name, azimuth,176);
4019 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(source,destination));
4021 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(source,destination));
4026 fprintf(fd,"set ylabel \"Normalized Height Referenced To LOS Path Between\\n%s and %s (meters)\"\n",destination.name,source.name);
4029 fprintf(fd,"set ylabel \"Normalized Height Referenced To LOS Path Between\\n%s and %s (feet)\"\n",destination.name,source.name);
4036 fprintf(fd,"set ylabel \"Height Referenced To LOS Path Between %s and %s (meters)\"\n",destination.name,source.name);
4039 fprintf(fd,"set ylabel \"Height Referenced To LOS Path Between %s and %s (feet)\"\n",destination.name,source.name);
4042 fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
4045 fprintf(fd,"plot \"profile.gp\" title \"Point-to-Point Profile\" with lines, \"reference.gp\" title \"Line of Sight Path\" with lines, \"curvature.gp\" axes x1y2 title \"Earth's Curvature Contour\" with lines, \"fresnel.gp\" axes x1y1 title \"First Fresnel Zone (%.3f MHz)\" with lines, \"fresnel_pt_6.gp\" title \"60%% of First Fresnel Zone\" with lines\n",f);
4048 fprintf(fd,"plot \"profile.gp\" title \"Point-to-Point Profile\" with lines, \"reference.gp\" title \"Line Of Sight Path\" with lines, \"curvature.gp\" axes x1y2 title \"Earth's Curvature Contour\" with lines\n");
4052 x=system("gnuplot splat.gp");
4057 unlink("profile.gp");
4058 unlink("reference.gp");
4059 unlink("curvature.gp");
4063 unlink("fresnel.gp");
4064 unlink("fresnel_pt_6.gp");
4067 fprintf(stdout," Done!\n");
4072 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
4075 void GraphLongley(struct site source, struct site destination, char *name)
4077 /* This function invokes gnuplot to generate an appropriate
4078 output file indicating the Longley-Rice model loss between
4079 the source and destination locations. "filename" is the
4080 name assigned to the output file generated by gnuplot.
4081 The filename extension is used to set gnuplot's terminal
4082 setting and output file type. If no extension is found,
4085 int x, y, z, errnum, errflag=0;
4086 char filename[255], term[30], ext[15], strmode[100],
4087 report_name[80], block=0;
4088 double maxloss=-100000.0, minloss=100000.0, loss, haavt,
4089 angle1, angle2, azimuth, pattern=1.0, patterndB=0.0,
4090 total_loss=0.0, cos_xmtr_angle, cos_test_angle=0.0,
4091 source_alt, test_alt, dest_alt, source_alt2, dest_alt2,
4092 distance, elevation, four_thirds_earth;
4093 FILE *fd=NULL, *fd2=NULL;
4095 sprintf(report_name,"%s-to-%s.lro",source.name,destination.name);
4097 four_thirds_earth=EARTHRADIUS*(4.0/3.0);
4099 for (x=0; report_name[x]!=0; x++)
4100 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
4103 fd2=fopen(report_name,"w");
4105 fprintf(fd2,"\n\t--==[ SPLAT! v%s Longley-Rice Model Path Loss Report ]==--\n\n",splat_version);
4106 fprintf(fd2,"Analysis of RF path conditions between %s and %s:\n",source.name, destination.name);
4107 fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4108 fprintf(fd2,"Transmitter site: %s\n",source.name);
4110 if (source.lat>=0.0)
4112 fprintf(fd2,"Site location: %.4f North / %.4f West",source.lat, source.lon);
4113 fprintf(fd2, " (%s N / ", dec2dms(source.lat));
4119 fprintf(fd2,"Site location: %.4f South / %.4f West",-source.lat, source.lon);
4120 fprintf(fd2, " (%s S / ", dec2dms(source.lat));
4123 fprintf(fd2, "%s W)\n", dec2dms(source.lon));
4127 fprintf(fd2,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(source));
4128 fprintf(fd2,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*source.alt,METERS_PER_FOOT*(source.alt+GetElevation(source)));
4133 fprintf(fd2,"Ground elevation: %.2f feet AMSL\n",GetElevation(source));
4134 fprintf(fd2,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",source.alt, source.alt+GetElevation(source));
4142 fprintf(fd2,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4144 fprintf(fd2,"Antenna height above average terrain: %.2f feet\n",haavt);
4147 azimuth=Azimuth(source,destination);
4148 angle1=ElevationAngle(source,destination);
4149 angle2=ElevationAngle2(source,destination,earthradius);
4151 if (got_azimuth_pattern || got_elevation_pattern)
4153 x=(int)rint(10.0*(10.0-angle2));
4155 if (x>=0 && x<=1000)
4156 pattern=(double)LR.antenna_pattern[(int)rint(azimuth)][x];
4158 patterndB=20.0*log10(pattern);
4160 fprintf(fd2,"Antenna pattern between %s and %s: %.3f (%.2f dB)\n",source.name, destination.name, pattern, patterndB);
4165 fprintf(fd2,"Distance to %s: %.2f kilometers\n",destination.name,METERS_PER_FOOT*Distance(source,destination));
4168 fprintf(fd2,"Distance to %s: %.2f miles.\n",destination.name,Distance(source,destination));
4170 fprintf(fd2,"Azimuth to %s: %.2f degrees\n",destination.name,azimuth);
4173 fprintf(fd2,"Elevation angle to %s: %+.4f degrees\n",destination.name,angle1);
4176 fprintf(fd2,"Depression angle to %s: %+.4f degrees\n",destination.name,angle1);
4181 fprintf(fd2,"Depression");
4183 fprintf(fd2,"Elevation");
4185 fprintf(fd2," angle to the first obstruction: %+.4f degrees\n",angle2);
4188 fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4192 fprintf(fd2,"Receiver site: %s\n",destination.name);
4194 if (destination.lat>=0.0)
4196 fprintf(fd2,"Site location: %.4f North / %.4f West",destination.lat, destination.lon);
4197 fprintf(fd2, " (%s N / ", dec2dms(destination.lat));
4202 fprintf(fd2,"Site location: %.4f South / %.4f West",-destination.lat, destination.lon);
4203 fprintf(fd2, " (%s S / ", dec2dms(destination.lat));
4206 fprintf(fd2, "%s W)\n", dec2dms(destination.lon));
4210 fprintf(fd2,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(destination));
4211 fprintf(fd2,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*destination.alt, METERS_PER_FOOT*(destination.alt+GetElevation(destination)));
4216 fprintf(fd2,"Ground elevation: %.2f feet AMSL\n",GetElevation(destination));
4217 fprintf(fd2,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",destination.alt, destination.alt+GetElevation(destination));
4220 haavt=haat(destination);
4225 fprintf(fd2,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4227 fprintf(fd2,"Antenna height above average terrain: %.2f feet\n",haavt);
4231 fprintf(fd2,"Distance to %s: %.2f kilometers\n",source.name,KM_PER_MILE*Distance(source,destination));
4234 fprintf(fd2,"Distance to %s: %.2f miles\n",source.name,Distance(source,destination));
4236 azimuth=Azimuth(destination,source);
4238 angle1=ElevationAngle(destination,source);
4239 angle2=ElevationAngle2(destination,source,earthradius);
4241 fprintf(fd2,"Azimuth to %s: %.2f degrees.\n",source.name,azimuth);
4244 fprintf(fd2,"Elevation angle to %s: %+.4f degrees\n",source.name,angle1);
4247 fprintf(fd2,"Depression angle to %s: %+.4f degrees\n",source.name,angle1);
4252 fprintf(fd2,"Depression");
4254 fprintf(fd2,"Elevation");
4256 fprintf(fd2," angle to the first obstruction: %+.4f degrees\n",angle2);
4259 fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4261 fprintf(fd2,"Longley-Rice path calculation parameters used in this analysis:\n\n");
4262 fprintf(fd2,"Earth's Dielectric Constant: %.3lf\n",LR.eps_dielect);
4263 fprintf(fd2,"Earth's Conductivity: %.3lf\n",LR.sgm_conductivity);
4264 fprintf(fd2,"Atmospheric Bending Constant (N): %.3lf\n",LR.eno_ns_surfref);
4265 fprintf(fd2,"Frequency: %.3lf (MHz)\n",LR.frq_mhz);
4266 fprintf(fd2,"Radio Climate: %d (",LR.radio_climate);
4268 switch (LR.radio_climate)
4271 fprintf(fd2,"Equatorial");
4275 fprintf(fd2,"Continental Subtropical");
4279 fprintf(fd2,"Maritime Subtropical");
4283 fprintf(fd2,"Desert");
4287 fprintf(fd2,"Continental Temperate");
4291 fprintf(fd2,"Martitime Temperate, Over Land");
4295 fprintf(fd2,"Maritime Temperate, Over Sea");
4299 fprintf(fd2,"Unknown");
4302 fprintf(fd2,")\nPolarization: %d (",LR.pol);
4305 fprintf(fd2,"Horizontal");
4308 fprintf(fd2,"Vertical");
4310 fprintf(fd2,")\nFraction of Situations: %.1lf%c\n",LR.conf*100.0,37);
4311 fprintf(fd2,"Fraction of Time: %.1lf%c\n",LR.rel*100.0,37);
4312 fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
4314 fprintf(fd2,"Analysis Results:\n\n");
4316 ReadPath(source, destination); /* source=TX, destination=RX */
4318 /* Copy elevations along path into the elev_l[] array. */
4320 for (x=0; x<path.length; x++)
4321 elev_l[x+2]=path.elevation[x]*METERS_PER_FOOT;
4324 fprintf(fd2,"Distance (km)");
4326 fprintf(fd2,"Distance (mi)");
4328 fprintf(fd2,"\tLoss (dB)\tErrnum\tComment\n\n");
4330 fd=fopen("profile.gp","w");
4332 azimuth=rint(Azimuth(source,destination));
4334 for (y=2; y<(path.length-1); y++) /* path.length-1 avoids LR error */
4336 distance=5280.0*path.distance[y];
4337 source_alt=four_thirds_earth+source.alt+path.elevation[0];
4338 dest_alt=four_thirds_earth+destination.alt+path.elevation[y];
4339 dest_alt2=dest_alt*dest_alt;
4340 source_alt2=source_alt*source_alt;
4342 /* Calculate the cosine of the elevation of
4343 the receiver as seen by the transmitter. */
4345 cos_xmtr_angle=((source_alt2)+(distance*distance)-(dest_alt2))/(2.0*source_alt*distance);
4347 if (got_elevation_pattern)
4349 /* If an antenna elevation pattern is available, the
4350 following code determines the elevation angle to
4351 the first obstruction along the path. */
4353 for (x=2, block=0; x<y && block==0; x++)
4355 distance=5280.0*(path.distance[y]-path.distance[x]);
4356 test_alt=four_thirds_earth+path.elevation[x];
4358 /* Calculate the cosine of the elevation
4359 angle of the terrain (test point)
4360 as seen by the transmitter. */
4362 cos_test_angle=((source_alt2)+(distance*distance)-(test_alt*test_alt))/(2.0*source_alt*distance);
4364 /* Compare these two angles to determine if
4365 an obstruction exists. Since we're comparing
4366 the cosines of these angles rather than
4367 the angles themselves, the sense of the
4368 following "if" statement is reversed from
4369 what it would be if the angles themselves
4372 if (cos_xmtr_angle>cos_test_angle)
4376 /* At this point, we have the elevation angle
4377 to the first obstruction (if it exists). */
4380 /* Determine path loss for each point along the
4381 path using Longley-Rice's point_to_point mode
4382 starting at x=2 (number_of_points = 1), the
4383 shortest distance terrain can play a role in
4386 elev_l[0]=y-1; /* (number of points - 1) */
4388 /* Distance between elevation samples */
4389 elev_l[1]=METERS_PER_MILE*(path.distance[y]-path.distance[y-1]);
4391 point_to_point(elev_l, source.alt*METERS_PER_FOOT,
4392 destination.alt*METERS_PER_FOOT, LR.eps_dielect,
4393 LR.sgm_conductivity, LR.eno_ns_surfref, LR.frq_mhz,
4394 LR.radio_climate, LR.pol, LR.conf, LR.rel, loss,
4398 elevation=((acos(cos_test_angle))/deg2rad)-90.0;
4400 elevation=((acos(cos_xmtr_angle))/deg2rad)-90.0;
4402 /* Integrate the antenna's radiation
4403 pattern into the overall path loss. */
4405 x=(int)rint(10.0*(10.0-elevation));
4407 if (x>=0 && x<=1000)
4409 pattern=(double)LR.antenna_pattern[(int)azimuth][x];
4412 patterndB=20.0*log10(pattern);
4418 total_loss=loss-patterndB;
4422 fprintf(fd,"%f\t%f\n",KM_PER_MILE*(path.distance[path.length-1]-path.distance[y]),total_loss);
4423 fprintf(fd2,"%7.2f\t\t%7.2f\t\t %d\t%s\n",KM_PER_MILE*path.distance[y],total_loss, errnum, strmode);
4428 fprintf(fd,"%f\t%f\n",path.distance[path.length-1]-path.distance[y],total_loss);
4429 fprintf(fd2,"%7.2f\t\t%7.2f\t\t %d\t%s\n",path.distance[y],total_loss, errnum, strmode);
4434 if (total_loss>maxloss)
4437 if (total_loss<minloss)
4443 fprintf(fd2,"\nLongley-Rice path loss between %s and %s is %.2f dB.\n",source.name, destination.name, loss);
4446 fprintf(fd2,"Total path loss including TX antenna pattern is %.2f dB.\n",total_loss);
4450 fprintf(fd2,"\nNotes on \"errnum\":\n\n");
4451 fprintf(fd2," 0: No error. :-)\n");
4452 fprintf(fd2," 1: Warning! Some parameters are nearly out of range.\n");
4453 fprintf(fd2," Results should be used with caution.\n");
4454 fprintf(fd2," 2: Note: Default parameters have been substituted for impossible ones.\n");
4455 fprintf(fd2," 3: Warning! A combination of parameters is out of range.\n");
4456 fprintf(fd2," Results are probably invalid.\n");
4457 fprintf(fd2," Other: Warning! Some parameters are out of range.\n");
4458 fprintf(fd2," Results are probably invalid.\n\nEnd of Report\n");
4461 fprintf(stdout,"Longley-Rice path loss between %s and %s is %.2f dB.\n",source.name, destination.name, loss);
4464 fprintf(stdout,"Total path loss including TX antenna pattern is %.2f dB.\n",total_loss);
4466 fprintf(stdout,"Path Loss Report written to: \"%s\"\n",report_name);
4473 /* Default filename and output file type */
4475 strncpy(filename,"loss\0",5);
4476 strncpy(term,"png\0",4);
4477 strncpy(ext,"png\0",4);
4482 /* Grab extension and terminal type from "name" */
4484 for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
4485 filename[x]=name[x];
4489 for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
4491 term[y]=tolower(name[x]);
4501 { /* No extension -- Default is png */
4504 strncpy(term,"png\0",4);
4505 strncpy(ext,"png\0",4);
4509 /* Either .ps or .postscript may be used
4510 as an extension for postscript output. */
4512 if (strncmp(term,"postscript",10)==0)
4513 strncpy(ext,"ps\0",3);
4515 else if (strncmp(ext,"ps",2)==0)
4516 strncpy(term,"postscript enhanced color\0",26);
4518 fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
4521 fd=fopen("splat.gp","w");
4523 fprintf(fd,"set grid\n");
4524 fprintf(fd,"set yrange [%2.3f to %2.3f]\n", minloss, maxloss);
4525 fprintf(fd,"set encoding iso_8859_1\n");
4526 fprintf(fd,"set term %s\n",term);
4527 fprintf(fd,"set title \"SPLAT! Loss Profile Along Path Between %s and %s (%.2f%c azimuth)\"\n",destination.name, source.name, Azimuth(destination,source),176);
4530 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f kilometers)\"\n",destination.name,source.name,KM_PER_MILE*Distance(destination,source));
4532 fprintf(fd,"set xlabel \"Distance Between %s and %s (%.2f miles)\"\n",destination.name,source.name,Distance(destination,source));
4534 if (got_azimuth_pattern || got_elevation_pattern)
4535 fprintf(fd,"set ylabel \"Total Path Loss (including TX antenna pattern) (dB)");
4537 fprintf(fd,"set ylabel \"Longley-Rice Path Loss (dB)");
4539 fprintf(fd,"\"\nset output \"%s.%s\"\n",filename,ext);
4540 fprintf(fd,"plot \"profile.gp\" title \"Path Loss\" with lines\n");
4544 x=system("gnuplot splat.gp");
4549 unlink("profile.gp");
4550 unlink("reference.gp");
4552 fprintf(stdout," Done!\n");
4557 fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
4560 void ObstructionReport(struct site xmtr, struct site rcvr, char report, double f)
4564 double h_r, h_t, h_x, h_r_orig, cos_tx_angle, cos_test_angle,
4565 cos_tx_angle_f1, cos_tx_angle_fpt6, haavt, d_tx, d_x,
4566 h_r_f1, h_r_fpt6, h_f, h_los, lambda=0.0, azimuth,
4567 pattern, patterndB, distance, angle1, angle2;
4568 char report_name[80], string[255], string_fpt6[255],
4572 sprintf(report_name,"%s-to-%s.txt",xmtr.name,rcvr.name);
4574 for (x=0; report_name[x]!=0; x++)
4575 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
4578 fd=fopen(report_name,"w");
4580 fprintf(fd,"\n\t\t--==[ SPLAT! v%s Obstruction Report ]==--\n\n",splat_version);
4581 fprintf(fd,"Analysis of great circle path between %s and %s:\n",xmtr.name, rcvr.name);
4582 fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
4583 fprintf(fd,"Transmitter site: %s\n",xmtr.name);
4588 fprintf(fd,"Site location: %.4f North / %.4f West",xmtr.lat, xmtr.lon);
4589 fprintf(fd, " (%s N / ", dec2dms(xmtr.lat));
4594 fprintf(fd,"Site location: %.4f South / %.4f West",-xmtr.lat, xmtr.lon);
4595 fprintf(fd, " (%s S / ", dec2dms(xmtr.lat));
4598 fprintf(fd, "%s W)\n", dec2dms(xmtr.lon));
4602 fprintf(fd,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(xmtr));
4603 fprintf(fd,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*xmtr.alt, METERS_PER_FOOT*(xmtr.alt+GetElevation(xmtr)));
4608 fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(xmtr));
4609 fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",xmtr.alt, xmtr.alt+GetElevation(xmtr));
4617 fprintf(fd,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4619 fprintf(fd,"Antenna height above average terrain: %.2f feet\n",haavt);
4624 distance=Distance(xmtr,rcvr);
4625 azimuth=Azimuth(xmtr,rcvr);
4626 angle1=ElevationAngle(xmtr,rcvr);
4627 angle2=ElevationAngle2(xmtr,rcvr,earthradius);
4629 if (got_azimuth_pattern || got_elevation_pattern)
4631 x=(int)rint(10.0*(10.0-angle2));
4633 if (x>=0 && x<=1000)
4634 pattern=(double)LR.antenna_pattern[(int)rint(azimuth)][x];
4638 fprintf(fd,"Antenna pattern toward %s: %.3f",rcvr.name,pattern);
4639 patterndB=20.0*log10(pattern);
4640 fprintf(fd," (%.2f dB)\n",patterndB);
4645 fprintf(fd,"Distance to %s: %.2f kilometers\n",rcvr.name,KM_PER_MILE*distance);
4648 fprintf(fd,"Distance to %s: %.2f miles\n",rcvr.name,distance);
4650 fprintf(fd,"Azimuth to %s: %.2f degrees\n",rcvr.name,azimuth);
4653 fprintf(fd,"Elevation angle to %s: %+.4f degrees\n",rcvr.name,angle1);
4656 fprintf(fd,"Depression angle to %s: %+.4f degrees\n",rcvr.name,angle1);
4661 fprintf(fd,"Depression");
4663 fprintf(fd,"Elevation");
4665 fprintf(fd," angle to the first obstruction: %+.4f degrees\n",angle2);
4668 fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
4672 fprintf(fd,"Receiver site: %s\n",rcvr.name);
4676 fprintf(fd,"Site location: %.4f North / %.4f West",rcvr.lat, rcvr.lon);
4677 fprintf(fd, " (%s N / ", dec2dms(rcvr.lat));
4682 fprintf(fd,"Site location: %.4f South / %.4f West",-rcvr.lat, rcvr.lon);
4683 fprintf(fd, " (%s S / ", dec2dms(rcvr.lat));
4686 fprintf(fd, "%s W)\n", dec2dms(rcvr.lon));
4690 fprintf(fd,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(rcvr));
4691 fprintf(fd,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*rcvr.alt, METERS_PER_FOOT*(rcvr.alt+GetElevation(rcvr)));
4696 fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(rcvr));
4697 fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",rcvr.alt, rcvr.alt+GetElevation(rcvr));
4705 fprintf(fd,"Antenna height above average terrain: %.2f meters\n",METERS_PER_FOOT*haavt);
4707 fprintf(fd,"Antenna height above average terrain: %.2f feet\n",haavt);
4710 azimuth=Azimuth(rcvr,xmtr);
4711 angle1=ElevationAngle(rcvr,xmtr);
4712 angle2=ElevationAngle2(rcvr,xmtr,earthradius);
4715 fprintf(fd,"Distance to %s: %.2f kilometers\n",xmtr.name,KM_PER_MILE*distance);
4717 fprintf(fd,"Distance to %s: %.2f miles\n",xmtr.name,distance);
4719 fprintf(fd,"Azimuth to %s: %.2f degrees\n",xmtr.name,azimuth);
4722 fprintf(fd,"Elevation to %s: %+.4f degrees\n",xmtr.name,angle1);
4725 fprintf(fd,"Depression angle to %s: %+.4f degrees\n",xmtr.name,angle1);
4730 fprintf(fd,"Depression");
4732 fprintf(fd,"Elevation");
4734 fprintf(fd," angle to the first obstruction: %+.4f degrees\n",angle2);
4738 fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
4742 /* Generate profile of the terrain. Create the path
4743 from transmitter to receiver because that's the
4744 way the original los() function did it, and going
4745 the other way can yield slightly different results. */
4747 ReadPath(xmtr,rcvr);
4748 h_r=GetElevation(rcvr)+rcvr.alt+earthradius;
4752 h_t=GetElevation(xmtr)+xmtr.alt+earthradius;
4753 d_tx=5280.0*Distance(rcvr,xmtr);
4754 cos_tx_angle=((h_r*h_r)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r*d_tx);
4755 cos_tx_angle_f1=cos_tx_angle;
4756 cos_tx_angle_fpt6=cos_tx_angle;
4759 lambda=9.8425e8/(f*1e6);
4761 /* At each point along the path calculate the cosine
4762 of a sort of "inverse elevation angle" at the receiver.
4763 From the antenna, 0 deg. looks at the ground, and 90 deg.
4764 is parallel to the ground.
4766 Start at the receiver. If this is the lowest antenna,
4767 then terrain obstructions will be nearest to it. (Plus,
4768 that's the way the original los() did it.)
4770 Calculate cosines only. That's sufficient to compare
4771 angles and it saves the extra computational burden of
4772 acos(). However, note the inverted comparison: if
4773 acos(A) > acos(B), then B > A. */
4775 for (x=path.length-1; x>0; x--)
4777 site_x.lat=path.lat[x];
4778 site_x.lon=path.lon[x];
4781 h_x=GetElevation(site_x)+earthradius;
4782 d_x=5280.0*Distance(rcvr,site_x);
4784 /* Deal with the LOS path first. */
4786 cos_test_angle=((h_r*h_r)+(d_x*d_x)-(h_x*h_x))/(2.0*h_r*d_x);
4788 if (cos_tx_angle>cos_test_angle)
4791 fprintf(fd,"SPLAT! detected obstructions at:\n\n");
4793 if (site_x.lat>=0.0)
4796 fprintf(fd,"\t%.4f N, %.4f W, %5.2f kilometers, %6.2f meters AMSL\n",site_x.lat, site_x.lon, KM_PER_MILE*(d_x/5280.0), METERS_PER_FOOT*(h_x-earthradius));
4798 fprintf(fd,"\t%.4f N, %.4f W, %5.2f miles, %6.2f feet AMSL\n",site_x.lat, site_x.lon, d_x/5280.0, h_x-earthradius);
4804 fprintf(fd,"\t%.4f S, %.4f W, %5.2f kilometers, %6.2f meters AMSL\n",-site_x.lat, site_x.lon, KM_PER_MILE*(d_x/5280.0), METERS_PER_FOOT*(h_x-earthradius));
4807 fprintf(fd,"\t%.4f S, %.4f W, %5.2f miles, %6.2f feet AMSL\n",-site_x.lat, site_x.lon, d_x/5280.0, h_x-earthradius);
4811 while (cos_tx_angle>cos_test_angle)
4814 cos_test_angle=((h_r*h_r)+(d_x*d_x)-(h_x*h_x))/(2.0*h_r*d_x);
4815 cos_tx_angle=((h_r*h_r)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r*d_tx);
4820 /* Now clear the first Fresnel zone, but don't
4821 clutter the obstruction report. */
4823 cos_tx_angle_f1=((h_r_f1*h_r_f1)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r_f1*d_tx);
4824 h_los=sqrt(h_r_f1*h_r_f1+d_x*d_x-2*h_r_f1*d_x*cos_tx_angle_f1);
4825 h_f=h_los-sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4830 cos_tx_angle_f1=((h_r_f1*h_r_f1)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r_f1*d_tx);
4831 h_los=sqrt(h_r_f1*h_r_f1+d_x*d_x-2*h_r_f1*d_x*cos_tx_angle_f1);
4832 h_f=h_los-sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4835 /* And clear the 60% F1 zone. */
4837 cos_tx_angle_fpt6=((h_r_fpt6*h_r_fpt6)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r_fpt6*d_tx);
4838 h_los=sqrt(h_r_fpt6*h_r_fpt6+d_x*d_x-2*h_r_fpt6*d_x*cos_tx_angle_fpt6);
4839 h_f=h_los-0.6*sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4844 cos_tx_angle_fpt6=((h_r_fpt6*h_r_fpt6)+(d_tx*d_tx)-(h_t*h_t))/(2.0*h_r_fpt6*d_tx);
4845 h_los=sqrt(h_r_fpt6*h_r_fpt6+d_x*d_x-2*h_r_fpt6*d_x*cos_tx_angle_fpt6);
4846 h_f=h_los-0.6*sqrt(lambda*d_x*(d_tx-d_x)/d_tx);
4854 sprintf(string,"\nAntenna at %s must be raised to at least %.2f meters AGL\nto clear all obstructions detected by SPLAT!\n",rcvr.name, METERS_PER_FOOT*(h_r-GetElevation(rcvr)-earthradius));
4856 sprintf(string,"\nAntenna at %s must be raised to at least %.2f feet AGL\nto clear all obstructions detected by SPLAT!\n",rcvr.name, h_r-GetElevation(rcvr)-earthradius);
4860 sprintf(string,"\nNo obstructions to LOS path due to terrain were detected by SPLAT!\n");
4864 if (h_r_fpt6>h_r_orig)
4867 sprintf(string_fpt6,"\nAntenna at %s must be raised to at least %.2f meters AGL\nto clear 60%c of the first Fresnel zone.\n",rcvr.name, METERS_PER_FOOT*(h_r_fpt6-GetElevation(rcvr)-earthradius),37);
4870 sprintf(string_fpt6,"\nAntenna at %s must be raised to at least %.2f feet AGL\nto clear 60%c of the first Fresnel zone.\n",rcvr.name, h_r_fpt6-GetElevation(rcvr)-earthradius,37);
4874 sprintf(string_fpt6,"\n60%c of the first Fresnel zone is clear.\n",37);
4876 if (h_r_f1>h_r_orig)
4879 sprintf(string_f1,"\nAntenna at %s must be raised to at least %.2f meters AGL\nto clear the first Fresnel zone.\n",rcvr.name, METERS_PER_FOOT*(h_r_f1-GetElevation(rcvr)-earthradius));
4882 sprintf(string_f1,"\nAntenna at %s must be raised to at least %.2f feet AGL\nto clear the first Fresnel zone.\n",rcvr.name, h_r_f1-GetElevation(rcvr)-earthradius);
4887 sprintf(string_f1,"\nThe first Fresnel zone is clear.\n\n");
4891 fprintf(fd,"%s",string);
4895 fprintf(fd,"%s",string_f1);
4896 fprintf(fd,"%s",string_fpt6);
4901 /* Display report summary on terminal */
4903 /* Line-of-sight status */
4905 fprintf(stdout,"%s",string);
4909 /* Fresnel zone status */
4911 fprintf(stdout,"%s",string_f1);
4912 fprintf(stdout,"%s",string_fpt6);
4915 fprintf(stdout, "\nObstruction report written to: \"%s\"\n",report_name);
4920 void SiteReport(struct site xmtr)
4922 char report_name[80];
4927 sprintf(report_name,"%s-site_report.txt",xmtr.name);
4929 for (x=0; report_name[x]!=0; x++)
4930 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
4933 fd=fopen(report_name,"w");
4935 fprintf(fd,"\n\t--==[ SPLAT! v%s Site Analysis Report For: %s ]==--\n\n",splat_version,xmtr.name);
4937 fprintf(fd,"---------------------------------------------------------------------------\n\n");
4941 fprintf(fd,"Site location: %.4f North / %.4f West",xmtr.lat, xmtr.lon);
4942 fprintf(fd, " (%s N / ",dec2dms(xmtr.lat));
4947 fprintf(fd,"Site location: %.4f South / %.4f West",-xmtr.lat, xmtr.lon);
4948 fprintf(fd, " (%s S / ",dec2dms(xmtr.lat));
4951 fprintf(fd, "%s W)\n",dec2dms(xmtr.lon));
4955 fprintf(fd,"Ground elevation: %.2f meters AMSL\n",METERS_PER_FOOT*GetElevation(xmtr));
4956 fprintf(fd,"Antenna height: %.2f meters AGL / %.2f meters AMSL\n",METERS_PER_FOOT*xmtr.alt, METERS_PER_FOOT*(xmtr.alt+GetElevation(xmtr)));
4961 fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(xmtr));
4962 fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",xmtr.alt, xmtr.alt+GetElevation(xmtr));
4967 if (terrain>-4999.0)
4970 fprintf(fd,"Antenna height above average terrain: %.2f meters\n\n",METERS_PER_FOOT*terrain);
4972 fprintf(fd,"Antenna height above average terrain: %.2f feet\n\n",terrain);
4974 /* Display the average terrain between 2 and 10 miles
4975 from the transmitter site at azimuths of 0, 45, 90,
4976 135, 180, 225, 270, and 315 degrees. */
4978 for (azi=0; azi<=315; azi+=45)
4980 fprintf(fd,"Average terrain at %3d degrees azimuth: ",azi);
4981 terrain=AverageTerrain(xmtr,(double)azi,2.0,10.0);
4983 if (terrain>-4999.0)
4986 fprintf(fd,"%.2f meters AMSL\n",METERS_PER_FOOT*terrain);
4988 fprintf(fd,"%.2f feet AMSL\n",terrain);
4992 fprintf(fd,"No terrain\n");
4996 fprintf(fd,"\n---------------------------------------------------------------------------\n\n");
4998 fprintf(stdout,"\nSite analysis report written to: \"%s\"\n",report_name);
5001 void LoadTopoData(int max_lon, int min_lon, int max_lat, int min_lat)
5003 /* This function loads the SDF files required
5004 to cover the limits of the region specified. */
5006 int x, y, width, ymin, ymax;
5008 width=ReduceAngle(max_lon-min_lon);
5010 if ((max_lon-min_lon)<=180.0)
5012 for (y=0; y<=width; y++)
5013 for (x=min_lat; x<=max_lat; x++)
5015 ymin=(int)(min_lon+(double)y);
5031 sprintf(string,"%d:%d:%d:%d",x, x+1, ymin, ymax);
5038 for (y=0; y<=width; y++)
5039 for (x=min_lat; x<=max_lat; x++)
5057 sprintf(string,"%d:%d:%d:%d",x, x+1, ymin, ymax);
5063 int LoadPLI(char *filename)
5065 int error=0, max_west, min_west, max_north, min_north;
5066 char string[80], *pointer=NULL;
5067 double latitude=0.0, longitude=0.0, azimuth=0.0, elevation=0.0,
5071 fd=fopen(filename,"r");
5075 fgets(string,78,fd);
5076 pointer=strchr(string,';');
5081 sscanf(string,"%d, %d",&max_west, &min_west);
5083 fgets(string,78,fd);
5084 pointer=strchr(string,';');
5089 sscanf(string,"%d, %d",&max_north, &min_north);
5091 fgets(string,78,fd);
5092 pointer=strchr(string,';');
5097 LoadTopoData(max_west-1, min_west, max_north-1, min_north);
5099 fprintf(stdout,"\nReading \"%s\"... ",filename);
5102 fscanf(fd,"%lf, %lf, %lf, %lf, %lf",&latitude, &longitude, &azimuth, &elevation, &loss);
5116 if (loss<=(double)maxdB)
5117 OrMask(latitude,longitude,((unsigned char)(loss))<<3);
5119 fscanf(fd,"%lf, %lf, %lf, %lf, %lf",&latitude, &longitude, &azimuth, &elevation, &loss);
5124 fprintf(stdout," Done!\n");
5134 void WriteKML(struct site source, struct site destination)
5137 char block, report_name[80];
5138 double distance, rx_alt, tx_alt, cos_xmtr_angle,
5139 azimuth, cos_test_angle, test_alt;
5142 ReadPath(source,destination);
5144 sprintf(report_name,"%s-to-%s.kml",source.name,destination.name);
5146 for (x=0; report_name[x]!=0; x++)
5147 if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
5150 fd=fopen(report_name,"w");
5152 fprintf(fd,"<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n");
5153 fprintf(fd,"<kml xmlns=\"http://earth.google.com/kml/2.0\">\n");
5154 fprintf(fd,"<!-- Generated by SPLAT! Version %s -->\n",splat_version);
5155 fprintf(fd,"<Folder>\n");
5156 fprintf(fd,"<name>SPLAT! Path</name>\n");
5157 fprintf(fd,"<open>1</open>\n");
5158 fprintf(fd,"<description>Path Between %s and %s</description>\n",source.name,destination.name);
5160 fprintf(fd,"<Placemark>\n");
5161 fprintf(fd," <name>%s</name>\n",source.name);
5162 fprintf(fd," <description>\n");
5163 fprintf(fd," Transmit Site\n");
5165 if (source.lat>=0.0)
5166 fprintf(fd," <BR>%s North</BR>\n",dec2dms(source.lat));
5168 fprintf(fd," <BR>%s South</BR>\n",dec2dms(source.lat));
5170 fprintf(fd," <BR>%s West</BR>\n",dec2dms(source.lon));
5172 azimuth=Azimuth(source,destination);
5173 distance=Distance(source,destination);
5176 fprintf(fd," <BR>%.2f km",distance*KM_PER_MILE);
5178 fprintf(fd," <BR>%.2f miles",distance);
5180 fprintf(fd," to %s</BR>\n <BR>toward an azimuth of %.2f%c</BR>\n",destination.name,azimuth,176);
5182 fprintf(fd," </description>\n");
5183 fprintf(fd," <visibility>1</visibility>\n");
5184 fprintf(fd," <Style>\n");
5185 fprintf(fd," <IconStyle>\n");
5186 fprintf(fd," <Icon>\n");
5187 fprintf(fd," <href>root://icons/palette-5.png</href>\n");
5188 fprintf(fd," <x>224</x>\n");
5189 fprintf(fd," <y>224</y>\n");
5190 fprintf(fd," <w>32</w>\n");
5191 fprintf(fd," <h>32</h>\n");
5192 fprintf(fd," </Icon>\n");
5193 fprintf(fd," </IconStyle>\n");
5194 fprintf(fd," </Style>\n");
5195 fprintf(fd," <Point>\n");
5196 fprintf(fd," <extrude>1</extrude>\n");
5197 fprintf(fd," <altitudeMode>relativeToGround</altitudeMode>\n");
5198 fprintf(fd," <coordinates>%f,%f,30</coordinates>\n",(source.lon<180.0?-source.lon:360.0-source.lon),source.lat);
5199 fprintf(fd," </Point>\n");
5200 fprintf(fd,"</Placemark>\n");
5202 fprintf(fd,"<Placemark>\n");
5203 fprintf(fd," <name>%s</name>\n",destination.name);
5204 fprintf(fd," <description>\n");
5205 fprintf(fd," Receive Site\n");
5207 if (destination.lat>=0.0)
5208 fprintf(fd," <BR>%s North</BR>\n",dec2dms(destination.lat));
5210 fprintf(fd," <BR>%s South</BR>\n",dec2dms(destination.lat));
5212 fprintf(fd," <BR>%s West</BR>\n",dec2dms(destination.lon));
5216 fprintf(fd," <BR>%.2f km",distance*KM_PER_MILE);
5218 fprintf(fd," <BR>%.2f miles",distance);
5220 fprintf(fd," to %s</BR>\n <BR>toward an azimuth of %.2f%c</BR>\n",source.name,Azimuth(destination,source),176);
5222 fprintf(fd," </description>\n");
5223 fprintf(fd," <visibility>1</visibility>\n");
5224 fprintf(fd," <Style>\n");
5225 fprintf(fd," <IconStyle>\n");
5226 fprintf(fd," <Icon>\n");
5227 fprintf(fd," <href>root://icons/palette-5.png</href>\n");
5228 fprintf(fd," <x>224</x>\n");
5229 fprintf(fd," <y>224</y>\n");
5230 fprintf(fd," <w>32</w>\n");
5231 fprintf(fd," <h>32</h>\n");
5232 fprintf(fd," </Icon>\n");
5233 fprintf(fd," </IconStyle>\n");
5234 fprintf(fd," </Style>\n");
5235 fprintf(fd," <Point>\n");
5236 fprintf(fd," <extrude>1</extrude>\n");
5237 fprintf(fd," <altitudeMode>relativeToGround</altitudeMode>\n");
5238 fprintf(fd," <coordinates>%f,%f,30</coordinates>\n",(destination.lon<180.0?-destination.lon:360.0-destination.lon),destination.lat);
5239 fprintf(fd," </Point>\n");
5240 fprintf(fd,"</Placemark>\n");
5242 fprintf(fd,"<Placemark>\n");
5243 fprintf(fd,"<name>Point-to-Point Path</name>\n");
5244 fprintf(fd," <visibility>1</visibility>\n");
5245 fprintf(fd," <open>0</open>\n");
5246 fprintf(fd," <Style>\n");
5247 fprintf(fd," <LineStyle>\n");
5248 fprintf(fd," <color>7fffffff</color>\n");
5249 fprintf(fd," </LineStyle>\n");
5250 fprintf(fd," <PolyStyle>\n");
5251 fprintf(fd," <color>7fffffff</color>\n");
5252 fprintf(fd," </PolyStyle>\n");
5253 fprintf(fd," </Style>\n");
5254 fprintf(fd," <LineString>\n");
5255 fprintf(fd," <extrude>1</extrude>\n");
5256 fprintf(fd," <tessellate>1</tessellate>\n");
5257 fprintf(fd," <altitudeMode>relativeToGround</altitudeMode>\n");
5258 fprintf(fd," <coordinates>\n");
5260 for (x=0; x<path.length; x++)
5261 fprintf(fd," %f,%f,5\n",(path.lon[x]<180.0?-path.lon[x]:360.0-path.lon[x]),path.lat[x]);
5263 fprintf(fd," </coordinates>\n");
5264 fprintf(fd," </LineString>\n");
5265 fprintf(fd,"</Placemark>\n");
5267 fprintf(fd,"<Placemark>\n");
5268 fprintf(fd,"<name>Line-of-Sight Path</name>\n");
5269 fprintf(fd," <visibility>1</visibility>\n");
5270 fprintf(fd," <open>0</open>\n");
5271 fprintf(fd," <Style>\n");
5272 fprintf(fd," <LineStyle>\n");
5273 fprintf(fd," <color>ff00ff00</color>\n");
5274 fprintf(fd," </LineStyle>\n");
5275 fprintf(fd," <PolyStyle>\n");
5276 fprintf(fd," <color>7f00ff00</color>\n");
5277 fprintf(fd," </PolyStyle>\n");
5278 fprintf(fd," </Style>\n");
5279 fprintf(fd," <LineString>\n");
5280 fprintf(fd," <extrude>1</extrude>\n");
5281 fprintf(fd," <tessellate>1</tessellate>\n");
5282 fprintf(fd," <altitudeMode>relativeToGround</altitudeMode>\n");
5283 fprintf(fd," <coordinates>\n");
5285 /* Walk across the "path", indentifying obstructions along the way */
5287 for (y=0; y<path.length; y++)
5289 distance=5280.0*path.distance[y];
5290 tx_alt=earthradius+source.alt+path.elevation[0];
5291 rx_alt=earthradius+destination.alt+path.elevation[y];
5293 /* Calculate the cosine of the elevation of the
5294 transmitter as seen at the temp rx point. */
5296 cos_xmtr_angle=((rx_alt*rx_alt)+(distance*distance)-(tx_alt*tx_alt))/(2.0*rx_alt*distance);
5298 for (x=y, block=0; x>=0 && block==0; x--)
5300 distance=5280.0*(path.distance[y]-path.distance[x]);
5301 test_alt=earthradius+path.elevation[x];
5303 cos_test_angle=((rx_alt*rx_alt)+(distance*distance)-(test_alt*test_alt))/(2.0*rx_alt*distance);
5305 /* Compare these two angles to determine if
5306 an obstruction exists. Since we're comparing
5307 the cosines of these angles rather than
5308 the angles themselves, the following "if"
5309 statement is reversed from what it would
5310 be if the actual angles were compared. */
5312 if (cos_xmtr_angle>cos_test_angle)
5317 fprintf(fd," %f,%f,-30\n",(path.lon[y]<180.0?-path.lon[y]:360.0-path.lon[y]),path.lat[y]);
5319 fprintf(fd," %f,%f,5\n",(path.lon[y]<180.0?-path.lon[y]:360.0-path.lon[y]),path.lat[y]);
5322 fprintf(fd," </coordinates>\n");
5323 fprintf(fd," </LineString>\n");
5324 fprintf(fd,"</Placemark>\n");
5326 fprintf(fd," <LookAt>\n");
5327 fprintf(fd," <longitude>%f</longitude>\n",(source.lon<180.0?-source.lon:360.0-source.lon));
5328 fprintf(fd," <latitude>%f</latitude>\n",source.lat);
5329 fprintf(fd," <range>300.0</range>\n");
5330 fprintf(fd," <tilt>45.0</tilt>\n");
5331 fprintf(fd," <heading>%f</heading>\n",azimuth);
5332 fprintf(fd," </LookAt>\n");
5334 fprintf(fd,"</Folder>\n");
5335 fprintf(fd,"</kml>\n");
5339 fprintf(stdout, "KML file written to: \"%s\"\n",report_name);
5344 int main(int argc, char *argv[])
5346 int x, y, z=0, min_lat, min_lon, max_lat, max_lon,
5347 rxlat, rxlon, txlat, txlon, west_min, west_max,
5348 north_min, north_max;
5350 unsigned char coverage=0, LRmap=0, ext[20], terrain_plot=0,
5351 elevation_plot=0, height_plot=0,
5352 longley_plot=0, cities=0, bfs=0, txsites=0,
5353 count, report='y', norm=0, topomap=0, geo=0,
5356 char mapfile[255], header[80], city_file[5][255],
5357 elevation_file[255], height_file[255],
5358 longley_file[255], terrain_file[255],
5359 string[255], rxfile[255], *env=NULL,
5360 txfile[255], map=0, boundary_file[5][255],
5361 udt_file[255], rxsite=0, plo_filename[255],
5362 pli_filename[255], nf=0;
5364 double altitude=0.0, altitudeLR=0.0, tx_range=0.0,
5365 rx_range=0.0, deg_range=0.0, deg_limit,
5366 deg_range_lon, er_mult, freq=0.0;
5368 struct site tx_site[4], rx_site;
5375 fprintf(stdout,"\n\t\t --==[ SPLAT! v%s Available Options... ]==--\n\n",splat_version);
5376 fprintf(stdout," -t txsite(s).qth (max of 4)\n");
5377 fprintf(stdout," -r rxsite.qth\n");
5378 fprintf(stdout," -c plot coverage of TX(s) with an RX antenna at X feet/meters AGL\n");
5379 fprintf(stdout," -L plot path loss map of TX based on an RX at X feet/meters AGL\n");
5380 fprintf(stdout," -s filename(s) of city/site file(s) to import (5 max)\n");
5381 fprintf(stdout," -b filename(s) of cartographic boundary file(s) to import (max of 5)\n");
5382 fprintf(stdout," -p filename of terrain profile graph to plot\n");
5383 fprintf(stdout," -e filename of terrain elevation graph to plot\n");
5384 fprintf(stdout," -h filename of terrain height graph to plot\n");
5385 fprintf(stdout," -H filename of normalized terrain height graph to plot\n");
5386 fprintf(stdout," -l filename of Longley-Rice graph to plot\n");
5387 fprintf(stdout," -o filename of topographic map to generate (.ppm)\n");
5388 fprintf(stdout," -u filename of user-defined terrain file to import\n");
5389 fprintf(stdout," -d sdf file directory path (overrides path in ~/.splat_path file)\n");
5390 fprintf(stdout," -n no analysis, brief report\n");
5391 fprintf(stdout," -N no analysis, no report\n");
5392 fprintf(stdout," -m earth radius multiplier\n");
5393 fprintf(stdout," -f frequency for Fresnel zone calculation (MHz)\n");
5394 fprintf(stdout," -R modify default range for -c or -L (miles/kilometers)\n");
5395 fprintf(stdout," -db maximum loss contour to display on path loss maps (80-230 dB)\n");
5396 fprintf(stdout," -nf do not plot Fresnel zones in height plots\n");
5397 fprintf(stdout," -plo filename of path-loss output file\n");
5398 fprintf(stdout," -pli filename of path-loss input file\n");
5399 fprintf(stdout," -udt filename of user defined terrain input file\n");
5400 fprintf(stdout," -geo generate an Xastir .geo georeference file (with .ppm output)\n");
5401 fprintf(stdout," -kml generate a Google Earth .kml file (for point-to-point links)\n");
5402 fprintf(stdout," -metric employ metric rather than imperial units for all user I/O\n\n");
5404 fprintf(stdout,"If that flew by too fast, consider piping the output through 'less':\n");
5405 fprintf(stdout,"\n\tsplat | less\n\n");
5406 fprintf(stdout,"Type 'man splat', or see the documentation for more details.\n\n");
5418 elevation_file[0]=0;
5428 earthradius=EARTHRADIUS;
5430 sprintf(header,"\n\t\t--==[ Welcome To SPLAT! v%s ]==--\n\n", splat_version);
5434 tx_site[x].lat=91.0;
5435 tx_site[x].lon=361.0;
5438 for (x=0; x<MAXSLOTS; x++)
5440 dem[x].min_el=32768;
5441 dem[x].max_el=-32768;
5442 dem[x].min_north=90;
5443 dem[x].max_north=-90;
5444 dem[x].min_west=360;
5448 /* Scan for command line arguments */
5450 for (x=1; x<=y; x++)
5452 if (strcmp(argv[x],"-R")==0)
5456 if (z<=y && argv[z][0] && argv[z][0]!='-')
5458 sscanf(argv[z],"%lf",&max_range);
5463 if (max_range>1000.0)
5468 if (strcmp(argv[x],"-m")==0)
5472 if (z<=y && argv[z][0] && argv[z][0]!='-')
5474 sscanf(argv[z],"%lf",&er_mult);
5482 earthradius*=er_mult;
5486 if (strcmp(argv[x],"-o")==0)
5490 if (z<=y && argv[z][0] && argv[z][0]!='-')
5491 strncpy(mapfile,argv[z],253);
5495 if (strcmp(argv[x],"-u")==0)
5499 if (z<=y && argv[z][0] && argv[z][0]!='-')
5500 strncpy(udt_file,argv[z],253);
5503 if (strcmp(argv[x],"-c")==0)
5507 if (z<=y && argv[z][0] && argv[z][0]!='-')
5509 sscanf(argv[z],"%lf",&altitude);
5514 if (strcmp(argv[x],"-db")==0 || strcmp(argv[x],"-dB")==0)
5518 if (z<=y && argv[z][0] && argv[z][0]!='-')
5520 sscanf(argv[z],"%d",&maxdB);
5532 if (strcmp(argv[x],"-p")==0)
5536 if (z<=y && argv[z][0] && argv[z][0]!='-')
5538 strncpy(terrain_file,argv[z],253);
5543 if (strcmp(argv[x],"-e")==0)
5547 if (z<=y && argv[z][0] && argv[z][0]!='-')
5549 strncpy(elevation_file,argv[z],253);
5554 if (strcmp(argv[x],"-h")==0 || strcmp(argv[x],"-H")==0)
5558 if (z<=y && argv[z][0] && argv[z][0]!='-')
5560 strncpy(height_file,argv[z],253);
5564 if (strcmp(argv[x],"-H")==0)
5570 if (strcmp(argv[x],"-n")==0)
5576 if (strcmp(argv[x],"-N")==0)
5582 if (strcmp(argv[x],"-metric")==0)
5585 if (strcmp(argv[x],"-geo")==0)
5588 if (strcmp(argv[x],"-kml")==0)
5591 if (strcmp(argv[x],"-nf")==0)
5594 if (strcmp(argv[x],"-d")==0)
5598 if (z<=y && argv[z][0] && argv[z][0]!='-')
5599 strncpy(sdf_path,argv[z],253);
5602 if (strcmp(argv[x],"-t")==0)
5604 /* Read Transmitter Location */
5608 while (z<=y && argv[z][0] && argv[z][0]!='-' && txsites<4)
5610 strncpy(txfile,argv[z],253);
5611 tx_site[txsites]=LoadQTH(txfile);
5619 if (strcmp(argv[x],"-L")==0)
5623 if (z<=y && argv[z][0] && argv[z][0]!='-')
5625 sscanf(argv[z],"%lf",&altitudeLR);
5628 fprintf(stdout,"c and L are exclusive options, ignoring L.\n");
5637 if (strcmp(argv[x],"-l")==0)
5641 if (z<=y && argv[z][0] && argv[z][0]!='-')
5643 strncpy(longley_file,argv[z],253);
5645 /* Doing this twice is harmless */
5650 if (strcmp(argv[x],"-r")==0)
5652 /* Read Receiver Location */
5656 if (z<=y && argv[z][0] && argv[z][0]!='-')
5658 strncpy(rxfile,argv[z],253);
5659 rx_site=LoadQTH(rxfile);
5664 if (strcmp(argv[x],"-s")==0)
5666 /* Read city file(s) */
5670 while (z<=y && argv[z][0] && argv[z][0]!='-' && cities<5)
5672 strncpy(city_file[cities],argv[z],253);
5680 if (strcmp(argv[x],"-b")==0)
5682 /* Read Boundary File(s) */
5686 while (z<=y && argv[z][0] && argv[z][0]!='-' && bfs<5)
5688 strncpy(boundary_file[bfs],argv[z],253);
5696 if (strcmp(argv[x],"-f")==0)
5700 if (z<=y && argv[z][0] && argv[z][0]!='-')
5702 sscanf(argv[z],"%lf",&freq);
5712 if (strcmp(argv[x],"-plo")==0)
5716 if (z<=y && argv[z][0] && argv[z][0]!='-')
5717 strncpy(plo_filename,argv[z],253);
5720 if (strcmp(argv[x],"-pli")==0)
5724 if (z<=y && argv[z][0] && argv[z][0]!='-')
5725 strncpy(pli_filename,argv[z],253);
5729 /* Perform some error checking on the arguments
5730 and switches parsed from the command-line.
5731 If an error is encountered, print a message
5732 and exit gracefully. */
5736 fprintf(stderr,"\n%c*** ERROR: No transmitter site(s) specified!\n\n",7);
5740 for (x=0, y=0; x<txsites; x++)
5742 if (tx_site[x].lat==91.0 && tx_site[x].lon==361.0)
5744 fprintf(stderr,"\n*** ERROR: Transmitter site #%d not found!",x+1);
5751 fprintf(stderr,"%c\n\n",7);
5755 if ((coverage+LRmap+pli_filename[0])==0 && rx_site.lat==91.0 && rx_site.lon==361.0)
5757 if (max_range!=0.0 && txsites!=0)
5759 /* Plot topographic map of radius "max_range" */
5768 fprintf(stderr,"\n%c*** ERROR: No receiver site found or specified!\n\n",7);
5773 /* No major errors were detected. Whew! :-) */
5775 /* Adjust input parameters if -metric option is used */
5779 altitudeLR/=METERS_PER_FOOT; /* meters --> feet */
5780 max_range/=KM_PER_MILE; /* kilometers --> miles */
5781 altitude/=METERS_PER_FOOT; /* meters --> feet */
5784 /* If no SDF path was specified on the command line (-d), check
5785 for a path specified in the $HOME/.splat_path file. If the
5786 file is not found, then sdf_path[] remains NULL, and the
5787 current working directory is assumed to contain the SDF
5793 sprintf(string,"%s/.splat_path",env);
5794 fd=fopen(string,"r");
5798 fgets(string,253,fd);
5800 /* Remove <CR> and/or <LF> from string */
5802 for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0 && x<253; x++);
5805 strncpy(sdf_path,string,253);
5811 /* Ensure a trailing '/' is present in sdf_path */
5817 if (sdf_path[x-1]!='/' && x!=0)
5824 fprintf(stdout,"%s",header);
5827 if (pli_filename[0])
5829 y=LoadPLI(pli_filename);
5831 for (x=0; x<txsites; x++)
5832 PlaceMarker(tx_site[x]);
5835 PlaceMarker(rx_site);
5839 for (x=0; x<bfs; x++)
5840 LoadBoundaries(boundary_file[x]);
5845 for (x=0; x<cities; x++)
5846 LoadCities(city_file[x]);
5849 WritePPMLR(mapfile,geo);
5860 min_lon=(int)floor(tx_site[0].lon);
5861 max_lon=(int)floor(tx_site[0].lon);
5863 for (y=0, z=0; z<txsites; z++)
5865 txlat=(int)floor(tx_site[z].lat);
5866 txlon=(int)floor(tx_site[z].lon);
5874 if (LonDiff(txlon,min_lon)<0.0)
5877 if (LonDiff(txlon,max_lon)>0.0)
5883 rxlat=(int)floor(rx_site.lat);
5884 rxlon=(int)floor(rx_site.lon);
5892 if (LonDiff(rxlon,min_lon)<0.0)
5895 if (LonDiff(rxlon,max_lon)>0.0)
5899 /* Load the required SDF files */
5901 LoadTopoData(max_lon, min_lon, max_lat, min_lat);
5903 if (coverage | LRmap | topomap)
5908 for (z=0; z<txsites; z++)
5910 /* "Ball park" estimates used to load any additional
5911 SDF files required to conduct this analysis. */
5913 tx_range=sqrt(1.5*(tx_site[z].alt+GetElevation(tx_site[z])));
5916 rx_range=sqrt(1.5*altitudeLR);
5918 rx_range=sqrt(1.5*altitude);
5920 /* deg_range determines the maximum
5921 amount of topo data we read */
5923 deg_range=(tx_range+rx_range)/69.0;
5925 /* max_range sets the maximum size of the
5926 analysis. A small, non-zero amount can
5927 be used to shrink the size of the analysis
5928 and limit the amount of topo data read by
5929 SPLAT! A very large number will only increase
5930 the width of the analysis, not the size of
5934 max_range=tx_range+rx_range;
5936 if (max_range<(tx_range+rx_range))
5937 deg_range=max_range/69.0;
5939 /* Prevent the demand for a really wide coverage
5940 from allocating more slots than are available
5945 case 2: deg_limit=0.25;
5948 case 4: deg_limit=0.5;
5951 case 9: deg_limit=1.0;
5954 case 16: deg_limit=2.0;
5957 case 25: deg_limit=3.0;
5960 if (tx_site[z].lat<70.0)
5961 deg_range_lon=deg_range/cos(deg2rad*tx_site[z].lat);
5963 deg_range_lon=deg_range/cos(deg2rad*70.0);
5965 /* Correct for squares in degrees not being square in miles */
5967 if (deg_range>deg_limit)
5968 deg_range=deg_limit;
5970 if (deg_range_lon>deg_limit)
5971 deg_range_lon=deg_limit;
5974 north_min=(int)floor(tx_site[z].lat-deg_range);
5975 north_max=(int)floor(tx_site[z].lat+deg_range);
5977 west_min=(int)floor(tx_site[z].lon-deg_range_lon);
5982 while (west_min>=360)
5985 west_max=(int)floor(tx_site[z].lon+deg_range_lon);
5990 while (west_max>=360)
5993 if (north_min<min_lat)
5996 if (north_max>max_lat)
5999 if (LonDiff(west_min,min_lon)<0.0)
6002 if (LonDiff(west_max,max_lon)>0.0)
6006 /* Load any additional SDF files, if required */
6008 LoadTopoData(max_lon, min_lon, max_lat, min_lat);
6014 if (mapfile[0] && topomap==0)
6017 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);