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 volatile 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
main(int argc,char * argv[])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
doit(double lat1,double long1,double lat2,double long2,double h,char * str)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
latlong(char * str,int islat)293 latlong(
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 */
373 static double
greatcircle(double lat1,double long1,double lat2,double long2)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
waveangle(double dg,double h,int n)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
propdelay(double dg,double h,int n)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
finddelay(double lat1,double long1,double lat2,double long2,double h,double * delay)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
satdoit(double lat1,double long1,double lat2,double long2,double lat3,double long3,char * str)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
satfinddelay(double lat1,double long1,double lat2,double long2,double * delay)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
satpropdelay(double dg)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