/*
 * propdelay - compute propagation delays
 *
 * cc -o propdelay propdelay.c -lm
 *
 * "Time and Frequency Users' Manual", NBS Technical Note 695 (1977).
 */

/*
 * This can be used to get a rough idea of the HF propagation delay
 * between two points (usually between you and the radio station).
 * The usage is
 *
 * propdelay latitudeA longitudeA latitudeB longitudeB
 *
 * where points A and B are the locations in question.  You obviously
 * need to know the latitude and longitude of each of the places.
 * The program expects the latitude to be preceded by an 'n' or 's'
 * and the longitude to be preceded by an 'e' or 'w'.  It understands
 * either decimal degrees or degrees:minutes:seconds.  Thus to compute
 * the delay between the WWVH (21:59:26N, 159:46:00W) and WWV (40:40:49N,
 * 105:02:27W) you could use:
 *
 * propdelay n21:59:26 w159:46 n40:40:49 w105:02:27
 *
 * By default it prints out a summer (F2 average virtual height 350 km) and
 * winter (F2 average virtual height 250 km) number.  The results will be
 * quite approximate but are about as good as you can do with HF time anyway.
 * You might pick a number between the values to use, or use the summer
 * value in the summer and switch to the winter value when the static
 * above 10 MHz starts to drop off in the fall.  You can also use the
 * -h switch if you want to specify your own virtual height.
 *
 * You can also do a
 *
 * propdelay -W n45:17:47 w75:45:22
 *
 * to find the propagation delays to WWV and WWVH (from CHU in this
 * case), and a
 *
 * propdelay -C n40:40:49 w105:02:27
 *
 * to find the delays to CHU.
 */

#include <stdio.h>
#include <strings.h>

#define	STREQ(a, b)	(*(a) == *(b) && strcmp((a), (b)) == 0)

/*
 * Program constants
 */
#define	EARTHRADIUS	(6370.0)	/* raduis of earth (km) */
#define	LIGHTSPEED	(299800.0)	/* speed of light, km/s */
#define	PI		(3.1415926536)
#define	RADPERDEG	(PI/180.0)	/* radians per degree */
#define MILE		(1.609344)      /* km in a mile */

#define	SUMMERHEIGHT	(350.0)		/* summer height in km */
#define	WINTERHEIGHT	(250.0)		/* winter height in km */

#define WWVLAT  "n40:40:49"
#define WWVLONG "w105:02:27"

#define WWVHLAT  "n21:59:26"
#define WWVHLONG "w159:46:00"

#define CHULAT	"n45:17:47"
#define	CHULONG	"w75:45:22"

char *wwvlat = WWVLAT;
char *wwvlong = WWVLONG;

char *wwvhlat = WWVHLAT;
char *wwvhlong = WWVHLONG;

char *chulat = CHULAT;
char *chulong = CHULONG;

int hflag = 0;
int Wflag = 0;
int Cflag = 0;
int height;

char *progname;
int debug;

/*
 * main - parse arguments and handle options
 */
main(argc, argv)
int argc;
char *argv[];
{
	int c;
	int errflg = 0;
	double lat1, long1;
	double lat2, long2;
	extern int optind;
	extern char *optarg;
	double latlong();
	void doit();
	extern double atof();

	progname = argv[0];
	while ((c = getopt(argc, argv, "dh:CW")) != EOF)
		switch (c) {
		case 'd':
			++debug;
			break;
		case 'h':
			hflag++;
			height = atof(optarg);
			if (height <= 0.0) {
				(void) fprintf(stderr, "height %s unlikely\n",
				    optarg);
				errflg++;
			}
			break;
		case 'C':
			Cflag++;
			break;
		case 'W':
			Wflag++;
			break;
		default:
			errflg++;
			break;
		}
	if (errflg || (!(Cflag || Wflag) && optind+4 != argc) || 
            ((Cflag || Wflag) && optind+2 != argc)) {
		(void) fprintf(stderr,
		    "usage: %s [-d] [-h height] lat1 long1 lat2 long2\n",
		    progname);
		(void) fprintf(stderr," - or -\n");
		(void) fprintf(stderr,
		    "usage: %s -CW [-d] lat long\n",
		    progname);
		exit(2);
	}

		   
	if (!(Cflag || Wflag)) {
		lat1 = latlong(argv[optind], 1);
		long1 = latlong(argv[optind + 1], 0);
		lat2 = latlong(argv[optind + 2], 1);
		long2 = latlong(argv[optind + 3], 0);
		if (hflag) {
			doit(lat1, long1, lat2, long2, height, "");
		} else {
			doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
			    "summer propagation, ");
			doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
			    "winter propagation, ");
		}
	} else if (Wflag) {
		/*
		 * Compute delay from WWV
	         */
		lat1 = latlong(argv[optind], 1);
		long1 = latlong(argv[optind + 1], 0);
		lat2 = latlong(wwvlat, 1);
		long2 = latlong(wwvlong, 0);
		if (hflag) {
			doit(lat1, long1, lat2, long2, height, "WWV  ");
		} else {
			doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
			    "WWV  summer propagation, ");
			doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
			    "WWV  winter propagation, ");
		}

		/*
		 * Compute delay from WWVH
	         */
		lat2 = latlong(wwvhlat, 1);
		long2 = latlong(wwvhlong, 0);
		if (hflag) {
			doit(lat1, long1, lat2, long2, height, "WWVH ");
		} else {
			doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
			    "WWVH summer propagation, ");
			doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
			    "WWVH winter propagation, ");
		}
	} else if (Cflag) {
		lat1 = latlong(argv[optind], 1);
		long1 = latlong(argv[optind + 1], 0);
		lat2 = latlong(chulat, 1);
		long2 = latlong(chulong, 0);
		if (hflag) {
			doit(lat1, long1, lat2, long2, height, "CHU ");
		} else {
			doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
			    "CHU summer propagation, ");
			doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
			    "CHU winter propagation, ");
		}
	}
	exit(0);
}


