+/***************************************************************************\
+* bearing: Determines distance and azimuth bearing between locations *
+* specified in a pair of .qth files *
+* Last update: 08-Feb-2009 *
+*****************************************************************************
+* Project started on December 7, 2007 by John A. Magliacane, KD2BD *
+*****************************************************************************
+* *
+* 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. *
+* *
+*****************************************************************************
+* gcc -Wall -O3 -s -lm -fomit-frame-pointer bearing.c -o bearing *
+\***************************************************************************/
+
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include <ctype.h>
+#include <unistd.h>
+
+char string[255];
+
+double TWOPI=6.283185307179586, PI=3.141592653589793,
+ deg2rad=1.74532925199e-02, KM_PER_MILE=1.609344;
+
+struct site { double lat;
+ double lon;
+ double azimuth;
+ char name[50];
+ } site;
+
+
+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;
+}
+
+char *dec2dms(double decimal)
+{
+ /* Converts decimal degrees to degrees, minutes, seconds,
+ (DMS) and returns the result as a character string. */
+
+ char sign;
+ int degrees, minutes, seconds;
+ double a, b, c, d;
+
+ if (decimal<0.0)
+ {
+ decimal=-decimal;
+ sign=-1;
+ }
+
+ else
+ sign=1;
+
+ 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*sign, 176, minutes, seconds);
+ return (string);
+}
+
+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 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 && input[a+1]!='\n' && 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.xx) */
+ {
+ sscanf(string,"%d %d %lf",°rees, &minutes, &seconds);
+
+ bearing=fabs((double)degrees);
+ bearing+=fabs(((double)minutes)/60.0);
+ bearing+=fabs(seconds/3600.0);
+
+ if ((degrees<0) || (minutes<0) || (seconds<0.0))
+ bearing=-bearing;
+ }
+
+ /* 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], *s=NULL;
+ 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.name[0]=0;
+ tempsite.azimuth=0.0;
+
+ fd=fopen(qthfile,"r");
+
+ if (fd!=NULL)
+ {
+ /* Site Name */
+ s=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 */
+ s=fgets(string,49,fd);
+ tempsite.lat=ReadBearing(string);
+
+ /* Site Longitude */
+ s=fgets(string,49,fd);
+ tempsite.lon=ReadBearing(string);
+
+ fclose(fd);
+ }
+
+ return tempsite;
+}
+
+int main(int argc, char *argv[])
+{
+ int x, y;
+ unsigned char sitenum, metric, error;
+ double distance, azimuth;
+ struct site point[2];
+
+ if (argc==1)
+ {
+ fprintf(stdout,"\nProvides distance and azimuth bearing between two site locations.\n");
+ fprintf(stdout,"\nUsage: bearing site1(.qth) site2(.qth)");
+ fprintf(stdout,"\n bearing site1(.qth) site2(.qth) -metric\n\n");
+ fflush(stdout);
+ return 1;
+ }
+
+ metric=0;
+ error=0;
+ y=argc-1;
+ sitenum=0;
+
+ point[0].lat=91.0;
+ point[0].lon=361.0;
+ point[1].lat=91.0;
+ point[1].lon=361.0;
+
+ /* Scan for command line arguments */
+
+ for (x=1; x<=y; x++)
+ {
+ if (strcmp(argv[x],"-metric")==0)
+ {
+ metric=1;
+ x++;
+ }
+
+ /* Read Receiver Location */
+
+ if (x<=y && argv[x][0] && argv[x][0]!='-' && sitenum<2)
+ {
+ point[sitenum]=LoadQTH(argv[x]);
+
+ if (point[sitenum].lat>90.0 || point[sitenum].lon>360.0)
+ fprintf(stderr,"\n%c*** \"%s\" is not a valid .qth file!",7,argv[x]);
+ else
+ sitenum++;
+ }
+ }
+
+ if (sitenum<2)
+ {
+ fprintf(stderr,"\n%c*** ERROR: Not enough valid sites specified!\n\n",7);
+ exit (-1);
+ }
+
+ for (x=0; x<sitenum; x++)
+ {
+ if (point[x].lat>90.0 && point[x].lon>360.0)
+ {
+ fprintf(stderr,"\n*** ERROR: site #%d not found!",x+1);
+ error=1;
+ }
+ }
+
+ if (error)
+ return -1;
+
+ else
+ {
+ /* Perform the calculations and display the results */
+
+ distance=Distance(point[0],point[1]);
+ azimuth=Azimuth(point[0],point[1]);
+
+ printf("\nThe distance between %s and %s is\n",point[0].name, point[1].name);
+
+ if (metric)
+ printf("%.2f kilometers",distance*KM_PER_MILE);
+ else
+ printf("%.2f miles",distance);
+
+ printf(" at a bearing of %.2f%c azimuth.\n\n",azimuth, 176);
+ }
+
+ return 0;
+}