/*
 *      gc.c
 *
 *      Great Circle.  This program is used to determine bearing
 *      and range to two stations given latitude and longitude.
 *
 *      Ver 1.07 By S. R. Sampson, N5OWK
 *      Public Domain (p) June 1993
 *
 *      Ref: "Air Navigation", Air Force Manual 51-40, 1 February 1987
 *      Ref: "ARRL Satellite Experimenters Handbook", August 1990
 *
 *      Usage examples:
 *
 *      gc n 35.19n97.27w 0s0e           (Moore to Prime/Equator)
 *      gc n 35.19N97.27W 38.51n77.02W   (Moore to Washington D.C., mixed case)
 *      gc n 33.56n118.24w 55.45n37.35e  (L.A. to Moscow)
 *	gc n 35N70W 35N71W		 (No decimal points used, all uppercase)
 *
 *	Modified the program to incorporate short and long path information
 *	from the Satellite Handbook.  This version also takes into consideration
 *	the two points being close enough to be in the near-field, and the
 *	antipodal points, which are easily calculated.  These last points were
 *	made in discussions with John Allison who makes the nice MAPIT program.
 *
 *	Compile GNU C with: cc -O gc.c -o gc -lm
 */

/* Includes */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>

/* Defines */

#define RADIAN  (180.0 / M_PI)

/* Globals */

struct	{
	double	miles;	/* arc length for 1 degree, various units of measure */
	char	*text;
} Units[] = {
	{ 60.0, "Nautical Miles"},
	{ 111.2, "kilometers"},
	{ 69.1, "Statute Miles"}
};

/* Simple Declare, No Prototypes */

/*
 *      Error routine
 */

void err(type)
int     type;
{
        switch(type)  {
        case 1:
                printf("\007Latitude Out of Range (90N to 90S)\n");
                break;
        case 2:
                printf("\007Longitude Out of Range (180W to 180E)\n");
                break;
        case 3:
                printf("\007Minutes Out of Range (0 to 59)\n");
        }

        exit(-1);
}

/*
 *      Convert Degrees and Minutes to Decimal
 */

double dm2dec(n)
double  n;
{
        double  t;

        t = (int)n;
        n -= t;
        n /= .60;

        if (n >= 1.0)
                err(3);

        return (n + t);
}

/*
 *      Parse the input line
 *
 *      dd(.mm)[NnSs]ddd(.mm)[EeWw]
 */

void parse(s, lat, lon)
char    *s;
double  *lat, *lon;
{
        register char   *i, *t, *e;

        e = s + strlen(s);
        for (i = s; i < e; ++i)  {
                switch (*i)  {
                case 'n':
                case 'N':
                        *i = '\0';
                        t = i + 1;
                        *lat = atof(s);
                        break;
                case 's':
                case 'S':
                        *i = '\0';
                        t = i + 1;
                        *lat = -atof(s);
                        break;
                case 'e':
                case 'E':
                        *i = '\0';
                        *lon = -atof(t);
                        break;
                case 'w':
                case 'W':
                        *i = '\0';
                        *lon = atof(t);
                }
        }

        *lat = dm2dec(*lat);
        *lon = dm2dec(*lon);

        if (*lat > 90.0 || *lat < -90.0)
                err(1);

        if (*lon > 180.0 || *lon < -180.0)
                err(2);

        /* Prevent ACOS() Domain Error */

        if (*lat == 90.0)
                *lat = 89.9;

        if (*lat == -90.0)
                *lat = -89.9;
}

void main(argc, argv)
int  argc;
char **argv;
{
	double	tmp, arc, cosaz, az, azsp, azlp, distsp, distlp;
	double  QTH_Lat, QTH_Long, DEST_Lat, DEST_Long, Delta_Long;
	int	units;

        if (argc != 4)  {
            fprintf(stderr, "\nUsage: gc units station1 station2\n\n" \
            "This program computes Great Circle Bearing and Range\n" \
            "given the latitude and longitude (degrees and minutes).\n\n" \
            "You must input the lat/long of the two stations.\n" \
            "The output will then be relative from station1 to station2.\n\n" \
            "Input the two station lat/longs using the following format:\n\n" \
            "\tdd.mmHddd.mmG  lead/lagging zeros can be left out.\n\n" \
            "d = Degrees, m = Minutes, H = Hemisphere (N or S), " \
	    "G = Greenwich (W or E)\n\n" \
            "units is 'n' for Nautical, 'k' for kilometers, and 's' for " \
            "Statute.\n\n");

                exit(1);
        }

        /* Process the command line data */

	switch (argv[1][0])  {
	case 'k':
	case 'K':
		units = 1;
		break;
	case 's':
	case 'S':
		units = 2;
		break;
	case 'n':
	case 'N':
	default:
		units = 0;
	}

        parse(argv[2], &QTH_Lat, &QTH_Long);
        parse(argv[3], &DEST_Lat, &DEST_Long);

        QTH_Lat    /= RADIAN;   /* Convert variables to Radians */
        QTH_Long   /= RADIAN;
        DEST_Lat   /= RADIAN;
        DEST_Long  /= RADIAN;

        Delta_Long = DEST_Long - QTH_Long;

        tmp = (sin(QTH_Lat) * sin(DEST_Lat)) +
                (cos(QTH_Lat) * cos(DEST_Lat) * cos(Delta_Long));

	if (tmp > .999999)  {
		printf("Station points coincide, use an Omni!\n\n");
		exit(0);
	} else if (tmp < -.999999)  {
		/*
		 * points are antipodal, he's straight down.
		 * So take 180 Degrees of arc times 60 nm,
		 * and you get 10800 nm, or whatever units...
		 */

		printf("Station is equal distance in all Azimuths " \
		       "(antipodal)\n%.0f %s\n\n",
		       (Units[units].miles) * 180.0,
		       Units[units].text);

		exit(0);
	} else  {
		arc = acos(tmp);

		/*
		 * One degree of arc is 60 Nautical miles
		 * at the surface of the earth, 111.2 km, or 69.1 sm
		 * This method is easier than the one in the handbook
		 */

		/* Short Path */

		distsp = (Units[units].miles) * (arc * RADIAN);

		/* Long Path */

		distlp = ((Units[units].miles) * 360.0) - distsp;

	}

	cosaz = (sin(DEST_Lat) - (sin(QTH_Lat) * cos(arc))) /
                (sin(arc) * cos(QTH_Lat));

	if (cosaz > .999999)
		az = 0.0;
	else if (cosaz < -.999999)
		az = 180.0;
	else
		az = acos(cosaz) * RADIAN;

	/*
	 * Handbook had the test ">= 0.0" which looks backwards??
	 */

        if (sin(Delta_Long) < 0.0)  {
		azsp = az;
		azlp = 180.0 + az;
	} else  {
		azsp = 360.0 - az;
		azlp = 180.0 - az;
	}

        /* Computations complete, show answer */

        printf("Short Path Bearing is %03.0f Degrees for %.0f %s\n",
	        azsp, distsp, Units[units].text);

        printf(" Long Path Bearing is %03.0f Degrees for %.0f %s\n",
	        azlp, distlp, Units[units].text);

        exit(0);
}