Imported Upstream version 1.1.1

[debian/splat] / splat.cpp

/****************************************************************************
*      SPLAT: An RF Signal Propagation Loss and Terrain Analysis Tool       *
*                         Last update: 31-Mar-2006                          *
*****************************************************************************
*            Project started in 1997 by John A. Magliacane, KD2BD           *
*****************************************************************************
*                                                                           *
*     Extensively modified by J. D. McDonald in Jan. 2004 to include        *
*    the Longley-Rice propagation model using C++ code from NTIA/ITS.       *
*                                                                           *
*              See: http://flattop.its.bldrdoc.gov/itm.html                 *
*                                                                           *
*****************************************************************************
*                                                                           *
* This program is free software; you can redistribute it and/or modify it   *
* under the terms of the GNU General Public License as published by the     *
* Free Software Foundation; either version 2 of the License or any later    *
* version.                                                                  *
*                                                                           *
* This program is distributed in the hope that it will useful, but WITHOUT  *
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or     *
* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License     *
* for more details.                                                         *
*                                                                           *
*****************************************************************************
 g++ -Wall -O3 -s -lm -lbz2 -fomit-frame-pointer itm.cpp splat.cpp -o splat 
*****************************************************************************/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <bzlib.h>
#include <unistd.h>
#include "fontdata.h"
#include "smallfont.h"

#define GAMMA 2.5
#define MAXSLOTS 9
#define BZBUFFER 65536

#if MAXSLOTS==4
#define ARRAYSIZE 4950
#endif

#if MAXSLOTS==9
#define ARRAYSIZE 10870
#endif

#if MAXSLOTS==16
#define ARRAYSIZE 19240
#endif

#if MAXSLOTS==25
#define ARRAYSIZE 30025
#endif

char    string[255], sdf_path[255], opened=0, *splat_version={"1.1.1"};

double  TWOPI=6.283185307179586, HALFPI=1.570796326794896,
        PI=3.141592653589793, deg2rad=1.74532925199e-02,
        EARTHRADIUS=20902230.97, METERS_PER_MILE=1609.344,
        METERS_PER_FOOT=0.3048, earthradius, max_range=0.0;

int     min_north=90, max_north=-90, min_west=360, max_west=-1,
        max_elevation=-32768, min_elevation=32768, bzerror, maxdB=230;

struct site { double lat;
              double lon;
              double alt;
              char name[50];
            } site;

struct  path { float lat[ARRAYSIZE];
               float lon[ARRAYSIZE];
               float elevation[ARRAYSIZE];
               float distance[ARRAYSIZE];
               int length;
             } path;

struct dem { int min_north;
             int max_north;
             int min_west;
             int max_west;
             int max_el;
             int min_el;
             short data[1200][1200];
             unsigned char mask[1200][1200];
           } dem[MAXSLOTS];

struct LR { double eps_dielect; 
            double sgm_conductivity; 
            double eno_ns_surfref;
            double frq_mhz; 
            double conf; 
            double rel;
            int radio_climate;  
            int pol;
          } LR;

double elev_l[ARRAYSIZE+10];

void point_to_point(double elev[], double tht_m, double rht_m,
          double eps_dielect, double sgm_conductivity, double eno_ns_surfref,
          double frq_mhz, int radio_climate, int pol, double conf,
          double rel, double &dbloss, char *strmode, int &errnum);

double arccos(double x, double y)
{
        /* This function implements the arc cosine function,
           returning a value between 0 and TWOPI. */

        double result=0.0;

        if (y>0.0)
                result=acos(x/y);

        if (y<0.0)
                result=PI+acos(x/y);

        return result;
}

int ReduceAngle(double angle)
{
        /* This function normalizes the argument to
           an integer angle between 0 and 180 degrees */

        double temp;

        temp=acos(cos(angle*deg2rad));

        return (int)rint(temp/deg2rad);
}

char *dec2dms(double decimal)
{
        /* Converts decimal degrees to degrees, minutes, seconds,
           (DMS) and returns the result as a character string. */
 
        int degrees, minutes, seconds;
        double a, b, c, d;

        if (decimal<0.0)
                decimal=-decimal;

        a=floor(decimal);
        b=60.0*(decimal-a);
        c=floor(b);
        d=60.0*(b-c);

        degrees=(int)a;
        minutes=(int)c;
        seconds=(int)d;

        if (seconds<0)
                seconds=0;

        if (seconds>59)
                seconds=59;

        string[0]=0;
        sprintf(string,"%d%c %d\' %d\"", degrees, 176, minutes, seconds);
        return (string);
}

int OrMask(double lat, double lon, int value)
{
        /* Lines, text, markings, and coverage areas are stored in a
           mask that is combined with topology data when topographic
           maps are generated by SPLAT!.  This function sets bits in
           the mask based on the latitude and longitude of the area
           pointed to. */

        int x, y, indx, minlat, minlon;
        char found;

        minlat=(int)floor(lat);
        minlon=(int)floor(lon);

        for (indx=0, found=0; indx<MAXSLOTS && found==0;)
                if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
                        found=1;
                else
                        indx++;

        if (found)
        {
                x=(int)(1199.0*(lat-floor(lat)));
                y=(int)(1199.0*(lon-floor(lon)));

                dem[indx].mask[x][y]|=value;

                return (dem[indx].mask[x][y]);
        }

        else
                return -1;
}

int GetMask(double lat, double lon)
{
        /* This function returns the mask bits based on the latitude
           and longitude given. */

        return (OrMask(lat,lon,0));
}

double GetElevation(struct site location)
{
        /* This function returns the elevation (in feet) of any location
           represented by the digital elevation model data in memory.
           Function returns -5000.0 for locations not found in memory. */

        char found;
        int x, y, indx, minlat, minlon;
        double elevation;

        elevation=-5000.0;

        minlat=(int)floor(location.lat);
        minlon=(int)floor(location.lon);

        x=(int)(1199.0*(location.lat-floor(location.lat)));
        y=(int)(1199.0*(location.lon-floor(location.lon)));

        for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
        {
                if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
                {

                        elevation=3.28084*dem[indx].data[x][y];
                        found=1;
                }
        }
        
        return elevation;
}

double Distance(struct site site1, struct site site2)
{
        /* This function returns the great circle distance
           in miles between any two site locations. */

        double lat1, lon1, lat2, lon2, distance;

        lat1=site1.lat*deg2rad;
        lon1=site1.lon*deg2rad;
        lat2=site2.lat*deg2rad;
        lon2=site2.lon*deg2rad;

        distance=3959.0*acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos((lon1)-(lon2)));

        return distance;
}

double Azimuth(struct site source, struct site destination)
{
        /* This function returns the azimuth (in degrees) to the
           destination as seen from the location of the source. */

        double dest_lat, dest_lon, src_lat, src_lon,
               beta, azimuth, diff, num, den, fraction;

        dest_lat=destination.lat*deg2rad;
        dest_lon=destination.lon*deg2rad;

        src_lat=source.lat*deg2rad;
        src_lon=source.lon*deg2rad;
                
        /* Calculate Surface Distance */

        beta=acos(sin(src_lat)*sin(dest_lat)+cos(src_lat)*cos(dest_lat)*cos(src_lon-dest_lon));

        /* Calculate Azimuth */

        num=sin(dest_lat)-(sin(src_lat)*cos(beta));
        den=cos(src_lat)*sin(beta);
        fraction=num/den;

        /* Trap potential problems in acos() due to rounding */

        if (fraction>=1.0)
                fraction=1.0;

        if (fraction<=-1.0)
                fraction=-1.0;

        /* Calculate azimuth */

        azimuth=acos(fraction);

        /* Reference it to True North */

        diff=dest_lon-src_lon;

        if (diff<=-PI)
                diff+=TWOPI;

        if (diff>=PI)
                diff-=TWOPI;

        if (diff>0.0)
                azimuth=TWOPI-azimuth;

        return (azimuth/deg2rad);               
}

double ElevationAngle(struct site local, struct site remote)
{
        /* This function returns the angle of elevation (in degrees)
           of the remote location as seen from the local site.
           A positive result represents an angle of elevation (uptilt),
           while a negative result represents an angle of depression
           (downtilt), as referenced to a normal to the center of
           the earth. */
           
        register double a, b, dx;

        a=GetElevation(remote)+remote.alt+earthradius;
        b=GetElevation(local)+local.alt+earthradius;

        dx=5280.0*Distance(local,remote);

        /* Apply the Law of Cosines */

        return ((180.0*(acos(((b*b)+(dx*dx)-(a*a))/(2.0*b*dx)))/PI)-90.0);
}

void ReadPath(struct site source, struct site destination)
{
        /* This function generates a sequence of latitude and
           longitude positions between a source location and
           a destination along a great circle path, and stores
           elevation and distance information for points along
           that path in the "path" structure for later use. */

        int c;
        double azimuth, distance, lat1, lon1, beta,
               den, num, lat2, lon2, total_distance;
        struct site tempsite;

        lat1=source.lat*deg2rad;
        lon1=source.lon*deg2rad;
        azimuth=Azimuth(source,destination)*deg2rad;
        total_distance=Distance(source,destination);

        for (distance=0, c=0; distance<=total_distance; distance+=0.04)
        {
                beta=distance/3959.0;
                lat2=asin(sin(lat1)*cos(beta)+cos(azimuth)*sin(beta)*cos(lat1));
                num=cos(beta)-(sin(lat1)*sin(lat2));
                den=cos(lat1)*cos(lat2);

                if (azimuth==0.0 && (beta>HALFPI-lat1))
                        lon2=lon1+PI;

                else if (azimuth==HALFPI && (beta>HALFPI+lat1))
                        lon2=lon1+PI;

                else if (fabs(num/den)>1.0)
                        lon2=lon1;

                else
                {
                        if ((PI-azimuth)>=0.0)
                                lon2=lon1-arccos(num,den);
                        else
                                lon2=lon1+arccos(num,den);
                }
        
                while (lon2<0.0)
                        lon2+=TWOPI;

                while (lon2>TWOPI)
                        lon2-=TWOPI;
 
                lat2=lat2/deg2rad;
                lon2=lon2/deg2rad;

                if (c<ARRAYSIZE)
                {
                        path.lat[c]=lat2;
                        path.lon[c]=lon2;
                        tempsite.lat=lat2;
                        tempsite.lon=lon2;
                        path.elevation[c]=GetElevation(tempsite);
                        path.distance[c]=distance;
                        c++;
                }
        }

        /* Make sure exact destination point is recorded at path.length-1 */

        if (c<ARRAYSIZE)
        {
                path.lat[c]=destination.lat;
                path.lon[c]=destination.lon;
                path.elevation[c]=GetElevation(destination);
                path.distance[c]=total_distance;
                c++;
        }

        if (c<ARRAYSIZE)
                path.length=c;
        else
                path.length=ARRAYSIZE-1;
}

double AverageTerrain(struct site source, double azimuthx, double start_distance, double end_distance)
{
        /* This function returns the average terrain calculated in
           the direction of "azimuth" (degrees) between "start_distance"
           and "end_distance" (miles) from the source location.  If
           the terrain is all water (non-critical error), -5000.0 is
           returned.  If not enough SDF data has been loaded into
           memory to complete the survey (critical error), then
           -9999.0 is returned. */
 
        int c, samples, endpoint;
        double beta, lat1, lon1, lat2, lon2, num, den, azimuth, terrain=0.0;
        struct site destination;

        lat1=source.lat*deg2rad;
        lon1=source.lon*deg2rad;

        /* Generate a path of elevations between the source
           location and the remote location provided. */

        beta=end_distance/3959.0;

        azimuth=deg2rad*azimuthx;

        lat2=asin(sin(lat1)*cos(beta)+cos(azimuth)*sin(beta)*cos(lat1));
        num=cos(beta)-(sin(lat1)*sin(lat2));
        den=cos(lat1)*cos(lat2);

        if (azimuth==0.0 && (beta>HALFPI-lat1))
                lon2=lon1+PI;

        else if (azimuth==HALFPI && (beta>HALFPI+lat1))
                lon2=lon1+PI;

        else if (fabs(num/den)>1.0)
                lon2=lon1;

        else
        {
                if ((PI-azimuth)>=0.0)
                        lon2=lon1-arccos(num,den);
                else
                        lon2=lon1+arccos(num,den);
        }
        
        while (lon2<0.0)
                lon2+=TWOPI;

        while (lon2>TWOPI)
                lon2-=TWOPI;
 
        lat2=lat2/deg2rad;
        lon2=lon2/deg2rad;

        destination.lat=lat2;
        destination.lon=lon2;

        /* If SDF data is missing for the endpoint of
           the radial, then the average terrain cannot
           be accurately calculated.  Return -9999.0 */

        if (GetElevation(destination)<-4999.0)
                return (-9999.0);
        else
        {
                ReadPath(source,destination);

                endpoint=path.length;

                /* Shrink the length of the radial if the
                   outermost portion is not over U.S. land. */

                for (c=endpoint-1; c>=0 && path.elevation[c]==0.0; c--);

                endpoint=c+1;

                for (c=0, samples=0; c<endpoint; c++)
                {
                        if (path.distance[c]>=start_distance)
                        {
                                terrain+=path.elevation[c];
                                samples++;
                        }
                }

                if (samples==0)
                        terrain=-5000.0;  /* No land */
                else
                        terrain=(terrain/(double)samples);

                return terrain;
        }
}

double haat(struct site antenna)
{
        /* This function returns the antenna's Height Above Average
           Terrain (HAAT) based on FCC Part 73.313(d).  If a critical
           error occurs, such as a lack of SDF data to complete the
           survey, -5000.0 is returned. */

        int azi, c;
        char error=0;
        double terrain, avg_terrain, haat, sum=0.0;

        /* Calculate the average terrain between 2 and 10 miles
           from the antenna site at azimuths of 0, 45, 90, 135,
           180, 225, 270, and 315 degrees. */

        for (c=0, azi=0; azi<=315 && error==0; azi+=45)
        {
                terrain=AverageTerrain(antenna, (double)azi, 2.0, 10.0);

                if (terrain<-9998.0)  /* SDF data is missing */
                        error=1;

                if (terrain>-4999.0)  /* It's land, not water */
                {
                        sum+=terrain;  /* Sum of averages */
                        c++;
                }
        }

        if (error)
                return -5000.0;
        else
        {
                avg_terrain=(sum/(double)c);
                haat=(antenna.alt+GetElevation(antenna))-avg_terrain;
                return haat;
        }
}

float LonDiff(float lon1, float lon2)
{
        /* This function returns the short path longitudinal
           difference between longitude1 and longitude2 
           as an angle between -180.0 and +180.0 degrees.
           If lon1 is west of lon2, the result is positive.
           If lon1 is east of lon2, the result is negative. */

        float diff;

        diff=lon1-lon2;

        if (diff<=-180.0)
                diff+=360.0;

        if (diff>=180.0)
                diff-=360.0;

        return diff;
}

