xref: /freebsd/contrib/ntp/clockstuff/propdelay.c (revision 1e413cf93298b5b97441a21d9a50fdcd0ee9945e)
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 
55 extern	double	sin	(double);
56 extern	double	cos	(double);
57 extern	double	acos	(double);
58 extern	double	tan	(double);
59 extern	double	atan	(double);
60 extern	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 
95 char *wwvlat = WWVLAT;
96 char *wwvlong = WWVLONG;
97 
98 char *wwvhlat = WWVHLAT;
99 char *wwvhlong = WWVHLONG;
100 
101 char *chulat = CHULAT;
102 char *chulong = CHULONG;
103 
104 char *goes_up_lat = GOES_UP_LAT;
105 char *goes_up_long = GOES_UP_LONG;
106 char *goes_east_long = GOES_EAST_LONG;
107 char *goes_stby_long = GOES_STBY_LONG;
108 char *goes_west_long = GOES_WEST_LONG;
109 char *goes_sat_lat = GOES_SAT_LAT;
110 
111 int hflag = 0;
112 int Wflag = 0;
113 int Cflag = 0;
114 int Gflag = 0;
115 int height;
116 
117 char *progname;
118 int debug;
119 
120 static	void	doit		(double, double, double, double, double, char *);
121 static	double	latlong		(char *, int);
122 static	double	greatcircle	(double, double, double, double);
123 static	double	waveangle	(double, double, int);
124 static	double	propdelay	(double, double, int);
125 static	int	finddelay	(double, double, double, double, double, double *);
126 static	void	satdoit		(double, double, double, double, double, double, char *);
127 static	void	satfinddelay	(double, double, double, double, double *);
128 static	double	satpropdelay	(double);
129 
130 /*
131  * main - parse arguments and handle options
132  */
133 int
134 main(
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  */
270 static void
271 doit(
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  */
292 static double
293 latlong(
294 	char *str,
295 	int islat
296 	)
297 {
298 	register char *cp;
299 	register char *bp;
300 	double arg;
301 	double div;
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 		div = 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) / div;
352 			div = 3600.0;
353 		}
354 		if (*cp != '\0')
355 		    arg += atof(cp) / div;
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  */
373 static double
374 greatcircle(
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  */
401 static double
402 waveangle(
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  */
423 static double
424 propdelay(
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  */
447 static int
448 finddelay(
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  */
482 static void
483 satdoit(
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  */
505 static void
506 satfinddelay(
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  */
525 static double
526 satpropdelay(
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 }
545