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
58 extern double sin (double);
59 extern double cos (double);
60 extern double acos (double);
61 extern double tan (double);
62 extern double atan (double);
63 extern 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
98 char *wwvlat = WWVLAT;
99 char *wwvlong = WWVLONG;
100
101 char *wwvhlat = WWVHLAT;
102 char *wwvhlong = WWVHLONG;
103
104 char *chulat = CHULAT;
105 char *chulong = CHULONG;
106
107 char *goes_up_lat = GOES_UP_LAT;
108 char *goes_up_long = GOES_UP_LONG;
109 char *goes_east_long = GOES_EAST_LONG;
110 char *goes_stby_long = GOES_STBY_LONG;
111 char *goes_west_long = GOES_WEST_LONG;
112 char *goes_sat_lat = GOES_SAT_LAT;
113
114 int hflag = 0;
115 int Wflag = 0;
116 int Cflag = 0;
117 int Gflag = 0;
118 int height;
119
120 char const *progname;
121
122 static void doit (double, double, double, double, double, char *);
123 static double latlong (char *, int);
124 static double greatcircle (double, double, double, double);
125 static double waveangle (double, double, int);
126 static double propdelay (double, double, int);
127 static int finddelay (double, double, double, double, double, double *);
128 static void satdoit (double, double, double, double, double, double, char *);
129 static void satfinddelay (double, double, double, double, double *);
130 static double satpropdelay (double);
131
132 /*
133 * main - parse arguments and handle options
134 */
135 int
main(int argc,char * argv[])136 main(
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 */
274 static void
doit(double lat1,double long1,double lat2,double long2,double h,char * str)275 doit(
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 */
296 static double
latlong(char * str,int islat)297 latlong(
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 */
377 static double
greatcircle(double lat1,double long1,double lat2,double long2)378 greatcircle(
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 */
405 static double
waveangle(double dg,double h,int n)406 waveangle(
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 */
427 static double
propdelay(double dg,double h,int n)428 propdelay(
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 */
451 static int
finddelay(double lat1,double long1,double lat2,double long2,double h,double * delay)452 finddelay(
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 */
486 static void
satdoit(double lat1,double long1,double lat2,double long2,double lat3,double long3,char * str)487 satdoit(
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 */
509 static void
satfinddelay(double lat1,double long1,double lat2,double long2,double * delay)510 satfinddelay(
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 */
529 static double
satpropdelay(double dg)530 satpropdelay(
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