void PlaceMarker(struct site location)
{
        /* This function places text and marker data in the mask array
           for illustration on topographic maps generated by SPLAT!.
           By default, SPLAT! centers text information BELOW the marker,
           but may move it above, to the left, or to the right of the
           marker depending on how much room is available on the map,
           or depending on whether the area is already occupied by
           another marker or label.  If no room or clear space is
           available on the map to place the marker and its associated
           text, then the marker and text are not written to the map. */
 
        int     a, b, c, byte;
        char    ok2print, occupied;
        double  x, y, lat, lon, textx=0.0, texty=0.0, xmin, xmax,
                ymin, ymax, p1, p3, p6, p8, p12, p16, p24, label_length;

        xmin=min_north;
        xmax=max_north;
        ymin=min_west;
        ymax=max_west;
        lat=location.lat;
        lon=location.lon;


        if (lat<xmax && lat>xmin && (LonDiff(lon,ymax)<0.0) && (LonDiff(lon,ymin)>0.0))
        {
                p1=1.0/1200.0;
                p3=3.0/1200.0;
                p6=6.0/1200.0;
                p8=8.0/1200.0;
                p12=12.0/1200.0;
                p16=16.0/1200.0;
                p24=24.0/1200.0;
                ok2print=0;
                occupied=0;

                /* Is Marker Position Clear Of Text Or Other Markers? */

                for (x=lat-p3; (x<=xmax && x>=xmin && x<=lat+p3); x+=p1)
                        for (y=lon-p3; (LonDiff(y,ymax)<=0.0) && (LonDiff(y,ymin)>=0.0) && (LonDiff(y,lon+p3)<=0.0); y+=p1)
                                occupied|=(GetMask(x,y)&2);

                if (occupied==0)
                {
                        /* Determine Where Text Can Be Positioned */

                        /* label_length=length in pixels.
                           Each character is 8 pixels wide. */

                        label_length=p1*(double)(strlen(location.name)<<3);

                        if ((LonDiff(lon+label_length,ymax)<=0.0) && (LonDiff(lon-label_length,ymin)>=0.0))
                        {
                                /* Default: Centered Text */

                                texty=lon+label_length/2.0;

                                if ((lat-p8)>=p16)
                                {
                                        /* Position Text Below The Marker */

                                        textx=lat-p8;

                                        x=textx;
                                        y=texty;
        
                                        /* Is This Position Clear Of
                                           Text Or Other Markers? */

                                        for (a=0, occupied=0; a<16; a++)
                                        {
                                                for (b=0; b<(int)strlen(location.name); b++)
                                                        for (c=0; c<8; c++, y-=p1)
                                                                occupied|=(GetMask(x,y)&2);
                                                x-=p1;
                                                y=texty;
                                        }

                                        x=textx;
                                        y=texty;

                                        if (occupied==0)
                                                ok2print=1;
                                }

                                else
                                {
                                        /* Position Text Above The Marker */

                                        textx=lat+p24;

                                        x=textx;
                                        y=texty;
        
                                        /* Is This Position Clear Of
                                           Text Or Other Markers? */

                                        for (a=0, occupied=0; a<16; a++)
                                        {
                                                for (b=0; b<(int)strlen(location.name); b++)
                                                        for (c=0; c<8; c++, y-=p1)
                                                                occupied|=(GetMask(x,y)&2);
                                                x-=p1;
                                                y=texty;
                                        }

                                        x=textx;
                                        y=texty;

                                        if (occupied==0)
                                                ok2print=1;
                                }
                        }

                        if (ok2print==0)
                        {
                                if (LonDiff(lon-label_length,ymin)>=0.0)
                                {
                                        /* Position Text To The
                                           Right Of The Marker */

                                        textx=lat+p6;
                                        texty=lon-p12;

                                        x=textx;
                                        y=texty;
        
                                        /* Is This Position Clear Of
                                           Text Or Other Markers? */

                                        for (a=0, occupied=0; a<16; a++)
                                        {
                                                for (b=0; b<(int)strlen(location.name); b++)
                                                        for (c=0; c<8; c++, y-=p1)
                                                                occupied|=(GetMask(x,y)&2);
                                                x-=p1;
                                                y=texty;
                                        }

                                        x=textx;
                                        y=texty;

                                        if (occupied==0)
                                                ok2print=1;
                                }

                                else
                                {
                                        /* Position Text To The
                                           Left Of The Marker */

                                        textx=lat+p6;
                                        texty=lon+p8+(label_length);

                                        x=textx;
                                        y=texty;
        
                                        /* Is This Position Clear Of
                                           Text Or Other Markers? */

                                        for (a=0, occupied=0; a<16; a++)
                                        {
                                                for (b=0; b<(int)strlen(location.name); b++)
                                                        for (c=0; c<8; c++, y-=p1)
                                                                occupied|=(GetMask(x,y)&2);
                                                x-=p1;
                                                y=texty;
                                        }

                                        x=textx;
                                        y=texty;

                                        if (occupied==0)
                                                ok2print=1;
                                }
                        }

                        /* textx and texty contain the latitude and longitude
                           coordinates that describe the placement of the text
                           on the map. */
        
                        if (ok2print)
                        {
                                /* Draw Text */

                                x=textx;
                                y=texty;

                                for (a=0; a<16 && ok2print; a++)
                                {
                                        for (b=0; b<(int)strlen(location.name); b++)
                                        {
                                                byte=fontdata[16*(location.name[b])+a];

                                                for (c=128; c>0; c=c>>1, y-=p1)
                                                        if (byte&c)
                                                                OrMask(x,y,2);
                                        }

                                        x-=p1;
                                        y=texty;
                                }
        
                                /* Draw Square Marker Centered
                                   On Location Specified */
        
                                for (x=lat-p3; (x<=xmax && x>=xmin && x<=lat+p3); x+=p1)
                                        for (y=lon-p3; (LonDiff(y,ymax)<=0.0) && (LonDiff(y,ymin)>=0.0) && (LonDiff(y,lon+p3)<=0.0); y+=p1)
                                                OrMask(x,y,2);
                        }
                }
        }
}

double ReadBearing(char *input)
{
        /* This function takes numeric input in the form of a character
           string, and returns an equivalent bearing in degrees as a
           decimal number (double).  The input may either be expressed
           in decimal format (40.139722) or degree, minute, second
           format (40 08 23).  This function also safely handles
           extra spaces found either leading, trailing, or
           embedded within the numbers expressed in the
           input string.  Decimal seconds are permitted. */
 
        double seconds, bearing=0.0;
        char string[20];
        int a, b, length, degrees, minutes;

        /* Copy "input" to "string", and ignore any extra
           spaces that might be present in the process. */

        string[0]=0;
        length=strlen(input);

        for (a=0, b=0; a<length && a<18; a++)
        {
                if ((input[a]!=32 && input[a]!='\n') || (input[a]==32 && input[a+1]!=32 && b!=0))
                {
                        string[b]=input[a];
                        b++;
                }        
        }

        string[b]=0;

        /* Count number of spaces in the clean string. */

        length=strlen(string);

        for (a=0, b=0; a<length; a++)
                if (string[a]==32)
                        b++;

        if (b==0)  /* Decimal Format (40.139722) */
                sscanf(string,"%lf",&bearing);

        if (b==2)  /* Degree, Minute, Second Format (40 08 23) */
        {
                sscanf(string,"%d %d %lf",&degrees, &minutes, &seconds);
                bearing=(double)degrees+((double)minutes/60)+(seconds/3600);
        }

        /* Anything else returns a 0.0 */

        if (bearing>360.0 || bearing<-90.0)
                bearing=0.0;

        return bearing;
}

struct site LoadQTH(char *filename)
{
        /* This function reads SPLAT! .qth (site location) files.
           The latitude and longitude may be expressed either in
           decimal degrees, or in degree, minute, second format.
           Antenna height is assumed to be expressed in feet above
           ground level (AGL), unless followed by the letter 'M',
           or 'm', or by the word "meters" or "Meters", in which
           case meters is assumed, and is handled accordingly. */

        int x;
        char string[50], qthfile[255];
        struct site tempsite;
        FILE *fd=NULL;

        for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
                qthfile[x]=filename[x];

        qthfile[x]='.';
        qthfile[x+1]='q';
        qthfile[x+2]='t';
        qthfile[x+3]='h';
        qthfile[x+4]=0;

        tempsite.lat=91.0;
        tempsite.lon=361.0;
        tempsite.alt=0.0;
        tempsite.name[0]=0;

        fd=fopen(qthfile,"r");

        if (fd!=NULL)
        {
                /* Site Name */
                fgets(string,49,fd);

                /* Strip <CR> and/or <LF> from end of site name */

                for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0; tempsite.name[x]=string[x], x++);

                tempsite.name[x]=0;

                /* Site Latitude */
                fgets(string,49,fd);
                tempsite.lat=ReadBearing(string);

                /* Site Longitude */
                fgets(string,49,fd);
                tempsite.lon=ReadBearing(string);

                /* Antenna Height */
                fgets(string,49,fd);
                fclose(fd);

                /* Remove <CR> and/or <LF> from antenna height string */
                for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0; x++);

                string[x]=0;

                /* Antenna height may either be in feet or meters.
                   If the letter 'M' or 'm' is discovered in
                   the string, then this is an indication that
                   the value given is expressed in meters, and
                   must be converted to feet before exiting. */

                for (x=0; string[x]!='M' && string[x]!='m' && string[x]!=0 && x<48; x++);
                if (string[x]=='M' || string[x]=='m')
                {
                        string[x]=0;
                        sscanf(string,"%lf",&tempsite.alt);
                        tempsite.alt*=3.28084;
                }

                else
                {
                        string[x]=0;
                        sscanf(string,"%lf",&tempsite.alt);
                }
        }

        return tempsite;
}

int LoadSDF_SDF(char *name)
{
        /* This function reads uncompressed SPLAT Data Files (.sdf)
           containing digital elevation model data into memory.
           Elevation data, maximum and minimum elevations, and
           quadrangle limits are stored in the first available
           dem[] structure. */

        int x, y, data, indx, minlat, minlon, maxlat, maxlon;
        char found, free_slot=0, line[20], sdf_file[255], path_plus_name[255];
        FILE *fd;

        for (x=0; name[x]!='.' && name[x]!=0 && x<250; x++)
                sdf_file[x]=name[x];

        sdf_file[x]=0;

        /* Parse filename for minimum latitude and longitude values */

        sscanf(sdf_file,"%d:%d:%d:%d",&minlat,&maxlat,&minlon,&maxlon);

        sdf_file[x]='.';
        sdf_file[x+1]='s';
        sdf_file[x+2]='d';
        sdf_file[x+3]='f';
        sdf_file[x+4]=0;

        /* Is it already in memory? */

        for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
        {
                if (minlat==dem[indx].min_north && minlon==dem[indx].min_west && maxlat==dem[indx].max_north && maxlon==dem[indx].max_west)
                        found=1;
        }

        /* Is room available to load it? */

        if (found==0)
        {       
                for (indx=0, free_slot=0; indx<MAXSLOTS && free_slot==0; indx++)
                        if (dem[indx].max_north==-90)
                                free_slot=1;
        }

        indx--;

        if (free_slot && found==0 && indx>=0 && indx<MAXSLOTS)
        {
                /* Search for SDF file in current working directory first */

                strncpy(path_plus_name,sdf_file,255);

                fd=fopen(path_plus_name,"rb");

                if (fd==NULL)
                {
                        /* Next, try loading SDF file from path specified
                           in $HOME/.splat_path file or by -d argument */

                        strncpy(path_plus_name,sdf_path,255);
                        strncat(path_plus_name,sdf_file,255);

                        fd=fopen(path_plus_name,"rb");
                }

                if (fd!=NULL)
                {
                        fprintf(stdout,"Loading \"%s\" into slot %d...",path_plus_name,indx+1);
                        fflush(stdout);

                        fgets(line,19,fd);
                        sscanf(line,"%d",&dem[indx].max_west);

                        fgets(line,19,fd);
                        sscanf(line,"%d",&dem[indx].min_north);

                        fgets(line,19,fd);
                        sscanf(line,"%d",&dem[indx].min_west);

                        fgets(line,19,fd);
                        sscanf(line,"%d",&dem[indx].max_north);

                        for (x=0; x<1200; x++)
                                for (y=0; y<1200; y++)
                                {
                                        fgets(line,19,fd);
                                        sscanf(line,"%d",&data);

                                        dem[indx].data[x][y]=data;

                                        if (data>dem[indx].max_el)
                                                dem[indx].max_el=data;

                                        if (data<dem[indx].min_el)
                                                dem[indx].min_el=data;
                                }

                        fclose(fd);

                        if (dem[indx].min_el<min_elevation)
                                min_elevation=dem[indx].min_el;

                        if (dem[indx].max_el>max_elevation)
                                max_elevation=dem[indx].max_el;

                        if (max_north==-90)
                                max_north=dem[indx].max_north;

                        else if (dem[indx].max_north>max_north)
                                max_north=dem[indx].max_north;

                        if (min_north==90)
                                min_north=dem[indx].min_north;

                        else if (dem[indx].min_north<min_north)
                                min_north=dem[indx].min_north;

                        if (max_west==-1)
                                max_west=dem[indx].max_west;

                        else
                        {
                                if (abs(dem[indx].max_west-max_west)<180)
                                {
                                        if (dem[indx].max_west>max_west)
                                                max_west=dem[indx].max_west;
                                }

                                else
                                {
                                        if (dem[indx].max_west<max_west)
                                                max_west=dem[indx].max_west;
                                }
                        }

                        if (min_west==360)
                                min_west=dem[indx].min_west;

                        else
                        {
                                if (abs(dem[indx].min_west-min_west)<180)
                                {
                                        if (dem[indx].min_west<min_west)
                                                min_west=dem[indx].min_west;
                                }

                                else
                                {
                                        if (dem[indx].min_west>min_west)
                                                min_west=dem[indx].min_west;
                                }
                        }

                        fprintf(stdout," Done!\n");
                        fflush(stdout);
                        return 1;
                }

                else
                        return -1;
        }

        else
                return 0;
}

char *BZfgets(BZFILE *bzfd, unsigned length)
{
        /* This function returns at most one less than 'length' number
           of characters from a bz2 compressed file whose file descriptor
           is pointed to by *bzfd.  In operation, a buffer is filled with
           uncompressed data (size = BZBUFFER), which is then parsed
           and doled out as NULL terminated character strings every time
           this function is invoked.  A NULL string indicates an EOF
           or error condition. */

        static int x, y, nBuf;
        static char buffer[BZBUFFER+1], output[BZBUFFER+1];
        char done=0;

        if (opened!=1 && bzerror==BZ_OK)
        {
                /* First time through.  Initialize everything! */

                x=0;
                y=0;
                nBuf=0;
                opened=1;
                output[0]=0;
        }

        do
        {
                if (x==nBuf && bzerror!=BZ_STREAM_END && bzerror==BZ_OK && opened)
                {
                        /* Uncompress data into a static buffer */

                        nBuf=BZ2_bzRead(&bzerror, bzfd, buffer, BZBUFFER);
                        buffer[nBuf]=0;
                        x=0;
                }

                /* Build a string from buffer contents */

                output[y]=buffer[x];

                if (output[y]=='\n' || output[y]==0 || y==(int)length-1)
                {
                        output[y+1]=0;
                        done=1;
                        y=0;
                }

                else
                        y++;
                x++;

        } while (done==0);

        if (output[0]==0)
                opened=0;

        return (output);
}

