![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to propdelay.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [ELWIX - Embedded LightWeight unIX -] / embedaddon / ntp / clockstuff|
Return to propdelay.c CVS log | Up to [ELWIX - Embedded LightWeight unIX -] / embedaddon / ntp / clockstuff |
1.1 misho 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
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 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
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: }