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
60 extern double sin (double);
61 extern double cos (double);
62 extern double acos (double);
63 extern double tan (double);
64 extern double atan (double);
65 extern 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
100 char *wwvlat = WWVLAT;
101 char *wwvlong = WWVLONG;
102
103 char *wwvhlat = WWVHLAT;
104 char *wwvhlong = WWVHLONG;
105
106 char *chulat = CHULAT;
107 char *chulong = CHULONG;
108
109 char *goes_up_lat = GOES_UP_LAT;
110 char *goes_up_long = GOES_UP_LONG;
111 char *goes_east_long = GOES_EAST_LONG;
112 char *goes_stby_long = GOES_STBY_LONG;
113 char *goes_west_long = GOES_WEST_LONG;
114 char *goes_sat_lat = GOES_SAT_LAT;
115
116 int hflag = 0;
117 int Wflag = 0;
118 int Cflag = 0;
119 int Gflag = 0;
120 int height;
121
122 char const *progname;
123
124 static void doit (double, double, double, double, double, char *);
125 static double latlong (char *, int);
126 static double greatcircle (double, double, double, double);
127 static double waveangle (double, double, int);
128 static double propdelay (double, double, int);
129 static int finddelay (double, double, double, double, double, double *);
130 static void satdoit (double, double, double, double, double, double, char *);
131 static void satfinddelay (double, double, double, double, double *);
132 static double satpropdelay (double);
133
134 /*
135 * main - parse arguments and handle options
136 */
137 int
main(int argc,char * argv[])138 main(
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 */
276 static void
doit(double lat1,double long1,double lat2,double long2,double h,char * str)277 doit(
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 */
298 static double
latlong(char * str,int islat)299 latlong(
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 */
379 static double
greatcircle(double lat1,double long1,double lat2,double long2)380 greatcircle(
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 */
407 static double
waveangle(double dg,double h,int n)408 waveangle(
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 */
429 static double
propdelay(double dg,double h,int n)430 propdelay(
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 */
453 static int
finddelay(double lat1,double long1,double lat2,double long2,double h,double * delay)454 finddelay(
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 */
488 static void
satdoit(double lat1,double long1,double lat2,double long2,double lat3,double long3,char * str)489 satdoit(
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 */
511 static void
satfinddelay(double lat1,double long1,double lat2,double long2,double * delay)512 satfinddelay(
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 */
531 static double
satpropdelay(double dg)532 satpropdelay(
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