int LoadSDF_BZ(char *name)
{
        /* This function reads .bz2 compressed SPLAT Data Files containing
           digital elevation model data into memory.  Elevation data,
           maximum and minimum elevations, and quadrangle limits are
           stored in the first available dem[] structure. */

        int x, y, data, indx, minlat, minlon, maxlat, maxlon;
        char found, free_slot=0, sdf_file[255], path_plus_name[255];
        FILE *fd;
        BZFILE *bzfd;

        for (x=0; name[x]!='.' && name[x]!=0 && x<247; x++)
                sdf_file[x]=name[x];

        sdf_file[x]=0;

        /* Parse sdf_file name for minimum latitude and longitude values */

        sscanf(sdf_file,"%d:%d:%d:%d",&minlat,&maxlat,&minlon,&maxlon);

        sdf_file[x]='.';
        sdf_file[x+1]='s';
        sdf_file[x+2]='d';
        sdf_file[x+3]='f';
        sdf_file[x+4]='.';
        sdf_file[x+5]='b';
        sdf_file[x+6]='z';
        sdf_file[x+7]='2';
        sdf_file[x+8]=0;

        /* Is it already in memory? */

        for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
        {
                if (minlat==dem[indx].min_north && minlon==dem[indx].min_west && maxlat==dem[indx].max_north && maxlon==dem[indx].max_west)
                        found=1;
        }

        /* Is room available to load it? */

        if (found==0)
        {       
                for (indx=0, free_slot=0; indx<MAXSLOTS && free_slot==0; indx++)
                        if (dem[indx].max_north==-90)
                                free_slot=1;
        }

        indx--;

        if (free_slot && found==0 && indx>=0 && indx<MAXSLOTS)
        {
                /* Search for SDF file in current working directory first */

                strncpy(path_plus_name,sdf_file,255);

                fd=fopen(path_plus_name,"rb");
                bzfd=BZ2_bzReadOpen(&bzerror,fd,0,0,NULL,0);

                if (fd==NULL || bzerror!=BZ_OK)
                {
                        /* Next, try loading SDF file from path specified
                           in $HOME/.splat_path file or by -d argument */

                        strncpy(path_plus_name,sdf_path,255);
                        strncat(path_plus_name,sdf_file,255);

                        fd=fopen(path_plus_name,"rb");
                        bzfd=BZ2_bzReadOpen(&bzerror,fd,0,0,NULL,0);
                }

                if (fd!=NULL && bzerror==BZ_OK)
                {
                        fprintf(stdout,"Loading \"%s\" into slot %d...",path_plus_name,indx+1);
                        fflush(stdout);

                        sscanf(BZfgets(bzfd,255),"%d",&dem[indx].max_west);
                        sscanf(BZfgets(bzfd,255),"%d",&dem[indx].min_north);
                        sscanf(BZfgets(bzfd,255),"%d",&dem[indx].min_west);
                        sscanf(BZfgets(bzfd,255),"%d",&dem[indx].max_north);
        
                        for (x=0; x<1200; x++)
                                for (y=0; y<1200; y++)
                                {
                                        sscanf(BZfgets(bzfd,20),"%d",&data);

                                        dem[indx].data[x][y]=data;

                                        if (data>dem[indx].max_el)
                                                dem[indx].max_el=data;

                                        if (data<dem[indx].min_el)
                                                dem[indx].min_el=data;
                                }

                        fclose(fd);

                        BZ2_bzReadClose(&bzerror,bzfd);

                        if (dem[indx].min_el<min_elevation)
                                min_elevation=dem[indx].min_el;
        
                        if (dem[indx].max_el>max_elevation)
                                max_elevation=dem[indx].max_el;

                        if (max_north==-90)
                                max_north=dem[indx].max_north;

                        else if (dem[indx].max_north>max_north)
                                max_north=dem[indx].max_north;

                        if (min_north==90)
                                min_north=dem[indx].min_north;

                        else if (dem[indx].min_north<min_north)
                                min_north=dem[indx].min_north;

                        if (max_west==-1)
                                max_west=dem[indx].max_west;

                        else
                        {
                                if (abs(dem[indx].max_west-max_west)<180)
                                {
                                        if (dem[indx].max_west>max_west)
                                                max_west=dem[indx].max_west;
                                }

                                else
                                {
                                        if (dem[indx].max_west<max_west)
                                                max_west=dem[indx].max_west;
                                }
                        }

                        if (min_west==360)
                                min_west=dem[indx].min_west;

                        else
                        {
                                if (abs(dem[indx].min_west-min_west)<180)
                                {
                                        if (dem[indx].min_west<min_west)
                                                min_west=dem[indx].min_west;
                                }

                                else
                                {
                                        if (dem[indx].min_west>min_west)
                                                min_west=dem[indx].min_west;
                                }
                        }

                        fprintf(stdout," Done!\n");
                        fflush(stdout);
                        return 1;
                }
                else
                        return -1;
        }
        else
                return 0;
}

char LoadSDF(char *name)
{
        /* This function loads the requested SDF file from the filesystem.
           It first tries to invoke the LoadSDF_SDF() function to load an
           uncompressed SDF file (since uncompressed files load slightly
           faster).  If that attempt fails, then it tries to load a
           compressed SDF file by invoking the LoadSDF_BZ() function.
           If that fails, then we can assume that no elevation data
           exists for the region requested, and that the region
           requested must be entirely over water. */

        int x, y, indx, minlat, minlon, maxlat, maxlon;
        char found, free_slot=0;
        int  return_value=-1;

        /* Try to load an uncompressed SDF first. */

        return_value=LoadSDF_SDF(name);

        /* If that fails, try loading a compressed SDF. */

        if (return_value==0 || return_value==-1)
                return_value=LoadSDF_BZ(name);

        /* If neither format can be found, then assume the area is water. */

        if (return_value==0 || return_value==-1)
        {
                /* Parse SDF name for minimum latitude and longitude values */

                sscanf(name,"%d:%d:%d:%d",&minlat,&maxlat,&minlon,&maxlon);

                /* Is it already in memory? */

                for (indx=0, found=0; indx<MAXSLOTS && found==0; indx++)
                {
                        if (minlat==dem[indx].min_north && minlon==dem[indx].min_west && maxlat==dem[indx].max_north && maxlon==dem[indx].max_west)
                                found=1;
                }

                /* Is room available to load it? */

                if (found==0)
                {       
                        for (indx=0, free_slot=0; indx<MAXSLOTS && free_slot==0; indx++)
                                if (dem[indx].max_north==-90)
                                        free_slot=1;
                }

                indx--;

                if (free_slot && found==0 && indx>=0 && indx<MAXSLOTS)
                {
                        fprintf(stdout,"Region  \"%s\" assumed as sea-level into slot %d...",name,indx+1);
                        fflush(stdout);

                        dem[indx].max_west=maxlon;
                        dem[indx].min_north=minlat;
                        dem[indx].min_west=minlon;
                        dem[indx].max_north=maxlat;

                        /* Fill DEM with sea-level topography */

                        for (x=0; x<1200; x++)
                                for (y=0; y<1200; y++)
                                {
                                        dem[indx].data[x][y]=0;

                                        if (dem[indx].min_el>0)
                                                dem[indx].min_el=0;
                                }

                        if (dem[indx].min_el<min_elevation)
                                min_elevation=dem[indx].min_el;

                        if (dem[indx].max_el>max_elevation)
                                max_elevation=dem[indx].max_el;

                        if (max_north==-90)
                                max_north=dem[indx].max_north;

                        else if (dem[indx].max_north>max_north)
                                max_north=dem[indx].max_north;

                        if (min_north==90)
                                min_north=dem[indx].min_north;

                        else if (dem[indx].min_north<min_north)
                                min_north=dem[indx].min_north;

                        if (max_west==-1)
                                max_west=dem[indx].max_west;

                        else
                        {
                                if (abs(dem[indx].max_west-max_west)<180)
                                {
                                        if (dem[indx].max_west>max_west)
                                                max_west=dem[indx].max_west;
                                }

                                else
                                {
                                        if (dem[indx].max_west<max_west)
                                                max_west=dem[indx].max_west;
                                }
                        }

                        if (min_west==360)
                                min_west=dem[indx].min_west;

                        else
                        {
                                if (abs(dem[indx].min_west-min_west)<180)
                                {
                                        if (dem[indx].min_west<min_west)
                                                min_west=dem[indx].min_west;
                                }

                                else
                                {
                                        if (dem[indx].min_west>min_west)
                                                min_west=dem[indx].min_west;
                                }
                        }

                        fprintf(stdout," Done!\n");
                        fflush(stdout);

                        return_value=1;
                }
        }

        return return_value;
}

void LoadCities(char *filename)
{
        /* This function reads SPLAT! city/site files, and plots
           the locations and names of the cities and site locations
           read on topographic maps generated by SPLAT! */

        int x, y, z;
        char input[80], str[3][80];
        struct site city_site;
        FILE *fd=NULL;

        fd=fopen(filename,"r");

        if (fd!=NULL)
        {
                fgets(input,78,fd);

                fprintf(stdout,"Reading \"%s\"... ",filename);
                fflush(stdout);

                while (fd!=NULL && feof(fd)==0)
                {
                        /* Parse line for name, latitude, and longitude */

                        for (x=0, y=0, z=0; x<78 && input[x]!=0 && z<3; x++)
                        {
                                if (input[x]!=',' && y<78)
                                {
                                        str[z][y]=input[x];
                                        y++;
                                }

                                else
                                {
                                        str[z][y]=0;
                                        z++;
                                        y=0;
                                }
                        }

                        strncpy(city_site.name,str[0],49);
                        city_site.lat=ReadBearing(str[1]);
                        city_site.lon=ReadBearing(str[2]);
                        city_site.alt=0.0;

                        PlaceMarker(city_site);

                        fgets(input,78,fd);
                }

                fclose(fd);
                fprintf(stdout,"Done!\n");
                fflush(stdout);
        }
        else
                fprintf(stderr,"*** ERROR: \"%s\": not found!\n",filename);
}

void LoadBoundaries(char *filename)
{
        /* This function reads Cartographic Boundary Files available from
           the U.S. Census Bureau, and plots the data contained in those
           files on the PPM Map generated by SPLAT!.  Such files contain
           the coordinates that describe the boundaries of cities,
           counties, and states. */

        int x;
        double lat0, lon0, lat1, lon1;
        char string[80];
        struct site source, destination;
        FILE *fd=NULL;

        fd=fopen(filename,"r");

        if (fd!=NULL)
        {
                fgets(string,78,fd);

                fprintf(stdout,"Reading \"%s\"... ",filename);
                fflush(stdout);

                do
                {
                        fgets(string,78,fd);
                        sscanf(string,"%lf %lf", &lon0, &lat0);
                        fgets(string,78,fd);

                        do
                        {
                                sscanf(string,"%lf %lf", &lon1, &lat1);

                                lon0=fabs(lon0);
                                lon1=fabs(lon1);

                                source.lat=lat0;
                                source.lon=lon0;
                                destination.lat=lat1;
                                destination.lon=lon1;

                                ReadPath(source,destination);

                                for (x=0; x<path.length; x++)
                                        OrMask(path.lat[x],path.lon[x],4);

                                lat0=lat1;
                                lon0=lon1;

                                fgets(string,78,fd);

                        } while (strncmp(string,"END",3)!=0 && feof(fd)==0);

                        fgets(string,78,fd);

                } while (strncmp(string,"END",3)!=0 && feof(fd)==0);

                fclose(fd);

                fprintf(stdout,"Done!\n");
                fflush(stdout);
        }

        else
                fprintf(stderr,"*** ERROR: \"%s\": not found!\n",filename);
}

void ReadLRParm(char *txsite_filename)
{
        /* This function reads Longley-Rice parameter data for the
           transmitter site.  The file name is the same as the txsite,
           except the filename extension is .lrp.  If the needed file
           is not found, then the file "splat.lrp" is read from the
           current working directory.  Failure to load this file will
           result in the default parameters hard coded into this
           being used and written to "splat.lrp". */

        double din;
        char filename[255], lookup[256], string[80];
        int iin, ok=0, x;
        FILE *fd=NULL, *outfile=NULL;

        /* Default parameters in case things go bad */

        LR.eps_dielect=15.0;
        LR.sgm_conductivity=0.005;
        LR.eno_ns_surfref=301.0;
        LR.frq_mhz=300.0;
        LR.radio_climate=5;
        LR.pol=0;
        LR.conf=0.50;
        LR.rel=0.50;

        /* Modify txsite filename to one with a .lrp extension. */

        strncpy(filename,txsite_filename,255);

        for (x=0; filename[x]!='.' && filename[x]!=0 && filename[x]!='\n' && x<249; x++);

        filename[x]='.';
        filename[x+1]='l';
        filename[x+2]='r';
        filename[x+3]='p';
        filename[x+4]=0;

        /* Small lookup table to parse file, pass
           numeric data, and ignore comments. */

        for (x=0; x<=255; x++)
                lookup[x]=0;

        /* Valid characters */

        for (x=48; x<=57; x++)
                lookup[x]=x;

        lookup[32]=32;
        lookup[46]='.';

        fd=fopen(filename,"r");

        if (fd==NULL)
        {
                /* Load default "splat.lrp" file */

                strncpy(filename,"splat.lrp\0",10);
                fd=fopen(filename,"r");
        }

        if (fd!=NULL)
        {
                fgets(string,80,fd);

                for (x=0; lookup[(int)string[x]] && x<20; x++)
                        string[x]=lookup[(int)string[x]];

                string[x]=0;

                ok=sscanf(string,"%lf", &din);

                if (ok)
                {
                        LR.eps_dielect=din;

                        fgets(string,80,fd);

                        for (x=0; lookup[(int)string[x]] && x<20; x++)
                                string[x]=lookup[(int)string[x]];

                        string[x]=0;

                        ok=sscanf(string,"%lf", &din);
                }

                if (ok)
                {
                        LR.sgm_conductivity=din;

                        fgets(string,80,fd);

                        for (x=0; lookup[(int)string[x]] && x<20; x++)
                                string[x]=lookup[(int)string[x]];

                        string[x]=0;

                        ok=sscanf(string,"%lf", &din);
                }

                if (ok)
                {
                        LR.eno_ns_surfref=din;

                        fgets(string,80,fd);

                        for (x=0; lookup[(int)string[x]] && x<20; x++)
                                string[x]=lookup[(int)string[x]];

                        string[x]=0;

                        ok=sscanf(string,"%lf", &din);
                }

                if (ok)
                {
                        LR.frq_mhz=din;

                        fgets(string,80,fd);

                        for (x=0; lookup[(int)string[x]] && x<20; x++)
                                string[x]=lookup[(int)string[x]];

                        string[x]=0;

                        ok=sscanf(string,"%d", &iin);
                }

                if (ok)
                {
                        LR.radio_climate=iin;

                        fgets(string,80,fd);

                        for (x=0; lookup[(int)string[x]] && x<20; x++)
                                string[x]=lookup[(int)string[x]];

                        string[x]=0;

                        ok=sscanf(string,"%d", &iin);
                }

                if (ok)
                {
                        LR.pol=iin;

                        fgets(string,80,fd);

                        for (x=0; lookup[(int)string[x]] && x<20; x++)
                                string[x]=lookup[(int)string[x]];

                        string[x]=0;

                        ok=sscanf(string,"%lf", &din);
                }

                if (ok)
                {
                        LR.conf=din;

                        fgets(string,80,fd);

                        for (x=0; lookup[(int)string[x]] && x<20; x++)
                                string[x]=lookup[(int)string[x]];

                        string[x]=0;

                        ok=sscanf(string,"%lf", &din);
                }

                if (ok)
                        LR.rel=din;

                fclose(fd);
        } 

        if (fd==NULL)
        {
                /* Create a "splat.lrp" file since one
                   could not be successfully loaded. */

                outfile=fopen("splat.lrp","w");

                fprintf(outfile,"%.3f\t; Earth Dielectric Constant (Relative permittivity)\n",LR.eps_dielect);

                fprintf(outfile,"%.3f\t; Earth Conductivity (Siemens per meter)\n", LR.sgm_conductivity);

                fprintf(outfile,"%.3f\t; Atmospheric Bending Constant (N-Units)\n",LR.eno_ns_surfref);

                fprintf(outfile,"%.3f\t; Frequency in MHz (20 MHz to 20 GHz)\n", LR.frq_mhz);

                fprintf(outfile,"%d\t; Radio Climate\n",LR.radio_climate);

                fprintf(outfile,"%d\t; Polarization (0 = Horizontal, 1 = Vertical)\n", LR.pol);

                fprintf(outfile,"%.2f\t; Fraction of situations\n",LR.conf);

                fprintf(outfile, "%.2f\t; Fraction of time\n",LR.rel);

                fprintf(outfile,"\nPlease consult SPLAT! documentation for the meaning and use of this data.\n");

                fclose(outfile);

                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);
        }

        if (fd==NULL || ok==0)
                fprintf(stderr,"Longley-Rice default parameters have been assumed for this analysis.\n");
}

