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