1/*	$NetBSD: propdelay.c,v 1.6 2020/05/25 20:47:19 christos Exp $	*/
2
3/* propdelay.c,v 3.1 1993/07/06 01:05:24 jbj Exp
4 * propdelay - compute propagation delays
5 *
6 * cc -o propdelay propdelay.c -lm
7 *
8 * "Time and Frequency Users' Manual", NBS Technical Note 695 (1977).
9 */
10
11/*
12 * This can be used to get a rough idea of the HF propagation delay
13 * between two points (usually between you and the radio station).
14 * The usage is
15 *
16 * propdelay latitudeA longitudeA latitudeB longitudeB
17 *
18 * where points A and B are the locations in question.  You obviously
19 * need to know the latitude and longitude of each of the places.
20 * The program expects the latitude to be preceded by an 'n' or 's'
21 * and the longitude to be preceded by an 'e' or 'w'.  It understands
22 * either decimal degrees or degrees:minutes:seconds.  Thus to compute
23 * the delay between the WWVH (21:59:26N, 159:46:00W) and WWV (40:40:49N,
24 * 105:02:27W) you could use:
25 *
26 * propdelay n21:59:26 w159:46 n40:40:49 w105:02:27
27 *
28 * By default it prints out a summer (F2 average virtual height 350 km) and
29 * winter (F2 average virtual height 250 km) number.  The results will be
30 * quite approximate but are about as good as you can do with HF time anyway.
31 * You might pick a number between the values to use, or use the summer
32 * value in the summer and switch to the winter value when the static
33 * above 10 MHz starts to drop off in the fall.  You can also use the
34 * -h switch if you want to specify your own virtual height.
35 *
36 * You can also do a
37 *
38 * propdelay -W n45:17:47 w75:45:22
39 *
40 * to find the propagation delays to WWV and WWVH (from CHU in this
41 * case), a
42 *
43 * propdelay -C n40:40:49 w105:02:27
44 *
45 * to find the delays to CHU, and a
46 *
47 * propdelay -G n52:03:17 w98:34:18
48 *
49 * to find delays to GOES via each of the three satellites.
50 */
51
52#ifdef HAVE_CONFIG_H
53# include <config.h>
54#endif
55#include <stdio.h>
56#include <string.h>
57
58#include "ntp_stdlib.h"
59
60extern	double	sin	(double);
61extern	double	cos	(double);
62extern	double	acos	(double);
63extern	double	tan	(double);
64extern	double	atan	(double);
65extern	double	sqrt	(double);
66
67#define	STREQ(a, b)	(*(a) == *(b) && strcmp((a), (b)) == 0)
68
69/*
70 * Program constants
71 */
72#define	EARTHRADIUS	(6370.0)	/* raduis of earth (km) */
73#define	LIGHTSPEED	(299800.0)	/* speed of light, km/s */
74#define	PI		(3.1415926536)
75#define	RADPERDEG	(PI/180.0)	/* radians per degree */
76#define MILE		(1.609344)      /* km in a mile */
77
78#define	SUMMERHEIGHT	(350.0)		/* summer height in km */
79#define	WINTERHEIGHT	(250.0)		/* winter height in km */
80
81#define SATHEIGHT	(6.6110 * 6378.0) /* geosync satellite height in km
82					     from centre of earth */
83
84#define WWVLAT  "n40:40:49"
85#define WWVLONG "w105:02:27"
86
87#define WWVHLAT  "n21:59:26"
88#define WWVHLONG "w159:46:00"
89
90#define CHULAT	"n45:17:47"
91#define	CHULONG	"w75:45:22"
92
93#define GOES_UP_LAT  "n37:52:00"
94#define GOES_UP_LONG "w75:27:00"
95#define GOES_EAST_LONG "w75:00:00"
96#define GOES_STBY_LONG "w105:00:00"
97#define GOES_WEST_LONG "w135:00:00"
98#define GOES_SAT_LAT "n00:00:00"
99
100char *wwvlat = WWVLAT;
101char *wwvlong = WWVLONG;
102
103char *wwvhlat = WWVHLAT;
104char *wwvhlong = WWVHLONG;
105
106char *chulat = CHULAT;
107char *chulong = CHULONG;
108
109char *goes_up_lat = GOES_UP_LAT;
110char *goes_up_long = GOES_UP_LONG;
111char *goes_east_long = GOES_EAST_LONG;
112char *goes_stby_long = GOES_STBY_LONG;
113char *goes_west_long = GOES_WEST_LONG;
114char *goes_sat_lat = GOES_SAT_LAT;
115
116int hflag = 0;
117int Wflag = 0;
118int Cflag = 0;
119int Gflag = 0;
120int height;
121
122char const *progname;
123
124static	void	doit		(double, double, double, double, double, char *);
125static	double	latlong		(char *, int);
126static	double	greatcircle	(double, double, double, double);
127static	double	waveangle	(double, double, int);
128static	double	propdelay	(double, double, int);
129static	int	finddelay	(double, double, double, double, double, double *);
130static	void	satdoit		(double, double, double, double, double, double, char *);
131static	void	satfinddelay	(double, double, double, double, double *);
132static	double	satpropdelay	(double);
133
134/*
135 * main - parse arguments and handle options
136 */
137int
138main(
139	int argc,
140	char *argv[]
141	)
142{
143	int c;
144	int errflg = 0;
145	double lat1, long1;
146	double lat2, long2;
147	double lat3, long3;
148
149	init_lib();
150
151	progname = argv[0];
152	while ((c = ntp_getopt(argc, argv, "dh:CWG")) != EOF)
153	    switch (c) {
154		case 'd':
155		    ++debug;
156		    break;
157		case 'h':
158		    hflag++;
159		    height = atof(ntp_optarg);
160		    if (height <= 0.0) {
161			    (void) fprintf(stderr, "height %s unlikely\n",
162					   ntp_optarg);
163			    errflg++;
164		    }
165		    break;
166		case 'C':
167		    Cflag++;
168		    break;
169		case 'W':
170		    Wflag++;
171		    break;
172		case 'G':
173		    Gflag++;
174		    break;
175		default:
176		    errflg++;
177		    break;
178	    }
179	if (errflg || (!(Cflag || Wflag || Gflag) && ntp_optind+4 != argc) ||
180	    ((Cflag || Wflag || Gflag) && ntp_optind+2 != argc)) {
181		(void) fprintf(stderr,
182			       "usage: %s [-d] [-h height] lat1 long1 lat2 long2\n",
183			       progname);
184		(void) fprintf(stderr," - or -\n");
185		(void) fprintf(stderr,
186			       "usage: %s -CWG [-d] lat long\n",
187			       progname);
188		exit(2);
189	}
190
191
192	if (!(Cflag || Wflag || Gflag)) {
193		lat1 = latlong(argv[ntp_optind], 1);
194		long1 = latlong(argv[ntp_optind + 1], 0);
195		lat2 = latlong(argv[ntp_optind + 2], 1);
196		long2 = latlong(argv[ntp_optind + 3], 0);
197		if (hflag) {
198			doit(lat1, long1, lat2, long2, height, "");
199		} else {
200			doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
201			     "summer propagation, ");
202			doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
203			     "winter propagation, ");
204		}
205	} else if (Wflag) {
206		/*
207		 * Compute delay from WWV
208		 */
209		lat1 = latlong(argv[ntp_optind], 1);
210		long1 = latlong(argv[ntp_optind + 1], 0);
211		lat2 = latlong(wwvlat, 1);
212		long2 = latlong(wwvlong, 0);
213		if (hflag) {
214			doit(lat1, long1, lat2, long2, height, "WWV  ");
215		} else {
216			doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
217			     "WWV  summer propagation, ");
218			doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
219			     "WWV  winter propagation, ");
220		}
221
222		/*
223		 * Compute delay from WWVH
224		 */
225		lat2 = latlong(wwvhlat, 1);
226		long2 = latlong(wwvhlong, 0);
227		if (hflag) {
228			doit(lat1, long1, lat2, long2, height, "WWVH ");
229		} else {
230			doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
231			     "WWVH summer propagation, ");
232			doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
233			     "WWVH winter propagation, ");
234		}
235	} else if (Cflag) {
236		lat1 = latlong(argv[ntp_optind], 1);
237		long1 = latlong(argv[ntp_optind + 1], 0);
238		lat2 = latlong(chulat, 1);
239		long2 = latlong(chulong, 0);
240		if (hflag) {
241			doit(lat1, long1, lat2, long2, height, "CHU ");
242		} else {
243			doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
244			     "CHU summer propagation, ");
245			doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
246			     "CHU winter propagation, ");
247		}
248	} else if (Gflag) {
249		lat1 = latlong(goes_up_lat, 1);
250		long1 = latlong(goes_up_long, 0);
251		lat3 = latlong(argv[ntp_optind], 1);
252		long3 = latlong(argv[ntp_optind + 1], 0);
253
254		lat2 = latlong(goes_sat_lat, 1);
255
256		long2 = latlong(goes_west_long, 0);
257		satdoit(lat1, long1, lat2, long2, lat3, long3,
258			"GOES Delay via WEST");
259
260		long2 = latlong(goes_stby_long, 0);
261		satdoit(lat1, long1, lat2, long2, lat3, long3,
262			"GOES Delay via STBY");
263
264		long2 = latlong(goes_east_long, 0);
265		satdoit(lat1, long1, lat2, long2, lat3, long3,
266			"GOES Delay via EAST");
267
268	}
269	exit(0);
270}
271
272
273/*
274 * doit - compute a delay and print it
275 */
276static void
277doit(
278	double lat1,
279	double long1,
280	double lat2,
281	double long2,
282	double h,
283	char *str
284	)
285{
286	int hops;
287	double delay;
288
289	hops = finddelay(lat1, long1, lat2, long2, h, &delay);
290	printf("%sheight %g km, hops %d, delay %g seconds\n",
291	       str, h, hops, delay);
292}
293
294
295/*
296 * latlong - decode a latitude/longitude value
297 */
298static double
299latlong(
300	char *str,
301	int islat
302	)
303{
304	register char *cp;
305	register char *bp;
306	double arg;
307	double divby;
308	int isneg;
309	char buf[32];
310	char *colon;
311
312	if (islat) {
313		/*
314		 * Must be north or south
315		 */
316		if (*str == 'N' || *str == 'n')
317		    isneg = 0;
318		else if (*str == 'S' || *str == 's')
319		    isneg = 1;
320		else
321		    isneg = -1;
322	} else {
323		/*
324		 * East is positive, west is negative
325		 */
326		if (*str == 'E' || *str == 'e')
327		    isneg = 0;
328		else if (*str == 'W' || *str == 'w')
329		    isneg = 1;
330		else
331		    isneg = -1;
332	}
333
334	if (isneg >= 0)
335	    str++;
336
337	colon = strchr(str, ':');
338	if (colon != NULL) {
339		/*
340		 * in hhh:mm:ss form
341		 */
342		cp = str;
343		bp = buf;
344		while (cp < colon)
345		    *bp++ = *cp++;
346		*bp = '\0';
347		cp++;
348		arg = atof(buf);
349		divby = 60.0;
350		colon = strchr(cp, ':');
351		if (colon != NULL) {
352			bp = buf;
353			while (cp < colon)
354			    *bp++ = *cp++;
355			*bp = '\0';
356			cp++;
357			arg += atof(buf) / divby;
358			divby = 3600.0;
359		}
360		if (*cp != '\0')
361		    arg += atof(cp) / divby;
362	} else {
363		arg = atof(str);
364	}
365
366	if (isneg == 1)
367	    arg = -arg;
368
369	if (debug > 2)
370	    (void) printf("latitude/longitude %s = %g\n", str, arg);
371
372	return arg;
373}
374
375
376/*
377 * greatcircle - compute the great circle distance in kilometers
378 */
379static double
380greatcircle(
381	double lat1,
382	double long1,
383	double lat2,
384	double long2
385	)
386{
387	double dg;
388	double l1r, l2r;
389
390	l1r = lat1 * RADPERDEG;
391	l2r = lat2 * RADPERDEG;
392	dg = EARTHRADIUS * acos(
393		(cos(l1r) * cos(l2r) * cos((long2-long1)*RADPERDEG))
394		+ (sin(l1r) * sin(l2r)));
395	if (debug >= 2)
396	    printf(
397		    "greatcircle lat1 %g long1 %g lat2 %g long2 %g dist %g\n",
398		    lat1, long1, lat2, long2, dg);
399	return dg;
400}
401
402
403/*
404 * waveangle - compute the wave angle for the given distance, virtual
405 *	       height and number of hops.
406 */
407static double
408waveangle(
409	double dg,
410	double h,
411	int n
412	)
413{
414	double theta;
415	double delta;
416
417	theta = dg / (EARTHRADIUS * (double)(2 * n));
418	delta = atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2)) - theta;
419	if (debug >= 2)
420	    printf("waveangle dist %g height %g hops %d angle %g\n",
421		   dg, h, n, delta / RADPERDEG);
422	return delta;
423}
424
425
426/*
427 * propdelay - compute the propagation delay
428 */
429static double
430propdelay(
431	double dg,
432	double h,
433	int n
434	)
435{
436	double phi;
437	double theta;
438	double td;
439
440	theta = dg / (EARTHRADIUS * (double)(2 * n));
441	phi = (PI/2.0) - atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2));
442	td = dg / (LIGHTSPEED * sin(phi));
443	if (debug >= 2)
444	    printf("propdelay dist %g height %g hops %d time %g\n",
445		   dg, h, n, td);
446	return td;
447}
448
449
450/*
451 * finddelay - find the propagation delay
452 */
453static int
454finddelay(
455	double lat1,
456	double long1,
457	double lat2,
458	double long2,
459	double h,
460	double *delay
461	)
462{
463	double dg;	/* great circle distance */
464	double delta;	/* wave angle */
465	int n;		/* number of hops */
466
467	dg = greatcircle(lat1, long1, lat2, long2);
468	if (debug)
469	    printf("great circle distance %g km %g miles\n", dg, dg/MILE);
470
471	n = 1;
472	while ((delta = waveangle(dg, h, n)) < 0.0) {
473		if (debug)
474		    printf("tried %d hop%s, no good\n", n, n>1?"s":"");
475		n++;
476	}
477	if (debug)
478	    printf("%d hop%s okay, wave angle is %g\n", n, n>1?"s":"",
479		   delta / RADPERDEG);
480
481	*delay = propdelay(dg, h, n);
482	return n;
483}
484
485/*
486 * satdoit - compute a delay and print it
487 */
488static void
489satdoit(
490	double lat1,
491	double long1,
492	double lat2,
493	double long2,
494	double lat3,
495	double long3,
496	char *str
497	)
498{
499	double up_delay,down_delay;
500
501	satfinddelay(lat1, long1, lat2, long2, &up_delay);
502	satfinddelay(lat3, long3, lat2, long2, &down_delay);
503
504	printf("%s, delay %g seconds\n", str, up_delay + down_delay);
505}
506
507/*
508 * satfinddelay - calculate the one-way delay time between a ground station
509 * and a satellite
510 */
511static void
512satfinddelay(
513	double lat1,
514	double long1,
515	double lat2,
516	double long2,
517	double *delay
518	)
519{
520	double dg;	/* great circle distance */
521
522	dg = greatcircle(lat1, long1, lat2, long2);
523
524	*delay = satpropdelay(dg);
525}
526
527/*
528 * satpropdelay - calculate the one-way delay time between a ground station
529 * and a satellite
530 */
531static double
532satpropdelay(
533	double dg
534	)
535{
536	double k1, k2, dist;
537	double theta;
538	double td;
539
540	theta = dg / (EARTHRADIUS);
541	k1 = EARTHRADIUS * sin(theta);
542	k2 = SATHEIGHT - (EARTHRADIUS * cos(theta));
543	if (debug >= 2)
544	    printf("Theta %g k1 %g k2 %g\n", theta, k1, k2);
545	dist = sqrt(k1*k1 + k2*k2);
546	td = dist / LIGHTSPEED;
547	if (debug >= 2)
548	    printf("propdelay dist %g height %g time %g\n", dg, dist, td);
549	return td;
550}
551