struct site los(struct site source, struct site destination)
{
        /* This function determines whether a line-of-sight path
           unobstructed by terrain exists between source (transmitter)
           and destination (receiver) based on the geographical
           locations of the two sites, their respective antenna
           heights above ground, and the terrain between them.
           A site structure is returned upon completion.  If the
           first character of site.name is ' ', then a clear path
           exists between source and destination.  If the first
           character is '*', then an obstruction exists, and the
           site.lat and site.lon elements of the structure provide
           the geographical location of the obstruction. */
           
        int x;
        char block;
        struct site test, blockage;
        register double distance, tx_alt, rx_alt,
                 cos_xmtr_angle, cos_test_angle, test_alt;

        ReadPath(source,destination);

        distance=5280.0*Distance(source,destination);
        tx_alt=earthradius+source.alt+GetElevation(source);
        rx_alt=earthradius+destination.alt+GetElevation(destination);

        /* Elevation angle of the xmtr (source) as seen by the rcvr */

        cos_xmtr_angle=((rx_alt*rx_alt)+(distance*distance)-(tx_alt*tx_alt))/(2.0*rx_alt*distance);

        /* Determine the elevation angle of each discrete location
           along the path between the receiver and transmitter.

           Since obstructions are more likely due to terrain effects
           closest to the receiver rather than farther away, we start
           looking for potential obstructions from the receiver's
           location, and work our way towards the transmitter.
           This loop is broken when the first obstruction is
           detected.  If we can travel all the way to the transmitter
           without detecting an obstruction, then we have a clear
           unobstructed path between transmitter and receiver. */

        for (x=path.length-1, block=0; x>0 && block==0; x--)
        {
                /* Build a structure for each test
                   point along the path to be surveyed. */

                test.lat=path.lat[x];
                test.lon=path.lon[x];

                /* Measure the distance between the
                   test point and the receiver locations */

                distance=5280.0*Distance(test,destination);
                test_alt=earthradius+path.elevation[x];

                /* Determine the cosine of the elevation of the test
                   point as seen from the location of the receiver */

                cos_test_angle=((rx_alt*rx_alt)+(distance*distance)-(test_alt*test_alt))/(2.0*rx_alt*distance);

                /* If the elevation angle to the test point (as seen from
                   the receiver) is greater than the elevation angle to the
                   transmitter (as seen by the receiver), then we have a
                   path obstructed by terrain.  Note: Since we're comparing
                   the cosines of these angles rather than the angles
                   themselves (eliminating the call to acos() saves
                   considerable time), the following "if" statement is
                   reversed from what it would normally be if the angles
                   were compared. */

                if (cos_xmtr_angle>cos_test_angle)
                {
                        block=1;
                        blockage.lat=path.lat[x];
                        blockage.lon=path.lon[x];
                        blockage.alt=path.elevation[x];
                        blockage.name[0]='*';
                }
        }

        if (block==0)
        {
                blockage.lat=0.0;
                blockage.lon=0.0;
                blockage.alt=0.0;
                blockage.name[0]=' ';
        }

        return blockage;
}

void PlotPath(struct site source, struct site destination, char mask_value)
{
        /* This function analyzes the path between the source and
           destination locations.  It determines which points along
           the path have line-of-sight visibility to the source.
           Points along with path having line-of-sight visibility
           to the source at an AGL altitude equal to that of the
           destination location are stored by setting bit 1 in the
           mask[][] array, which are displayed in green when PPM
           maps are later generated by SPLAT!. */

        char block;
        int x, y;
        register double cos_xmtr_angle, cos_test_angle, test_alt;
        double distance, rx_alt, tx_alt;

        ReadPath(source,destination);

        for (y=0; y<path.length; y++)
        {
                /* Test this point only if it hasn't been already
                   tested and found to be free of obstructions. */

                if ((GetMask(path.lat[y],path.lon[y])&mask_value)==0)
                {
                        distance=5280.0*path.distance[y];
                        tx_alt=earthradius+source.alt+path.elevation[0];
                        rx_alt=earthradius+destination.alt+path.elevation[y];

                        /* Calculate the cosine of the elevation of the
                           transmitter as seen at the temp rx point. */

                        cos_xmtr_angle=((rx_alt*rx_alt)+(distance*distance)-(tx_alt*tx_alt))/(2.0*rx_alt*distance);

                        for (x=y, block=0; x>=0 && block==0; x--)
                        {
                                distance=5280.0*(path.distance[y]-path.distance[x]);
                                test_alt=earthradius+path.elevation[x];

                                cos_test_angle=((rx_alt*rx_alt)+(distance*distance)-(test_alt*test_alt))/(2.0*rx_alt*distance);

                                /* Compare these two angles to determine if
                                   an obstruction exists.  Since we're comparing
                                   the cosines of these angles rather than
                                   the angles themselves, the following "if"
                                   statement is reversed from what it would
                                   be if the actual angles were compared. */

                                if (cos_xmtr_angle>cos_test_angle)
                                        block=1;
                        }

                        if (block==0)
                                OrMask(path.lat[y],path.lon[y],mask_value);
                }
        }
}

void PlotLRPath(struct site source, struct site destination)
{
        /* This function plots the RF signal path loss
           between source and destination points based
           on the Longley-Rice propagation model. */

        char strmode[100];
        int x, y, errnum;
        double loss;

        ReadPath(source,destination);
        elev_l[1]=0.04*METERS_PER_MILE;

        for (x=0; x<path.length; x++)
                elev_l[x+2]=path.elevation[x]*METERS_PER_FOOT;  
         
        for (y=0; y<path.length; y++)
        {
                /* Test this point only if it has not already been tested. */

                if (GetMask(path.lat[y],path.lon[y])==0 && 0.04*y<=max_range)
                {
                        elev_l[0]=y+1;

                        point_to_point(elev_l,source.alt*METERS_PER_FOOT, 
                        destination.alt*METERS_PER_FOOT,
                        LR.eps_dielect, LR.sgm_conductivity, LR.eno_ns_surfref,
                        LR.frq_mhz, LR.radio_climate, LR.pol, LR.conf, LR.rel,
                        loss, strmode, errnum);

                        /* Note: PASS BY REFERENCE ... loss and errnum are
                           passed by reference, and are only used in this
                           file by this function */

                        if (loss>225.0)
                                loss=225.0;

                        if (loss<75.0)
                                loss=75.0;

                        loss-=75.0;
                        loss/=10.0;
                        loss+=1.0;
                
                        OrMask(path.lat[y],path.lon[y],((unsigned char)(loss))<<3);
                }

                else if (GetMask(path.lat[y],path.lon[y])==0 && 0.04*y>max_range)
                        OrMask(path.lat[y],path.lon[y],1);
        }
}

void PlotCoverage(struct site source, double altitude)
{
        /* This function performs a 360 degree sweep around the
           transmitter site (source location), and plots the
           line-of-sight coverage of the transmitter on the SPLAT!
           generated topographic map based on a receiver located
           at the specified altitude (in feet AGL).  Results
           are stored in memory, and written out in the form
           of a topographic map when the WritePPM() function
           is later invoked. */

        float lat, lon, one_pixel;
        static unsigned char mask_value;
        int z, count;
        struct site edge;
        unsigned char symbol[4], x;

        /* Initialize mask_value */

        if (mask_value!=8 && mask_value!=16 && mask_value!=32)
                mask_value=1;

        one_pixel=1.0/1200.0;

        symbol[0]='.';
        symbol[1]='o';
        symbol[2]='O';
        symbol[3]='o';

        count=0;        

        fprintf(stdout,"\nComputing line-of-sight coverage of %s with an RX antenna\nat %.2f feet AGL:\n\n 0%c to  25%c ",source.name,altitude,37,37);
        fflush(stdout);

        /* 18.75=1200 pixels/degree divided by 64 loops
           per progress indicator symbol (.oOo) printed. */

        z=(int)(18.75*ReduceAngle(max_west-min_west));

        for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
        {
                if (lon>=360.0)
                        lon-=360.0;

                edge.lat=max_north;
                edge.lon=lon;
                edge.alt=altitude;

                PlotPath(source,edge,mask_value);
                count++;

                if (count==z) 
                {
                        fprintf(stdout,"%c",symbol[x]);
                        fflush(stdout);
                        count=0;

                        if (x==3)
                                x=0;
                        else
                                x++;
                }
        }

        count=0;
        fprintf(stdout,"\n25%c to  50%c ",37,37);
        fflush(stdout);
        
        z=(int)(18.75*(max_north-min_north));

        for (lat=max_north, x=0; lat>=min_north; lat-=one_pixel)
        {
                edge.lat=lat;
                edge.lon=min_west;
                edge.alt=altitude;

                PlotPath(source,edge,mask_value);
                count++;

                if (count==z) 
                {
                        fprintf(stdout,"%c",symbol[x]);
                        fflush(stdout);
                        count=0;

                        if (x==3)
                                x=0;
                        else
                                x++;
                }
        }

        count=0;
        fprintf(stdout,"\n50%c to  75%c ",37,37);
        fflush(stdout);

        z=(int)(18.75*ReduceAngle(max_west-min_west));

        for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
        {
                if (lon>=360.0)
                        lon-=360.0;

                edge.lat=min_north;
                edge.lon=lon;
                edge.alt=altitude;

                PlotPath(source,edge,mask_value);
                count++;

                if (count==z)
                {
                        fprintf(stdout,"%c",symbol[x]);
                        fflush(stdout);
                        count=0;

                        if (x==3)
                                x=0;
                        else
                                x++;
                }
        }

        count=0;
        fprintf(stdout,"\n75%c to 100%c ",37,37);
        fflush(stdout);
        
        z=(int)(18.75*(max_north-min_north));

        for (lat=min_north, x=0; lat<=max_north; lat+=one_pixel)
        {
                edge.lat=lat;
                edge.lon=max_west;
                edge.alt=altitude;

                PlotPath(source,edge,mask_value);
                count++;

                if (count==z)
                {
                        fprintf(stdout,"%c",symbol[x]);
                        fflush(stdout);
                        count=0;

                        if (x==3)
                                x=0;
                        else
                                x++;
                }
        }

        fprintf(stdout,"\nDone!\n");
        fflush(stdout);

        /* Assign next mask value */

        switch (mask_value)
        {
                case 1:
                        mask_value=8;
                        break;

                case 8:
                        mask_value=16;
                        break;

                case 16:
                        mask_value=32;
        }
}

void PlotLRMap(struct site source, double altitude)
{
        /* This function performs a 360 degree sweep around the
           transmitter site (source location), and plots the
           Longley-Rice attenuation on the SPLAT! generated
           topographic map based on a receiver located at
           the specified altitude (in feet AGL).  Results
           are stored in memory, and written out in the form
           of a topographic map when the WritePPMLR() function
           is later invoked. */

        int z, count;
        struct site edge;
        float lat, lon, one_pixel;
        unsigned char symbol[4], x;

        one_pixel=1.0/1200.0;

        symbol[0]='.';
        symbol[1]='o';
        symbol[2]='O';
        symbol[3]='o';

        count=0;

        fprintf(stdout,"\nComputing Longley-Rice coverage of %s ", source.name);
        fprintf(stdout,"out to a radius\nof %.2f miles with an RX antenna at %.2f feet AGL:\n\n 0%c to  25%c ",max_range,altitude,37,37);
        fflush(stdout);

        /* 18.75=1200 pixels/degree divided by 64 loops
           per progress indicator symbol (.oOo) printed. */

        z=(int)(18.75*ReduceAngle(max_west-min_west));

        for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
        {
                if (lon>=360.0)
                        lon-=360.0;

                edge.lat=max_north;
                edge.lon=lon;
                edge.alt=altitude;

                PlotLRPath(source,edge);
                count++;

                if (count==z) 
                {
                        fprintf(stdout,"%c",symbol[x]);
                        fflush(stdout);
                        count=0;

                        if (x==3)
                                x=0;
                        else
                                x++;
                }
        }

        count=0;
        fprintf(stdout,"\n25%c to  50%c ",37,37);
        fflush(stdout);
        
        z=(int)(18.75*(max_north-min_north));

        for (lat=max_north, x=0; lat>=min_north; lat-=one_pixel)
        {
                edge.lat=lat;
                edge.lon=min_west;
                edge.alt=altitude;

                PlotLRPath(source,edge);
                count++;

                if (count==z) 
                {
                        fprintf(stdout,"%c",symbol[x]);
                        fflush(stdout);
                        count=0;

                        if (x==3)
                                x=0;
                        else
                                x++;
                }
        }

        count=0;
        fprintf(stdout,"\n50%c to  75%c ",37,37);
        fflush(stdout);

        z=(int)(18.75*ReduceAngle(max_west-min_west));

        for (lon=min_west, x=0; (LonDiff(lon,max_west)<=0.0); lon+=one_pixel)
        {
                if (lon>=360.0)
                        lon-=360.0;

                edge.lat=min_north;
                edge.lon=lon;
                edge.alt=altitude;

                PlotLRPath(source,edge);
                count++;

                if (count==z)
                {
                        fprintf(stdout,"%c",symbol[x]);
                        fflush(stdout);
                        count=0;

                        if (x==3)
                                x=0;
                        else
                                x++;
                }
        }

        count=0;
        fprintf(stdout,"\n75%c to 100%c ",37,37);
        fflush(stdout);
        
        z=(int)(18.75*(max_north-min_north));

        for (lat=min_north, x=0; lat<=max_north; lat+=one_pixel)
        {
                edge.lat=lat;
                edge.lon=max_west;
                edge.alt=altitude;

                PlotLRPath(source,edge);
                count++;

                if (count==z)
                {
                        fprintf(stdout,"%c",symbol[x]);
                        fflush(stdout);
                        count=0;

                        if (x==3)
                                x=0;
                        else
                                x++;
                }
        }

        fprintf(stdout,"\nDone!\n");
        fflush(stdout);
}

