/*
 *      gc.c
 *
 *      Great Circle.  This program is used to determine bearing
 *      and range to a station given latitude and longitude.
 *
 *      Ver 1.03 By S. R. Sampson, N5OWK
 *      Public Domain (p) November 1989
 *
 *      Ref: Air Force Manual 51-40, "Air Navigation", 1 February 1987
 *
 *      Usage examples:
 *
 *      gc 35.19n97.27w 0s0e            (Moore to Prime/Equator)
 *      gc 35.19N97.27W 38.51n77.02W    (Moore to Washington D.C.)
 *      gc 33.56n118.24w 55.45n37.35e   (L.A. to Moscow U.S.S.R.)
 *	gc 35N70W 35N71W		(No decimal points used)
 */

/* Includes */

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

/* Defines */

#define RADIAN  (180.0 / M_PI)

/* Globals */

double  tmp,
        dist,
        range,
        bearing,
        QTH_Lat,
        QTH_Long,
        DEST_Lat,
        DEST_Long,
        Delta_Long;

/* 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;
        int             l;

        l = strlen(s);
        for (i = s; i < (s + l); ++i)  {
                switch(toupper(*i))  {
                case 'N':
                        *i = '\0';
                        t = i + 1;
                        *lat = atof(s);
                        break;
                case 'S':
                        *i = '\0';
                        t = i + 1;
                        *lat = -atof(s);
                        break;
                case 'E':
                        *i = '\0';
                        *lon = -atof(t);
                        break;
                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;
}

main(argc, argv)
register int  argc;
register char **argv;
{
        if (argc != 3)  {
                printf("Usage: gc station1 station2\n\n");
                printf("This program computes Great Circle Bearing and Range\n");
                printf("given the latitude and longitude (degrees and minutes).\n\n");
                printf("You must input the lat/long of the two stations.\n");
                printf("The output will then be relative from station1 to station2.\n\n");
                printf("Input the two station lat/longs using the following format:\n\n");
                printf("\tdd.mmHddd.mmG  lead/lagging zeros can be left out.\n\n");
                printf("d = Degrees, m = Minutes, H = Hemisphere (N or S), G = Greenwich (W or E)\n");

                exit(1);
        }

        /* Process the command line data */

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

        /* Compute the Bearing and Range, From the Formula in Chapter 23 */

        Delta_Long = DEST_Long - QTH_Long;

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

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

        dist = acos(tmp);
        range = 60.0 * (dist * RADIAN);

        tmp = (sin(DEST_Lat) - (sin(QTH_Lat) * cos(dist))) /
                (sin(dist) * cos(QTH_Lat));

        bearing = acos(tmp) * RADIAN;

        if (Delta_Long > 0.0)
                bearing = 360.0 - bearing;

        /* Computations complete, show answer */

        printf("\nBearing is %.0f Degrees for %.0f Nautical Miles\n",
                bearing, range);

        exit(0);
}

/* EOF */