/contrib/ntp/clockstuff/propdelay.c

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