void WritePPM(char *filename)
{
        /* This function generates a topographic map in Portable Pix Map
           (PPM) format based on logarithmically scaled topology data,
           as well as the content of flags held in the mask[][] array.
           The image created is rotated counter-clockwise 90 degrees
           from its representation in dem[][] so that north points
           up and east points right in the image generated. */

        char mapfile[255];
        unsigned char found, mask;
        unsigned width, height, output;
        int indx, x, y, x0=0, y0=0, minlat, minlon;
        float lat, lon, one_pixel, conversion, one_over_gamma;
        FILE *fd;

        one_pixel=1.0/1200.0;
        one_over_gamma=1.0/GAMMA;
        conversion=255.0/pow((double)(max_elevation-min_elevation),one_over_gamma);

        width=(unsigned)(1200*ReduceAngle(max_west-min_west));
        height=(unsigned)(1200*ReduceAngle(max_north-min_north));

        if (filename[0]==0)
                strncpy(mapfile, "map.ppm\0",8);
        else
        {
                for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
                        mapfile[x]=filename[x];

                mapfile[x]='.';
                mapfile[x+1]='p';
                mapfile[x+2]='p';
                mapfile[x+3]='m';
                mapfile[x+4]=0;
        }

        fd=fopen(mapfile,"wb");

        fprintf(fd,"P6\n%u %u\n255\n",width,height);

        fprintf(stdout,"\nWriting \"%s\" (%ux%u pixmap image)... ",mapfile,width,height);
        fflush(stdout);

        for (y=0, lat=((double)max_north)-one_pixel; y<(int)height; y++, lat-=one_pixel)
        {
                minlat=(int)floor(lat);

                for (x=0, lon=((double)max_west)-one_pixel; x<(int)width; x++, lon-=one_pixel)
                {
                        if (lon<0.0)
                                lon+=360.0;

                        minlon=(int)floor(lon);

                        for (indx=0, found=0; indx<MAXSLOTS && found==0;)
                                if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
                                        found=1;
                                else
                                        indx++;

                        if (found)
                        {
                                x0=(int)(1199.0*(lat-floor(lat)));
                                y0=(int)(1199.0*(lon-floor(lon)));

                                mask=dem[indx].mask[x0][y0];

                                if (mask&2)
                                        /* Text Labels: Red */
                                        fprintf(fd,"%c%c%c",255,0,0);

                                else if (mask&4)
                                        /* County Boundaries: Light Cyan */
                                        fprintf(fd,"%c%c%c",128,128,255);

                                else switch (mask&57)
                                {
                                        case 1:
                                        /* TX1: Green */
                                        fprintf(fd,"%c%c%c",0,255,0);
                                        break;

                                        case 8:
                                        /* TX2: Cyan */
                                        fprintf(fd,"%c%c%c",0,255,255);
                                        break;

                                        case 9:
                                        /* TX1 + TX2: Yellow */
                                        fprintf(fd,"%c%c%c",255,255,0);
                                        break;

                                        case 16:
                                        /* TX3: Medium Violet */
                                        fprintf(fd,"%c%c%c",147,112,219);
                                        break;

                                        case 17:
                                        /* TX1 + TX3: Pink */
                                        fprintf(fd,"%c%c%c",255,192,203);
                                        break;

                                        case 24:
                                        /* TX2 + TX3: Orange */
                                        fprintf(fd,"%c%c%c",255,165,0);
                                        break;

                                        case 25:
                                        /* TX1 + TX2 + TX3: Dark Green */
                                        fprintf(fd,"%c%c%c",0,100,0);
                                        break;

                                        case 32:
                                        /* TX4: Sienna 1 */
                                        fprintf(fd,"%c%c%c",255,130,71);
                                        break;

                                        case 33:
                                        /* TX1 + TX4: Green Yellow */
                                        fprintf(fd,"%c%c%c",173,255,47);
                                        break;

                                        case 40:
                                        /* TX2 + TX4: Dark Sea Green 1 */
                                        fprintf(fd,"%c%c%c",193,255,193);
                                        break;

                                        case 41:
                                        /* TX1 + TX2 + TX4: Blanched Almond */
                                        fprintf(fd,"%c%c%c",255,235,205);
                                        break;

                                        case 48:
                                        /* TX3 + TX4: Dark Turquoise */
                                        fprintf(fd,"%c%c%c",0,206,209);
                                        break;

                                        case 49:
                                        /* TX1 + TX3 + TX4: Medium Spring Green */
                                        fprintf(fd,"%c%c%c",0,250,154);
                                        break;

                                        case 56:
                                        /* TX2 + TX3 + TX4: Tan */
                                        fprintf(fd,"%c%c%c",210,180,140);
                                        break;

                                        case 57:
                                        /* TX1 + TX2 + TX3 + TX4: Gold2 */
                                        fprintf(fd,"%c%c%c",238,201,0);
                                        break;

                                        default:
                                        /* Water: Medium Blue */
                                        if (dem[indx].data[x0][y0]==0)
                                                fprintf(fd,"%c%c%c",0,0,170);
                                        else
                                        {
                                                /* Elevation: Greyscale */
                                                output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
                                                fprintf(fd,"%c%c%c",output,output,output);
                                        }
                                }
                        }

                        else
                        {
                                /* We should never get here, but if */
                                /* we do, display the region as black */

                                fprintf(fd,"%c%c%c",0,0,0);
                        }
                }
        }

        fclose(fd);
        fprintf(stdout,"Done!\n");
        fflush(stdout);
}

void WritePPMLR(char *filename)
{
        /* This function generates a topographic map in Portable Pix Map
           (PPM) format based on the content of flags held in the mask[][] 
           array (only).  The image created is rotated counter-clockwise
           90 degrees from its representation in dem[][] so that north
           points up and east points right in the image generated. */

        char mapfile[255];
        unsigned width, height, output;
        unsigned char found, mask, cityorcounty;
        int indx, x, y, t, t2, x0, y0, minlat, minlon, loss;
        float lat, lon, one_pixel, conversion, one_over_gamma;
        FILE *fd;

        one_pixel=1.0/1200.0;
        one_over_gamma=1.0/GAMMA;
        conversion=255.0/pow((double)(max_elevation-min_elevation),one_over_gamma);

        width=(unsigned)(1200*ReduceAngle(max_west-min_west));
        height=(unsigned)(1200*ReduceAngle(max_north-min_north));

        if (filename[0]==0)
                strncpy(mapfile, "map.ppm\0",8);
        else
        {
                for (x=0; filename[x]!='.' && filename[x]!=0 && x<250; x++)
                        mapfile[x]=filename[x];

                mapfile[x]='.';
                mapfile[x+1]='p';
                mapfile[x+2]='p';
                mapfile[x+3]='m';
                mapfile[x+4]=0;
        }

        fd=fopen(mapfile,"wb");

        fprintf(fd,"P6\n%u %u\n255\n",width,height+30);

        fprintf(stdout,"\nWriting \"%s\" (%ux%u pixmap image)... ",mapfile,width,height+30);
        fflush(stdout);

        for (y=0, lat=((double)max_north)-one_pixel; y<(int)height; y++, lat-=one_pixel)
        {
                minlat=(int)floor(lat);

                for (x=0, lon=((double)max_west)-one_pixel; x<(int)width; x++, lon-=one_pixel)
                {
                        if (lon<0.0)
                                lon+=360.0;

                        minlon=(int)floor(lon);

                        for (indx=0, found=0; indx<MAXSLOTS && found==0;)
                                if (minlat==dem[indx].min_north && minlon==dem[indx].min_west)
                                        found=1;
                                else
                                        indx++;
                        if (found)
                        {
                                x0=(int)(1199.0*(lat-floor(lat)));
                                y0=(int)(1199.0*(lon-floor(lon)));

                                mask=dem[indx].mask[x0][y0];
                                loss=70+(10*(int)((mask&248)>>3));
                                cityorcounty=0;

                                if (mask&2)
                                {
                                        /* Text Labels - Black or Red */

                                        /* if ((mask&120) && (loss<=maxdB)) */
                                        if ((mask&120) && (loss<=90))
                                                fprintf(fd,"%c%c%c",0,0,0);
                                        else
                                                fprintf(fd,"%c%c%c",255,0,0);

                                        cityorcounty=1;
                                }

                                else if (mask&4)
                                {
                                        /* County Boundaries: Black */

                                        fprintf(fd,"%c%c%c",0,0,0);

                                        cityorcounty=1;
                                }

                                if (cityorcounty==0)
                                {
                                        if (loss>maxdB)

                                        { /* Display land or sea elevation */

                                                if (dem[indx].data[x0][y0]==0)
                                                        fprintf(fd,"%c%c%c",0,0,170);
                                                else
                                                {
                                                        /* Elevation: Greyscale */
                                                        output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
                                                        fprintf(fd,"%c%c%c",output,output,output);
                                                }
                                        }

                                        else switch (loss)
                                        {
                                                /* Plot signal loss in color */

                                                case 80:
                                                fprintf(fd,"%c%c%c",255,0,0);
                                                break;

                                                case 90:
                                                fprintf(fd,"%c%c%c",255,128,0);
                                                break;

                                                case 100:
                                                fprintf(fd,"%c%c%c",255,165,0);
                                                break;

                                                case 110:
                                                fprintf(fd,"%c%c%c",255,206,0);
                                                break;

                                                case 120:
                                                fprintf(fd,"%c%c%c",255,255,0);
                                                break;

                                                case 130:
                                                fprintf(fd,"%c%c%c",184,255,0);
                                                break;

                                                case 140:
                                                fprintf(fd,"%c%c%c",0,255,0);
                                                break;

                                                case 150:
                                                fprintf(fd,"%c%c%c",0,208,0);
                                                break;

                                                case 160:
                                                fprintf(fd,"%c%c%c",0,196,196);
                                                break;

                                                case 170:
                                                fprintf(fd,"%c%c%c",0,148,255);
                                                break;

                                                case 180:
                                                fprintf(fd,"%c%c%c",80,80,255);
                                                break;

                                                case 190:
                                                fprintf(fd,"%c%c%c",0,38,255);
                                                break;

                                                case 200:
                                                fprintf(fd,"%c%c%c",142,63,255);
                                                break;

                                                case 210:
                                                fprintf(fd,"%c%c%c",196,54,255);
                                                break;

                                                case 220:
                                                fprintf(fd,"%c%c%c",255,0,255);
                                                break;

                                                case 230:
                                                fprintf(fd,"%c%c%c",255,194,204);
                                                break;

                                                default:

                                                if (dem[indx].data[x0][y0]==0)
                                                        fprintf(fd,"%c%c%c",0,0,170);
                                                else
                                                {
                                                        /* Elevation: Greyscale */
                                                        output=(unsigned)(0.5+pow((double)(dem[indx].data[x0][y0]-min_elevation),one_over_gamma)*conversion);
                                                        fprintf(fd,"%c%c%c",output,output,output);
                                                }
                                        }
                                }
                        }

                        else
                        {
                                /* We should never get here, but if */
                                /* we do, display the region as black */

                                fprintf(fd,"%c%c%c",0,0,0);
                        }
                }
        }

        /* Display legend along bottom of image */

        x0=width/16;

        for (y0=0; y0<30; y0++)
        {
                for (indx=0; indx<16; indx++)
                {
                        for (x=0; x<x0; x++)
                        {
                                t=indx;  
                                t2=indx+8;

                                if (y0>=10 && y0<=18)
                                {  
                                        if (t2>9)
                                        {
                                                if (x>=11 && x<=17)     
                                                        if (smallfont[t2/10][y0-10][x-11])
                                                                t=255; 
                                        }

                                        if (x>=19 && x<=25)     
                                                if (smallfont[t2%10][y0-10][x-19])
                                                        t=255;
 
                                        if (x>=27 && x<=33)
                                                if (smallfont[0][y0-10][x-27])
                                                        t=255; 
                                }

                                switch (t)
                                {
                                        case 0:
                                        fprintf(fd,"%c%c%c",255,0,0);
                                        break;

                                        case 1:
                                        fprintf(fd,"%c%c%c",255,128,0);
                                        break;

                                        case 2:
                                        fprintf(fd,"%c%c%c",255,165,0);
                                        break;

                                        case 3:
                                        fprintf(fd,"%c%c%c",255,206,0);
                                        break;

                                        case 4:
                                        fprintf(fd,"%c%c%c",255,255,0);
                                        break;

                                        case 5:
                                        fprintf(fd,"%c%c%c",184,255,0);
                                        break;

                                        case 6:
                                        fprintf(fd,"%c%c%c",0,255,0);
                                        break;

                                        case 7:
                                        fprintf(fd,"%c%c%c",0,208,0);
                                        break;

                                        case 8:
                                        fprintf(fd,"%c%c%c",0,196,196);
                                        break;

                                        case 9:
                                        fprintf(fd,"%c%c%c",0,148,255);
                                        break;

                                        case 10:
                                        fprintf(fd,"%c%c%c",80,80,255);
                                        break;

                                        case 11:
                                        fprintf(fd,"%c%c%c",0,38,255);
                                        break;

                                        case 12:
                                        fprintf(fd,"%c%c%c",142,63,255);
                                        break;

                                        case 13:
                                        fprintf(fd,"%c%c%c",196,54,255);
                                        break;

                                        case 14:
                                        fprintf(fd,"%c%c%c",255,0,255);
                                        break;

                                        case 255:
                                        /* Black */
                                        fprintf(fd,"%c%c%c",0,0,0);
                                        break;

                                        default:
                                        fprintf(fd,"%c%c%c",255,194,204);
                                }
                        } 
                }
        }

        fclose(fd);
        fprintf(stdout,"Done!\n");
        fflush(stdout);
}

void GraphTerrain(struct site source, struct site destination, char *name)
{
        /* This function invokes gnuplot to generate an appropriate
           output file indicating the terrain profile between the source
           and destination locations.  "filename" is the name assigned
           to the output file generated by gnuplot.  The filename extension
           is used to set gnuplot's terminal setting and output file type.
           If no extension is found, .gif is assumed.  */

        int x, y, z;
        char filename[255], term[15], ext[15];
        FILE *fd=NULL;

        ReadPath(destination,source);

        fd=fopen("profile.gp","wb");

        for (x=0; x<path.length; x++)
                fprintf(fd,"%f\t%f\n",path.distance[x],path.elevation[x]);

        fclose(fd);

        if (name[0]==0)
        {
                /* Default filename and output file type */

                strncpy(filename,"profile\0",8);
                strncpy(term,"gif\0",4);
                strncpy(ext,"gif\0",4);
        }

        else
        {
                /* Grab extension and terminal type from "name" */

                for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
                        filename[x]=name[x];

                if (name[x]=='.')
                {
                        for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
                        {
                                term[y]=tolower(name[x]);
                                ext[y]=term[y];
                        }

                        ext[y]=0;
                        term[y]=0;
                        filename[z]=0;
                }

                else
                {       /* No extension -- Default is gif */

                        filename[x]=0;
                        strncpy(term,"gif\0",4);
                        strncpy(ext,"gif\0",4);
                }
        }

        /* Either .ps or .postscript may be used
           as an extension for postscript output. */

        if (strncmp(term,"postscript",10)==0)
                strncpy(ext,"ps\0",3);

        else if (strncmp(ext,"ps",2)==0)
                strncpy(term,"postscript\0",11);

        fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
        fflush(stdout);

        fd=fopen("splat.gp","w");
        fprintf(fd,"set grid\n");
        fprintf(fd,"set autoscale\n");
        fprintf(fd,"set term %s\n",term);
        fprintf(fd,"set title \"SPLAT! Terrain Profile\"\n");
        fprintf(fd,"set xlabel \"Distance Between %s and %s (miles)\"\n",destination.name,source.name);
        fprintf(fd,"set ylabel \"Ground Elevation Above Sea Level (feet)\"\n");
        fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
        fprintf(fd,"plot \"profile.gp\" title \"\" with lines\n");
        fclose(fd);
                        
        x=system("gnuplot splat.gp");

        if (x!=-1)
        {
                unlink("splat.gp");
                unlink("profile.gp");
                fprintf(stdout," Done!\n");
                fflush(stdout);
        }

        else
                fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
}