/*
 * doit - compute a delay and print it
 */
void
doit(lat1, long1, lat2, long2, h, str)
	double lat1;
	double long1;
	double lat2;
	double long2;
	double h;
	char *str;
{
	int hops;
	double delay;
	int finddelay();

	hops = finddelay(lat1, long1, lat2, long2, h, &delay);
	printf("%sheight %g km, hops %d, delay %g seconds\n",
	    str, h, hops, delay);
}


/*
 * latlong - decode a latitude/longitude value
 */
double
latlong(str, islat)
	char *str;
	int islat;
{
	register char *cp;
	register char *bp;
	double arg;
	double div;
	int isneg;
	char buf[32];
	char *colon;
	extern double atof();

	if (islat) {
		/*
		 * Must be north or south
		 */
		if (*str == 'N' || *str == 'n')
			isneg = 0;
		else if (*str == 'S' || *str == 's')
			isneg = 1;
		else
			isneg = -1;
	} else {
		/*
		 * East is positive, west is negative
		 */
		if (*str == 'E' || *str == 'e')
			isneg = 0;
		else if (*str == 'W' || *str == 'w')
			isneg = 1;
		else
			isneg = -1;
	}

	if (isneg >= 0)
		str++;

	colon = index(str, ':');
	if (colon != NULL) {
		/*
		 * in hhh:mm:ss form
		 */
		cp = str;
		bp = buf;
		while (cp < colon)
			*bp++ = *cp++;
		*bp = '\0';
		cp++;
		arg = atof(buf);
		div = 60.0;
		colon = index(cp, ':');
		if (colon != NULL) {
			bp = buf;
			while (cp < colon)
				*bp++ = *cp++;
			*bp = '\0';
			cp++;
			arg += atof(buf) / div;
			div = 3600.0;
		}
		if (*cp != '\0')
			arg += atof(cp) / div;
	} else {
		arg = atof(str);
	}

	if (isneg == 1)
		arg = -arg;

	if (debug > 2)
		(void) printf("latitude/longitude %s = %g\n", str, arg);

	return arg;
}


/*
 * greatcircle - compute the great circle distance in kilometers
 */
double
greatcircle(lat1, long1, lat2, long2)
	double lat1;
	double long1;
	double lat2;
	double long2;
{
	double dg;
	double l1r, l2r;
	extern double sin();
	extern double cos();
	extern double acos();

	l1r = lat1 * RADPERDEG;
	l2r = lat2 * RADPERDEG;
	dg = EARTHRADIUS * acos(
	    (cos(l1r) * cos(l2r) * cos((long2-long1)*RADPERDEG))
	    + (sin(l1r) * sin(l2r)));
	if (debug >= 2)
		printf(
		    "greatcircle lat1 %g long1 %g lat2 %g long2 %g dist %g\n",
		    lat1, long1, lat2, long2, dg);
	return dg;
}


/*
 * waveangle - compute the wave angle for the given distance, virtual
 *	       height and number of hops.
 */
double
waveangle(dg, h, n)
	double dg;
	double h;
	int n;
{
	double theta;
	double delta;
	extern double tan();
	extern double sin();
	extern double atan();

	theta = dg / (EARTHRADIUS * (double)(2 * n));
	delta = atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2)) - theta;
	if (debug >= 2)
		printf("waveangle dist %g height %g hops %d angle %g\n",
		    dg, h, n, delta / RADPERDEG);
	return delta;
}


/*
 * propdelay - compute the propagation delay
 */
double
propdelay(dg, h, n)
	double dg;
	double h;
	int n;
{
	double phi;
	double theta;
	double td;
	extern double sin();
	extern double tan();
	extern double atan();

	theta = dg / (EARTHRADIUS * (double)(2 * n));
	phi = (PI/2.0) - atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2));
	td = dg / (LIGHTSPEED * sin(phi));
	if (debug >= 2)
		printf("propdelay dist %g height %g hops %d time %g\n",
		    dg, h, n, td);
	return td;
}


/*
 * finddelay - find the propagation delay
 */
int
finddelay(lat1, long1, lat2, long2, h, delay)
	double lat1;
	double long1;
	double lat2;
	double long2;
	double h;
	double *delay;
{
	double dg;	/* great circle distance */
	double delta;	/* wave angle */
	int n;		/* number of hops */

	dg = greatcircle(lat1, long1, lat2, long2);
	if (debug)
		printf("great circle distance %g km %g miles\n", dg, dg/MILE);
	
	n = 1;
	while ((delta = waveangle(dg, h, n)) < 0.0) {
		if (debug)
			printf("tried %d hop%s, no good\n", n, n>1?"s":"");
		n++;
	}
	if (debug)
		printf("%d hop%s okay, wave angle is %g\n", n, n>1?"s":"",
		    delta / RADPERDEG);

	*delay = propdelay(dg, h, n);
	return n;
}