void GraphElevation(struct site source, struct site destination, char *name)
{
        /* This function invokes gnuplot to generate an appropriate
           output file indicating the terrain profile between the source
           and destination locations.  "filename" is the name assigned
           to the output file generated by gnuplot.  The filename extension
           is used to set gnuplot's terminal setting and output file type.
           If no extension is found, .gif is assumed. */

        int x, y, z;
        char filename[255], term[15], ext[15];
        double angle, refangle, maxangle=-90.0;
        struct site remote;
        FILE *fd=NULL, *fd2=NULL;

        ReadPath(destination,source);  /* destination=RX, source=TX */
        refangle=ElevationAngle(destination,source);

        fd=fopen("profile.gp","wb");
        fd2=fopen("reference.gp","wb");

        for (x=1; x<path.length-1; x++)
        {
                remote.lat=path.lat[x];
                remote.lon=path.lon[x];
                remote.alt=0.0;
                angle=ElevationAngle(destination,remote);
                fprintf(fd,"%f\t%f\n",path.distance[x],angle);
                fprintf(fd2,"%f\t%f\n",path.distance[x],refangle);

                if (angle>maxangle)
                        maxangle=angle;
        }

        fprintf(fd,"%f\t%f\n",path.distance[path.length-1],refangle);
        fprintf(fd2,"%f\t%f\n",path.distance[path.length-1],refangle);

        fclose(fd);
        fclose(fd2);

        if (name[0]==0)
        {
                /* Default filename and output file type */

                strncpy(filename,"profile\0",8);
                strncpy(term,"gif\0",4);
                strncpy(ext,"gif\0",4);
        }

        else
        {
                /* Grab extension and terminal type from "name" */

                for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
                        filename[x]=name[x];

                if (name[x]=='.')
                {
                        for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
                        {
                                term[y]=tolower(name[x]);
                                ext[y]=term[y];
                        }

                        ext[y]=0;
                        term[y]=0;
                        filename[z]=0;
                }

                else
                {       /* No extension -- Default is gif */

                        filename[x]=0;
                        strncpy(term,"gif\0",4);
                        strncpy(ext,"gif\0",4);
                }
        }

        /* Either .ps or .postscript may be used
           as an extension for postscript output. */

        if (strncmp(term,"postscript",10)==0)
                strncpy(ext,"ps\0",3);

        else if (strncmp(ext,"ps",2)==0)
                strncpy(term,"postscript\0",11);

        fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
        fflush(stdout);

        fd=fopen("splat.gp","w");

        fprintf(fd,"set grid\n");
        fprintf(fd,"set yrange [%2.3f to %2.3f]\n", (-fabs(refangle)-0.25), maxangle+0.25);
        fprintf(fd,"set term %s\n",term);
        fprintf(fd,"set title \"SPLAT! Elevation Profile\"\n");
        fprintf(fd,"set xlabel \"Distance Between %s and %s (miles)\"\n",destination.name,source.name);
        fprintf(fd,"set ylabel \"Elevation Angle Along Path Between %s and %s (degrees)\"\n",destination.name,source.name);
        fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
        fprintf(fd,"plot \"profile.gp\" title \"Real Earth Profile\" with lines, \"reference.gp\" title \"Line Of Sight Path\" with lines\n");

        fclose(fd);
                        
        x=system("gnuplot splat.gp");

        if (x!=-1)
        {
                unlink("splat.gp");
                unlink("profile.gp");
                unlink("reference.gp"); 

                fprintf(stdout," Done!\n");
                fflush(stdout);
        }

        else
                fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
}

void GraphHeight(struct site source, struct site destination, char *name)
{
        /* This function invokes gnuplot to generate an appropriate
           output file indicating the terrain profile between the source
           and destination locations.  What is plotted is the height of
           land above or below a straight line between the receibe and
           transmit sites.  "filename" is the name assigned to the output
           file generated by gnuplot.  The filename extension is used
           to set gnuplot's terminal setting and output file type.
           If no extension is found, .gif is assumed. */

        int x, y, z;
        char filename[255], term[15], ext[15];
        double a, b, c, height, refangle, cangle, maxheight=-100000.0,
               minheight=100000.0;
        struct site remote;
        FILE *fd=NULL, *fd2=NULL;

        ReadPath(destination,source);  /* destination=RX, source=TX */
        refangle=ElevationAngle(destination,source);
        b=GetElevation(destination)+destination.alt+earthradius;

        fd=fopen("profile.gp","wb");
        fd2=fopen("reference.gp","wb");

        for (x=1; x<path.length-1; x++)
        {
                remote.lat=path.lat[x];
                remote.lon=path.lon[x];
                remote.alt=0.0;

                a=GetElevation(remote)+earthradius;

                cangle=5280.0*Distance(destination,remote)/earthradius;

                c=b*sin(refangle*deg2rad+HALFPI)/sin(HALFPI-refangle*deg2rad-cangle);

                height=a-c;

                fprintf(fd,"%f\t%f\n",path.distance[x],height);
                fprintf(fd2,"%f\t%f\n",path.distance[x],0.);

                if (height>maxheight)
                        maxheight=height;

                if (height<minheight)
                        minheight=height;
        }

        fprintf(fd,"%f\t%f\n",path.distance[path.length-1],0.0);
        fprintf(fd2,"%f\t%f\n",path.distance[path.length-1],0.0);

        fclose(fd);
        fclose(fd2);

        if (name[0]==0)
        {
                /* Default filename and output file type */

                strncpy(filename,"height\0",8);
                strncpy(term,"gif\0",4);
                strncpy(ext,"gif\0",4);
        }

        else
        {
                /* Grab extension and terminal type from "name" */

                for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
                        filename[x]=name[x];

                if (name[x]=='.')
                {
                        for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
                        {
                                term[y]=tolower(name[x]);
                                ext[y]=term[y];
                        }

                        ext[y]=0;
                        term[y]=0;
                        filename[z]=0;
                }

                else
                {       /* No extension -- Default is gif */

                        filename[x]=0;
                        strncpy(term,"gif\0",4);
                        strncpy(ext,"gif\0",4);
                }
        }

        /* Either .ps or .postscript may be used
           as an extension for postscript output. */

        if (strncmp(term,"postscript",10)==0)
                strncpy(ext,"ps\0",3);

        else if (strncmp(ext,"ps",2)==0)
                strncpy(term,"postscript\0",11);

        fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
        fflush(stdout);

        fd=fopen("splat.gp","w");

        minheight-=20.0;
        maxheight+=20.0;

        if (maxheight<20.0)
                maxheight=20.0;

        fprintf(fd,"set grid\n");
        fprintf(fd,"set yrange [%2.3f to %2.3f]\n", minheight, maxheight);
        fprintf(fd,"set term %s\n",term);
        fprintf(fd,"set title \"SPLAT! Height Profile\"\n");
        fprintf(fd,"set xlabel \"Distance Between %s and %s (miles)\"\n",destination.name,source.name);
        fprintf(fd,"set ylabel \"Ground Height Above Path Between %s and %s (feet)\"\n",destination.name,source.name);
        fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
        fprintf(fd,"plot \"profile.gp\" title \"Real Earth Profile\" with lines, \"reference.gp\" title \"Line Of Sight Path\" with lines\n");

        fclose(fd);

        x=system("gnuplot splat.gp");

        if (x!=-1)
        {
                unlink("splat.gp");
                unlink("profile.gp");
                unlink("reference.gp"); 
                fprintf(stdout," Done!\n");
                fflush(stdout);
        }

        else
                fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
}

void GraphLongley(struct site source, struct site destination, char *name)
{
        /* This function invokes gnuplot to generate an appropriate
           output file indicating the Longley-Rice model loss between
           the source and destination locations.   "filename" is the
           name assigned to the output file generated by gnuplot.
           The filename extension is used to set gnuplot's terminal
           setting and output file type.  If no extension is found,
           .gif is assumed.  */

        int x, y, z, errnum, errflag=0;
        char filename[255], term[15], ext[15], strmode[100], report_name[80];
        double maxloss=-100000.0, minloss=100000.0, loss, haavt, angle;
        FILE *fd=NULL, *fd2=NULL;

        sprintf(report_name,"%s-to-%s.lro",source.name,destination.name);

        for (x=0; report_name[x]!=0; x++)
                if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
                        report_name[x]='_';     

        fd2=fopen(report_name,"w");

        fprintf(fd2,"\n\t--==[ SPLAT! v%s Longley-Rice Model Path Loss Report ]==--\n\n",splat_version);
        fprintf(fd2,"Analysis of RF path conditions between %s and %s:\n",source.name, destination.name);
        fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");
        fprintf(fd2,"Transmitter site: %s\n",source.name);

        if (source.lat>=0.0)
        {
                fprintf(fd2,"Site location: %.4f North / %.4f West",source.lat, source.lon);
                fprintf(fd2, " (%s N / ", dec2dms(source.lat));
        }

        else
        {

                fprintf(fd2,"Site location: %.4f South / %.4f West",-source.lat, source.lon);
                fprintf(fd2, " (%s S / ", dec2dms(source.lat));
        }
        
        fprintf(fd2, "%s W)\n", dec2dms(source.lon));
        fprintf(fd2,"Ground elevation: %.2f feet AMSL\n",GetElevation(source));
        fprintf(fd2,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",source.alt, source.alt+GetElevation(source));

        haavt=haat(source);

        if (haavt>-4999.0)
                fprintf(fd2,"Antenna height above average terrain: %.2f feet\n",haavt);

        fprintf(fd2,"Distance to %s: %.2f miles.\n",destination.name,Distance(source,destination));
        fprintf(fd2,"Azimuth to %s: %.2f degrees.\n",destination.name,Azimuth(source,destination));

        angle=ElevationAngle(source,destination);

        if (angle>=0.0)
                fprintf(fd2,"Angle of elevation between %s and %s: %+.4f degrees.\n",source.name,destination.name,angle);

        if (angle<0.0)
                fprintf(fd2,"Angle of depression between %s and %s: %+.4f degrees.\n",source.name,destination.name,angle);

        fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");

        /* Receiver */

        fprintf(fd2,"Receiver site: %s\n",destination.name);

        if (destination.lat>=0.0)
        {
                fprintf(fd2,"Site location: %.4f North / %.4f West",destination.lat, destination.lon);
                fprintf(fd2, " (%s N / ", dec2dms(destination.lat));
        }

        else
        {
                fprintf(fd2,"Site location: %.4f South / %.4f West",-destination.lat, destination.lon);
                fprintf(fd2, " (%s S / ", dec2dms(destination.lat));
        }

        fprintf(fd2, "%s W)\n", dec2dms(destination.lon));
        fprintf(fd2,"Ground elevation: %.2f feet AMSL\n",GetElevation(destination));
        fprintf(fd2,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",destination.alt, destination.alt+GetElevation(destination));

        haavt=haat(destination);

        if (haavt>-4999.0)
                fprintf(fd2,"Antenna height above average terrain: %.2f feet\n",haavt);

        fprintf(fd2,"Distance to %s: %.2f miles.\n",source.name,Distance(source,destination));
        fprintf(fd2,"Azimuth to %s: %.2f degrees.\n",source.name,Azimuth(destination,source));

        angle=ElevationAngle(destination,source);

        if (angle>=0.0)
                fprintf(fd2,"Angle of elevation between %s and %s: %+.4f degrees.\n",destination.name,source.name,angle);

        if (angle<0.0)
                fprintf(fd2,"Angle of depression between %s and %s: %+.4f degrees.\n",destination.name,source.name,angle);

        fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");

        fprintf(fd2,"Longley-Rice path calculation parameters used in this analysis:\n\n");
        fprintf(fd2,"Earth's Dielectric Constant: %.3lf\n",LR.eps_dielect);
        fprintf(fd2,"Earth's Conductivity: %.3lf\n",LR.sgm_conductivity);
        fprintf(fd2,"Atmospheric Bending Constant (N): %.3lf\n",LR.eno_ns_surfref);
        fprintf(fd2,"Frequency: %.3lf (MHz)\n",LR.frq_mhz);
        fprintf(fd2,"Radio Climate: %d (",LR.radio_climate);

        switch (LR.radio_climate)
        {
                case 1:
                fprintf(fd2,"Equatorial");
                break;

                case 2:
                fprintf(fd2,"Continental Subtropical");
                break;

                case 3:
                fprintf(fd2,"Maritime Subtropical");
                break;

                case 4:
                fprintf(fd2,"Desert");
                break;

                case 5:
                fprintf(fd2,"Continental Temperate");
                break;

                case 6:
                fprintf(fd2,"Martitime Temperate, Over Land");
                break;

                case 7:
                fprintf(fd2,"Maritime Temperate, Over Sea");
                break;

                default:
                fprintf(fd2,"Unknown");
        }

        fprintf(fd2,")\nPolarization: %d (",LR.pol);

        if (LR.pol==0)
                fprintf(fd2,"Horizontal");

        if (LR.pol==1)
                fprintf(fd2,"Vertical");

        fprintf(fd2,")\nFraction of Situations: %.1lf%c\n",LR.conf*100.0,37);
        fprintf(fd2,"Fraction of Time: %.1lf%c\n",LR.rel*100.0,37);

        fprintf(fd2,"\n-------------------------------------------------------------------------\n\n");

        fprintf(fd2,"Analysis Results:\n\n");

        ReadPath(source, destination);  /* destination=RX, source=TX */

        elev_l[1]=0.04*METERS_PER_MILE;

        for (x=1; x<path.length; x++)
                elev_l[x+1]=path.elevation[x]*METERS_PER_FOOT;  

        fprintf(fd2,"Distance (mi)\tLoss (dB)\tErrnum\tComment\n\n"); 

        fd=fopen("profile.gp","w");

        for (x=2; x<path.length; x++)
        {
                elev_l[0]=x-1;

                point_to_point(elev_l, source.alt*METERS_PER_FOOT, 
                        destination.alt*METERS_PER_FOOT,
                        LR.eps_dielect, LR.sgm_conductivity, LR.eno_ns_surfref,
                        LR.frq_mhz, LR.radio_climate, LR.pol, LR.conf, LR.rel,
                        loss, strmode, errnum);

                        /* Note: PASS BY REFERENCE ... loss and errnum are
                           passed by reference, and are only used in this
                           file by this function */


                fprintf(fd,"%f\t%f\n",path.distance[path.length-1]-path.distance[x],loss);
                fprintf(fd2,"%7.2f\t\t%7.2f\t\t  %d\t%s\n",path.distance[x],loss, errnum, strmode);

                errflag|=errnum;
                  
                if (loss>maxloss)
                        maxloss=loss;

                if (loss<minloss)
                        minloss=loss;
        }

        fclose(fd);

        if (errflag)
        {
                fprintf(fd2,"\nNotes on \"errnum\"...\n\n");
                fprintf(fd2,"  0: No error.  :-)\n");
                fprintf(fd2,"  1: Warning!  Some parameters are nearly out of range.\n");
                fprintf(fd2,"     Results should be used with caution.\n");
                fprintf(fd2,"  2: Note: Default parameters have been substituted for impossible ones.\n");
                fprintf(fd2,"  3: Warning!  A combination of parameters is out of range.\n");
                fprintf(fd2,"     Results are probably invalid.\n");
                fprintf(fd2,"  Other: Warning!  Some parameters are out of range.\n");
                fprintf(fd2,"    Results are probably invalid.\n\nEnd of Report\n");
        }

        fprintf(stdout,"Longley-Rice Path Loss between %s and %s is %.2f db\n",source.name, destination.name, loss);
        fprintf(stdout,"Path Loss Report written to: \"%s\"\n",report_name);
        fflush(stdout);

        fclose(fd2);

        if (name[0]==0)
        {
                /* Default filename and output file type */

                strncpy(filename,"loss\0",5);
                strncpy(term,"gif\0",4);
                strncpy(ext,"gif\0",4);
        }

        else
        {
                /* Grab extension and terminal type from "name" */

                for (x=0; name[x]!='.' && name[x]!=0 && x<254; x++)
                        filename[x]=name[x];

                if (name[x]=='.')
                {
                        for (y=0, z=x, x++; name[x]!=0 && x<254 && y<14; x++, y++)
                        {
                                term[y]=tolower(name[x]);
                                ext[y]=term[y];
                        }

                        ext[y]=0;
                        term[y]=0;
                        filename[z]=0;
                }

                else
                {       /* No extension -- Default is gif */

                        filename[x]=0;
                        strncpy(term,"gif\0",4);
                        strncpy(ext,"gif\0",4);
                }
        }

        /* Either .ps or .postscript may be used
           as an extension for postscript output. */

        if (strncmp(term,"postscript",10)==0)
                strncpy(ext,"ps\0",3);

        else if (strncmp(ext,"ps",2)==0)
                strncpy(term,"postscript\0",11);

        fprintf(stdout,"Writing \"%s.%s\"...",filename,ext);
        fflush(stdout);

        fd=fopen("splat.gp","w");

        fprintf(fd,"set grid\n");
        fprintf(fd,"set yrange [%2.3f to %2.3f]\n", minloss, maxloss);
        fprintf(fd,"set term %s\n",term);
        fprintf(fd,"set title \"SPLAT! Loss Profile\"\n");
        fprintf(fd,"set xlabel \"Distance Between %s and %s (miles)\"\n",destination.name,source.name);
        fprintf(fd,"set ylabel \"Longley-Rice Loss (dB)\"\n");
        fprintf(fd,"set output \"%s.%s\"\n",filename,ext);
        fprintf(fd,"plot \"profile.gp\" title \"Longley-Rice Loss\" with lines\n");

        fclose(fd);
                        
        x=system("gnuplot splat.gp");

        if (x!=-1)
        {
                unlink("splat.gp");
                unlink("profile.gp");
                unlink("reference.gp"); 

                fprintf(stdout," Done!\n");
                fflush(stdout);
        }

        else
                fprintf(stderr,"\n*** ERROR: Error occurred invoking gnuplot!\n");
}

void ObstructionReport(struct site xmtr, struct site rcvr, char report)
{
        struct site result, result2, new_site;
        double angle, haavt;
        unsigned char block;
        char report_name[80], string[255];
        int x;
        FILE *fd;

        sprintf(report_name,"%s-to-%s.txt",xmtr.name,rcvr.name);

        for (x=0; report_name[x]!=0; x++)
                if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
                        report_name[x]='_';     

        fd=fopen(report_name,"w");

        fprintf(fd,"\n\t\t--==[ SPLAT! v%s Obstruction Report ]==--\n\n",splat_version);
        fprintf(fd,"Analysis of line-of-sight path conditions between %s and %s:\n",xmtr.name, rcvr.name);
        fprintf(fd,"\n-------------------------------------------------------------------------\n\n");
        fprintf(fd,"Transmitter site: %s\n",xmtr.name);


        if (xmtr.lat>=0.0)
        {
                fprintf(fd,"Site location: %.4f North / %.4f West",xmtr.lat, xmtr.lon);
                fprintf(fd, " (%s N / ", dec2dms(xmtr.lat));
        }

        else
        {
                fprintf(fd,"Site location: %.4f South / %.4f West",-xmtr.lat, xmtr.lon);
                fprintf(fd, " (%s S / ", dec2dms(xmtr.lat));
        }

        fprintf(fd, "%s W)\n", dec2dms(xmtr.lon));
        fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(xmtr));
        fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",xmtr.alt, xmtr.alt+GetElevation(xmtr));

        haavt=haat(xmtr);

        if (haavt>-4999.0)
                fprintf(fd,"Antenna height above average terrain: %.2f feet\n",haavt);

        fprintf(fd,"Distance to %s: %.2f miles.\n",rcvr.name,Distance(xmtr,rcvr));
        fprintf(fd,"Azimuth to %s: %.2f degrees.\n",rcvr.name,Azimuth(xmtr,rcvr));

        angle=ElevationAngle(xmtr,rcvr);

        if (angle>=0.0)
                fprintf(fd,"Angle of elevation between %s and %s: %+.4f degrees.\n",xmtr.name,rcvr.name,angle);

        if (angle<0.0)
                fprintf(fd,"Angle of depression between %s and %s: %+.4f degrees.\n",xmtr.name,rcvr.name,angle);

        fprintf(fd,"\n-------------------------------------------------------------------------\n\n");

        /* Receiver */

        fprintf(fd,"Receiver site: %s\n",rcvr.name);

        if (rcvr.lat>=0.0)
        {
                fprintf(fd,"Site location: %.4f North / %.4f West",rcvr.lat, rcvr.lon);
                fprintf(fd, " (%s N / ", dec2dms(rcvr.lat));
        }

        else
        {
                fprintf(fd,"Site location: %.4f South / %.4f West",-rcvr.lat, rcvr.lon);
                fprintf(fd, " (%s S / ", dec2dms(rcvr.lat));
        }

        fprintf(fd, "%s W)\n", dec2dms(rcvr.lon));
        fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(rcvr));
        fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",rcvr.alt, rcvr.alt+GetElevation(rcvr));

        haavt=haat(rcvr);

        if (haavt>-4999.0)
                fprintf(fd,"Antenna height above average terrain: %.2f feet\n",haavt);

        fprintf(fd,"Distance to %s: %.2f miles.\n",xmtr.name,Distance(xmtr,rcvr));
        fprintf(fd,"Azimuth to %s: %.2f degrees.\n",xmtr.name,Azimuth(rcvr,xmtr));

        angle=ElevationAngle(rcvr,xmtr);

        if (angle>=0.0)
                fprintf(fd,"Angle of elevation between %s and %s: %+.4f degrees.\n",rcvr.name,xmtr.name,angle);

        if (angle<0.0)
                fprintf(fd,"Angle of depression between %s and %s: %+.4f degrees.\n",rcvr.name,xmtr.name,angle);

        fprintf(fd,"\n-------------------------------------------------------------------------\n\n");

        if (report=='y')
        {
                /* Write an Obstruction Report */

                new_site=rcvr;
                result=los(xmtr,rcvr);
                result2=result;
                result2.alt-=1;
                block=result.name[0];

                if (block=='*')
                        fprintf(fd,"SPLAT! detected obstructions at:\n\n");

                while (block=='*')
                {
                        if (result.lat!=result2.lat || result.lon!=result2.lon || result.alt!=result2.alt)
                        {
                                if (result.lat>=0.0)
                                        fprintf(fd,"\t%.4f N, %.4f W, %5.2f miles, %6.2f feet AMSL.\n",result.lat, result.lon, Distance(rcvr,result), result.alt);
                                else

                                        fprintf(fd,"\t%.4f S, %.4f W, %5.2f miles, %6.2f feet AMSL.\n",-result.lat, result.lon, Distance(rcvr,result), result.alt);
                        }

                        result2=result;
                        new_site.alt+=1.0;

                        /* Can you hear me now? :-) */

                        result=los(xmtr,new_site);
                        block=result.name[0];
                }

                if (new_site.alt!=rcvr.alt)
                        sprintf(string,"\nAntenna at %s must be raised to at least %.2f feet AGL\nto clear all obstructions detected by SPLAT!\n\n",rcvr.name, new_site.alt);
                else
                        sprintf(string,"\nNo obstructions due to terrain were detected by SPLAT!\n\n");
        }

        fprintf(fd,"%s",string);

        fclose(fd);

        /* Display LOS status to terminal */

        fprintf(stdout,"%sObstruction report written to: \"%s\"\n",string,report_name);
        fflush(stdout);
}

void SiteReport(struct site xmtr)
{
        char report_name[80];
        double terrain;
        int x, azi;
        FILE *fd;

        sprintf(report_name,"%s-site_report.txt",xmtr.name);

        for (x=0; report_name[x]!=0; x++)
                if (report_name[x]==32 || report_name[x]==17 || report_name[x]==92 || report_name[x]==42 || report_name[x]==47)
                        report_name[x]='_';     

        fd=fopen(report_name,"w");

        fprintf(fd,"\n\t--==[ SPLAT! v%s Site Analysis Report For: %s ]==--\n\n",splat_version,xmtr.name);

        fprintf(fd,"---------------------------------------------------------------------------\n\n");

        if (xmtr.lat>=0.0)
        {
                fprintf(fd,"Site location: %.4f North / %.4f West",xmtr.lat, xmtr.lon);
                fprintf(fd, " (%s N / ",dec2dms(xmtr.lat));
        }

        else
        {
                fprintf(fd,"Site location: %.4f South / %.4f West",-xmtr.lat, xmtr.lon);
                fprintf(fd, " (%s S / ",dec2dms(xmtr.lat));
        }

        fprintf(fd, "%s W)\n",dec2dms(xmtr.lon));
        fprintf(fd,"Ground elevation: %.2f feet AMSL\n",GetElevation(xmtr));
        fprintf(fd,"Antenna height: %.2f feet AGL / %.2f feet AMSL\n",xmtr.alt, xmtr.alt+GetElevation(xmtr));

        terrain=haat(xmtr);

        if (terrain>-4999.0)
        {
                fprintf(fd,"Antenna height above average terrain: %.2f feet\n\n",terrain);

                /* Display the average terrain between 2 and 10 miles
                   from the transmitter site at azimuths of 0, 45, 90,
                   135, 180, 225, 270, and 315 degrees. */

                for (azi=0; azi<=315; azi+=45)
                {
                        fprintf(fd,"Average terrain at %3d degrees azimuth: ",azi);
                        terrain=AverageTerrain(xmtr,(double)azi,2.0,10.0);

                        if (terrain>-4999.0)
                                fprintf(fd,"%.2f feet AMSL\n",terrain);
                        else
                                fprintf(fd,"No terrain\n");
                }
        }

        fprintf(fd,"\n---------------------------------------------------------------------------\n\n");
        fclose(fd);
        fprintf(stdout,"\nSite analysis report written to: \"%s\"\n",report_name);
}

int main(char argc, char *argv[])
{
        int             x, y, ymin, ymax, width, z=0, min_lat, min_lon,
                        max_lat, max_lon, rxlat, rxlon, txlat, txlon,
                        west_min, west_max, north_min, north_max;

        unsigned char   coverage=0, LRmap=0, ext[20], terrain_plot=0,
                        elevation_plot=0, height_plot=0, 
                        longley_plot=0, cities=0, bfs=0, txsites=0,
                        count, report='y';
 
        char            mapfile[255], header[80], city_file[5][255], 
                        elevation_file[255], height_file[255], 
                        longley_file[255], terrain_file[255],
                        string[255], rxfile[255], *env=NULL,
                        txfile[255], map=0, boundary_file[5][255],
                        rxsite=0;

        double          altitude=0.0, altitudeLR=0.0, tx_range=0.0,
                        rx_range=0.0, deg_range=0.0, deg_limit,
                        deg_range_lon, er_mult;

        struct          site tx_site[4], rx_site;

        FILE            *fd;

        sprintf(header,"\n\t\t--==[ Welcome To SPLAT! v%s ]==--\n\n", splat_version);

        if (argc==1)
        {
                fprintf(stdout, "%sAvailable Options...\n\n\t -t txsite(s).qth (max of 4)\n\t -r rxsite.qth\n",header);
                fprintf(stdout,"\t -c plot coverage area(s) of TX(s) based on an RX antenna at X feet AGL\n");
                fprintf(stdout,"\t -L plot path loss map of TX based on an RX antenna at X feet AGL\n");
                fprintf(stdout,"\t -s filename(s) of city/site file(s) to import (max of 5)\n");
                fprintf(stdout,"\t -b filename(s) of cartographic boundary file(s) to import (max of 5)\n");
                fprintf(stdout,"\t -p filename of terrain profile graph to plot\n");
                fprintf(stdout,"\t -e filename of terrain elevation graph to plot\n");
                fprintf(stdout,"\t -h filename of terrain height graph to plot\n");
                fprintf(stdout,"\t -l filename of Longley-Rice graph to plot\n");
                fprintf(stdout,"\t -o filename of topographic map to generate (.ppm)\n");
                fprintf(stdout,"\t -d sdf file directory path (overrides path in ~/.splat_path file)\n");
                fprintf(stdout,"\t -n no analysis, brief report\n\t -N no analysis, no report\n");
                fprintf(stdout,"\t -m earth radius multiplier\n");
                fprintf(stdout,"\t -R modify default range for -c or -L (miles)\n");
                fprintf(stdout,"\t-db maximum loss contour to display on path loss maps (80-230 dB)\n\n");

                fprintf(stdout,"Type 'man splat', or see the documentation for more details.\n\n");
                fflush(stdout);
                return 1;
        }

        y=argc-1;

        rxfile[0]=0;
        txfile[0]=0;
        string[0]=0;
        mapfile[0]=0;
        elevation_file[0]=0;
        terrain_file[0]=0;
        sdf_path[0]=0;
        path.length=0;
        rx_site.lat=91.0;
        rx_site.lon=361.0;
        earthradius=EARTHRADIUS;

        for (x=0; x<4; x++)
        {
                tx_site[x].lat=91.0;
                tx_site[x].lon=361.0;
        }

        for (x=0; x<MAXSLOTS; x++)
        {
                dem[x].min_el=32768;
                dem[x].max_el=-32768;
                dem[x].min_north=90;
                dem[x].max_north=-90;
                dem[x].min_west=360;
                dem[x].max_west=-1;
        }


        /* Scan for command line arguments */

        for (x=1; x<=y; x++)
        {
                if (strcmp(argv[x],"-R")==0)
                {
                        z=x+1;

                        if (z<=y && argv[z][0] && argv[z][0]!='-')
                        {
                                sscanf(argv[z],"%lf",&max_range);

                                if (max_range<0.0)
                                        max_range=0.0;

                                if (max_range>1000.0)
                                        max_range=1000.0;
                        }                        
                }

                if (strcmp(argv[x],"-m")==0)
                {
                        z=x+1;

                        if (z<=y && argv[z][0] && argv[z][0]!='-')
                        {
                                sscanf(argv[z],"%lf",&er_mult);

                                if (er_mult<0.1)
                                        er_mult=1.0;

                                if (er_mult>1.0e6)
                                        er_mult=1.0e6;

                                earthradius*=er_mult;
                        }                        
                }

                if (strcmp(argv[x],"-o")==0)
                {
                        z=x+1;

                        if (z<=y && argv[z][0] && argv[z][0]!='-')
                                strncpy(mapfile,argv[z],253);
                        map=1;
                }

                if (strcmp(argv[x],"-c")==0)
                {
                        z=x+1;

                        if (z<=y && argv[z][0] && argv[z][0]!='-')
                        {
                                sscanf(argv[z],"%lf",&altitude);
                                coverage=1;
                        }
                }

                if (strcmp(argv[x],"-db")==0)
                {
                        z=x+1;

                        if (z<=y && argv[z][0] && argv[z][0]!='-')
                        {
                                sscanf(argv[z],"%d",&maxdB);

                                maxdB=abs(maxdB);

                                if (maxdB<80)
                                        maxdB=80;

                                if (maxdB>230)
                                        maxdB=230;
                        }                        
                }

                if (strcmp(argv[x],"-p")==0)
                { 
                        z=x+1;

                        if (z<=y && argv[z][0] && argv[z][0]!='-')
                        {
                                strncpy(terrain_file,argv[z],253);
                                terrain_plot=1;
                        }
                }

                if (strcmp(argv[x],"-e")==0)
                {
                        z=x+1;

                        if (z<=y && argv[z][0] && argv[z][0]!='-')
                        {
                                strncpy(elevation_file,argv[z],253);
                                elevation_plot=1;
                        }
                }

                if (strcmp(argv[x],"-h")==0)
                {
                        z=x+1;

                        if (z<=y && argv[z][0] && argv[z][0]!='-')
                        {
                                strncpy(height_file,argv[z],253);
                                height_plot=1;
                        }
                }

                if (strcmp(argv[x],"-n")==0)
                {
                        if (z<=y && argv[z][0] && argv[z][0]!='-')
                        {
                                report='n';
                                map=1;
                        }
                }

                if (strcmp(argv[x],"-N")==0)
                {
                        if (z<=y && argv[z][0] && argv[z][0]!='-');
                        {
                                report='N';
                                map=1;
                        }
                }

                if (strcmp(argv[x],"-d")==0)
                {
                        z=x+1;

                        if (z<=y && argv[z][0] && argv[z][0]!='-')
                                strncpy(sdf_path,argv[z],253);
                }

                if (strcmp(argv[x],"-t")==0)
                {
                        /* Read Transmitter Location */

                        z=x+1;

                        while (z<=y && argv[z][0] && argv[z][0]!='-' && txsites<4)
                        {
                                strncpy(txfile,argv[z],253);
                                tx_site[txsites]=LoadQTH(txfile);
                                txsites++;
                                z++;
                        }
                        z--;
                }

                if (strcmp(argv[x],"-L")==0)
                {
                        z=x+1;

                        if (z<=y && argv[z][0] && argv[z][0]!='-')
                        {
                                sscanf(argv[z],"%lf",&altitudeLR);

                                if (coverage)
                                        fprintf(stdout,"c and L are exclusive options, ignoring L.\n");
                                else
                                {
                                        LRmap=1;
                                        ReadLRParm(txfile);
                                } 
                        }
                }

                if (strcmp(argv[x],"-l")==0)
                {
                        z=x+1;

                        if (z<=y && argv[z][0] && argv[z][0]!='-')
                        {
                                strncpy(longley_file,argv[z],253);
                                longley_plot=1;
                                /* Doing this twice is harmless */
                                ReadLRParm(txfile);
                        }
                }

                if (strcmp(argv[x],"-r")==0)
                {
                        /* Read Receiver Location */

                        z=x+1;

                        if (z<=y && argv[z][0] && argv[z][0]!='-')
                        {
                                strncpy(rxfile,argv[z],253);
                                rx_site=LoadQTH(rxfile);
                                rxsite=1;
                        }
                }

                if (strcmp(argv[x],"-s")==0)
                {
                        /* Read city file(s) */

                        z=x+1;

                        while (z<=y && argv[z][0] && argv[z][0]!='-' && cities<5)
                        {
                                strncpy(city_file[cities],argv[z],253);
                                cities++;
                                z++;
                        }
                        z--;
                }

                if (strcmp(argv[x],"-b")==0)
                {
                        /* Read Boundary File(s) */

                        z=x+1;

                        while (z<=y && argv[z][0] && argv[z][0]!='-' && bfs<5)
                        {
                                strncpy(boundary_file[bfs],argv[z],253);
                                bfs++;
                                z++;
                        }
                        z--;
                }
        }

        /* Perform some error checking on the arguments
           and switches parsed from the command-line.
           If an error is encountered, print a message
           and exit gracefully. */

        if (txsites==0)
        {
                fprintf(stderr,"\n%c*** ERROR: No transmitter site(s) specified!\n\n",7);
                exit (-1);
        }

        for (x=0, y=0; x<txsites; x++)
        {
                if (tx_site[x].lat==91.0 && tx_site[x].lon==361.0)
                {
                        fprintf(stderr,"\n*** ERROR: Transmitter site #%d not found!",x+1);
                        y++;
                }
        }

        if (y)
        {
                fprintf(stderr,"%c\n\n",7);
                exit (-1);
        }

        if ((coverage+LRmap)==0 && rx_site.lat==91.0 && rx_site.lon==361.0)
        {
                fprintf(stderr,"\n%c*** ERROR: No receiver site found or specified!\n\n",7);
                exit (-1);
        }

        /* No errors were detected.  Whew!  :-) */

        /* If no SDF path was specified on the command line (-d), check
           for a path specified in the $HOME/.splat_path file.  If the
           file is not found, then sdf_path[] remains NULL, and the
           current working directory is assumed to contain the SDF
           files. */

        if (sdf_path[0]==0)
        {
                env=getenv("HOME");
                sprintf(string,"%s/.splat_path",env);
                fd=fopen(string,"r");

                if (fd!=NULL)
                {
                        fgets(string,253,fd);

                        /* Remove <CR> and/or <LF> from string */

                        for (x=0; string[x]!=13 && string[x]!=10 && string[x]!=0 && x<253; x++);
                        string[x]=0;

                        strncpy(sdf_path,string,253);

                        fclose(fd);
                }
        }

        /* Ensure a trailing '/' is present in sdf_path */

        if (sdf_path[0])
        {
                x=strlen(sdf_path);

                if (sdf_path[x-1]!='/' && x!=0)
                {
                        sdf_path[x]='/';
                        sdf_path[x+1]=0;
                }
        }

        fprintf(stdout,"%s",header);
        fflush(stdout);

        x=0;
        y=0;

        min_lat=90;
        max_lat=-90;

        min_lon=(int)floor(tx_site[0].lon);
        max_lon=(int)floor(tx_site[0].lon);

        for (y=0, z=0; z<txsites; z++)
        {
                txlat=(int)floor(tx_site[z].lat);
                txlon=(int)floor(tx_site[z].lon);

                if (txlat<min_lat)
                        min_lat=txlat;

                if (txlat>max_lat)
                        max_lat=txlat;

                if (LonDiff(txlon,min_lon)<0.0)
                        min_lon=txlon;

                if (LonDiff(txlon,max_lon)>0.0)
                        max_lon=txlon;
        }

        if (rxsite)
        {
                rxlat=(int)floor(rx_site.lat);
                rxlon=(int)floor(rx_site.lon);

                if (rxlat<min_lat)
                        min_lat=rxlat;

                if (rxlat>max_lat)
                        max_lat=rxlat;

                if (LonDiff(rxlon,min_lon)<0.0)
                        min_lon=rxlon;

                if (LonDiff(rxlon,max_lon)>0.0)
                        max_lon=rxlon;
        }


        /* Load the required SDF files */ 

        width=ReduceAngle(max_lon-min_lon);

        if ((max_lon-min_lon)<=180.0)
        {
                for (y=0; y<=width; y++)
                        for (x=min_lat; x<=max_lat; x++)
                        {
                                ymin=(int)(min_lon+(double)y);

                                while (ymin<0)
                                        ymin+=360;

                                while (ymin>=360)
                                        ymin-=360;

                                ymax=ymin+1;

                                while (ymax<0)
                                        ymax+=360;

                                while (ymax>=360)
                                        ymax-=360;

                                sprintf(string,"%d:%d:%d:%d",x, x+1, ymin, ymax);
                                LoadSDF(string);
                        }
        }

        else
        {
                for (y=0; y<=width; y++)
                        for (x=min_lat; x<=max_lat; x++)
                        {
                                ymin=max_lon+y;

                                while (ymin<0)
                                        ymin+=360;

                                while (ymin>=360)
                                        ymin-=360;
                                        
                                ymax=ymin+1;

                                while (ymax<0)
                                        ymax+=360;

                                while (ymax>=360)
                                        ymax-=360;

                                sprintf(string,"%d:%d:%d:%d",x, x+1, ymin, ymax);
                                LoadSDF(string);
                        }
        }

        if (coverage | LRmap)
        {
                if (LRmap)
                        txsites=1;

                for (z=0; z<txsites; z++)
                {
                        /* "Ball park" estimates used to load any additional
                           SDF files required to conduct this analysis. */

                        tx_range=sqrt(1.5*(tx_site[z].alt+GetElevation(tx_site[z])));

                        if (LRmap)
                                rx_range=sqrt(1.5*altitudeLR);
                        else
                                rx_range=sqrt(1.5*altitude);

                        /* deg_range determines the maximum
                           amount of topo data we read */

                        deg_range=(tx_range+rx_range)/69.0;

                        /* max_range sets the maximum size of the
                           analysis.  A small, non-zero amount can
                           be used to shrink the size of the analysis
                           and limit the amount of topo data read by
                           SPLAT!  A very large number will only increase
                           the width of the analysis, not the size of
                           the map. */

                        if (max_range==0.0)
                                max_range=tx_range+rx_range;

                        if (max_range<(tx_range+rx_range))
                                deg_range=max_range/69.0;

                        /* Prevent the demand for a really wide coverage
                           from allocating more slots than are available
                           in memory. */

                        switch (MAXSLOTS)
                        {
                                case 2: deg_limit=0.25;
                                        break;

                                case 4: deg_limit=0.5;
                                        break;

                                case 9: deg_limit=1.0;
                                        break;

                                case 16: deg_limit=2.0;
                                        break;

                                case 25: deg_limit=3.0;
                        }

                        if (tx_site[z].lat<70.0)
                                deg_range_lon=deg_range/cos(deg2rad*tx_site[z].lat);
                        else
                                deg_range_lon=deg_range/cos(deg2rad*70.0);

                        /* Correct for squares in degrees not being square in miles */  

                        if (deg_range>deg_limit)
                                deg_range=deg_limit;

                        if (deg_range_lon>deg_limit)
                                deg_range_lon=deg_limit;


                        north_min=(int)floor(tx_site[z].lat-deg_range);
                        north_max=(int)floor(tx_site[z].lat+deg_range);

                        west_min=(int)floor(tx_site[z].lon-deg_range_lon);

                        while (west_min<0)
                                west_min+=360;

                        while (west_min>=360)
                                west_min-=360;

                        west_max=(int)floor(tx_site[z].lon+deg_range_lon);

                        while (west_max<0)
                                west_max+=360;

                        while (west_max>=360)
                                west_max-=360;

                        if (north_min<min_lat)
                                min_lat=north_min;

                        if (north_max>max_lat)
                                max_lat=north_max;

                        if (LonDiff(west_min,min_lon)<0.0)
                                min_lon=west_min;

                        if (LonDiff(west_max,max_lon)>0.0)
                                max_lon=west_max;
                }


                /* Load the required SDF files */ 

                width=ReduceAngle(max_lon-min_lon);

                if ((max_lon-min_lon)<=180.0)
                {
                        for (y=0; y<=width; y++)
                                for (x=min_lat; x<=max_lat; x++)
                                {
                                        ymin=(int)(min_lon+(double)y);

                                        while (ymin<0)
                                                ymin+=360;

                                        while (ymin>=360)
                                                ymin-=360;

                                        ymax=ymin+1;

                                        while (ymax<0)
                                                ymax+=360;

                                        while (ymax>=360)
                                                ymax-=360;

                                        sprintf(string,"%d:%d:%d:%d",x, x+1, ymin, ymax);
                                        LoadSDF(string);
                                }
                }

                else
                {
                        for (y=0; y<=width; y++)
                                for (x=min_lat; x<=max_lat; x++)
                                {
                                        ymin=(int)(max_lon+(double)y);

                                        while (ymin<0)
                                                ymin+=360;

                                        while (ymin>=360)
                                                ymin-=360;
                                        
                                        ymax=ymin+1;

                                        while (ymax<0)
                                                ymax+=360;

                                        while (ymax>=360)
                                                ymax-=360;

                                        sprintf(string,"%d:%d:%d:%d",x, x+1, ymin, ymax);
                                        LoadSDF(string);
                                }
                }
        }


        if (mapfile[0])
                map=1;

        if (coverage | LRmap)
        {
                for (x=0; x<txsites; x++)
                {
                        if (coverage)
                                PlotCoverage(tx_site[x],altitude);

                        if (LRmap)
                                PlotLRMap(tx_site[x],altitudeLR);

                        PlaceMarker(tx_site[x]);

                        if (report!='N')
                                SiteReport(tx_site[x]);
                }

                map=1;
        }

        if (coverage==0 && LRmap==0)       
        {
                PlaceMarker(rx_site);

                for (x=0; x<txsites; x++)
                {
                        PlaceMarker(tx_site[x]);

                        if (report=='y')
                        {
                                switch (x)
                                {
                                        case 0:
                                                PlotPath(tx_site[x],rx_site,1);
                                                break;

                                        case 1:
                                                PlotPath(tx_site[x],rx_site,8);
                                                break;

                                        case 2:
                                                PlotPath(tx_site[x],rx_site,16);
                                                break;

                                        case 3:
                                                PlotPath(tx_site[x],rx_site,32);
                                }
                        }

                        if (report!='N')
                                ObstructionReport(tx_site[x],rx_site,report);
                }
        }

        if (map)
        {
                if (bfs)
                {
                        for (x=0; x<bfs; x++)
                                LoadBoundaries(boundary_file[x]);
                }

                if (cities)
                {
                        for (x=0; x<cities; x++)
                                LoadCities(city_file[x]);
                }
                                
                if (!LRmap)
                        WritePPM(mapfile);
                else
                        WritePPMLR(mapfile);
        }

        if (terrain_plot)
        {
                if (txsites>1)
                {
                        for (x=0; terrain_file[x]!='.' && terrain_file[x]!=0 && x<80; x++);

                        if (terrain_file[x]=='.')  /* extension */
                        {
                                ext[0]='.';
                                for (y=1, z=x, x++; terrain_file[x]!=0 && x<253 && y<14; x++, y++)
                                        ext[y]=terrain_file[x];

                                ext[y]=0;
                                terrain_file[z]=0;
                        }

                        else
                        {
                                ext[0]=0;  /* No extension */
                                terrain_file[x]=0;
                        }

                        for (count=0; count<txsites; count++)
                        {
                                sprintf(string,"%s-%c%s%c",terrain_file,'1'+count,ext,0);
                                GraphTerrain(tx_site[count],rx_site,string);
                        }
                }

                else
                        GraphTerrain(tx_site[0],rx_site,terrain_file);
        }

        if (elevation_plot)
        {
                if (txsites>1)
                {
                        for (x=0; elevation_file[x]!='.' && elevation_file[x]!=0 && x<80; x++);

                        if (elevation_file[x]=='.')  /* extension */
                        {
                                ext[0]='.';
                                for (y=1, z=x, x++; elevation_file[x]!=0 && x<253 && y<14; x++, y++)
                                        ext[y]=elevation_file[x];

                                ext[y]=0;
                                elevation_file[z]=0;
                        }

                        else
                        {
                                ext[0]=0;  /* No extension */
                                elevation_file[x]=0;
                        }

                        for (count=0; count<txsites; count++)
                        {
                                sprintf(string,"%s-%c%s%c",elevation_file,'1'+count,ext,0);
                                GraphElevation(tx_site[count],rx_site,string);
                        }
                }

                else
                        GraphElevation(tx_site[0],rx_site,elevation_file);
        }

        if (height_plot)
        {
                if (txsites>1)
                {
                        for (x=0; height_file[x]!='.' && height_file[x]!=0 && x<80; x++);

                        if (height_file[x]=='.')  /* extension */
                        {
                                ext[0]='.';
                                for (y=1, z=x, x++; height_file[x]!=0 && x<253 && y<14; x++, y++)
                                        ext[y]=height_file[x];

                                ext[y]=0;
                                height_file[z]=0;
                        }

                        else
                        {
                                ext[0]=0;  /* No extension */
                                height_file[x]=0;
                        }

                        for (count=0; count<txsites; count++)
                        {
                                sprintf(string,"%s-%c%s%c",height_file,'1'+count,ext,0);
                                GraphHeight(tx_site[count],rx_site,string);
                        }
                }

                else
                        GraphHeight(tx_site[0],rx_site,height_file);
        }

        if (longley_plot)
        {
                if (txsites>1)
                {
                        for (x=0; longley_file[x]!='.' && longley_file[x]!=0 && x<80; x++);

                        if (longley_file[x]=='.')  /* extension */
                        {
                                ext[0]='.';
                                for (y=1, z=x, x++; longley_file[x]!=0 && x<253 && y<14; x++, y++)
                                        ext[y]=longley_file[x];

                                ext[y]=0;
                                longley_file[z]=0;
                        }

                        else
                        {
                                ext[0]=0;  /* No extension */
                                longley_file[x]=0;
                        }

                        for (count=0; count<txsites; count++)
                        {
                                sprintf(string,"%s-%c%s%c",longley_file,'1'+count,ext,0);
                                GraphLongley(tx_site[count],rx_site,string);
                        }
                }

                else
                        GraphLongley(tx_site[0],rx_site,longley_file);
        }

        return 0;
}