1 /* propdelay.c,v 3.1 1993/07/06 01:05:24 jbj Exp
2 * propdelay - compute propagation delays
4 * cc -o propdelay propdelay.c -lm
6 * "Time and Frequency Users' Manual", NBS Technical Note 695 (1977).
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).
14 * propdelay latitudeA longitudeA latitudeB longitudeB
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:
24 * propdelay n21:59:26 w159:46 n40:40:49 w105:02:27
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.
36 * propdelay -W n45:17:47 w75:45:22
38 * to find the propagation delays to WWV and WWVH (from CHU in this
41 * propdelay -C n40:40:49 w105:02:27
43 * to find the delays to CHU, and a
45 * propdelay -G n52:03:17 w98:34:18
47 * to find delays to GOES via each of the three satellites.
53 #include "ntp_stdlib.h"
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);
62 #define STREQ(a, b) (*(a) == *(b) && strcmp((a), (b)) == 0)
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 */
73 #define SUMMERHEIGHT (350.0) /* summer height in km */
74 #define WINTERHEIGHT (250.0) /* winter height in km */
76 #define SATHEIGHT (6.6110 * 6378.0) /* geosync satellite height in km
77 from centre of earth */
79 #define WWVLAT "n40:40:49"
80 #define WWVLONG "w105:02:27"
82 #define WWVHLAT "n21:59:26"
83 #define WWVHLONG "w159:46:00"
85 #define CHULAT "n45:17:47"
86 #define CHULONG "w75:45:22"
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"
95 char *wwvlat = WWVLAT;
96 char *wwvlong = WWVLONG;
98 char *wwvhlat = WWVHLAT;
99 char *wwvhlong = WWVHLONG;
101 char *chulat = CHULAT;
102 char *chulong = CHULONG;
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;
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);
131 * main - parse arguments and handle options
146 while ((c = ntp_getopt(argc, argv, "dh:CWG")) != EOF)
153 height = atof(ntp_optarg);
155 (void) fprintf(stderr, "height %s unlikely\n",
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",
178 (void) fprintf(stderr," - or -\n");
179 (void) fprintf(stderr,
180 "usage: %s -CWG [-d] lat long\n",
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);
192 doit(lat1, long1, lat2, long2, height, "");
194 doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
195 "summer propagation, ");
196 doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
197 "winter propagation, ");
201 * Compute delay from WWV
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);
208 doit(lat1, long1, lat2, long2, height, "WWV ");
210 doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
211 "WWV summer propagation, ");
212 doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
213 "WWV winter propagation, ");
217 * Compute delay from WWVH
219 lat2 = latlong(wwvhlat, 1);
220 long2 = latlong(wwvhlong, 0);
222 doit(lat1, long1, lat2, long2, height, "WWVH ");
224 doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
225 "WWVH summer propagation, ");
226 doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
227 "WWVH winter propagation, ");
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);
235 doit(lat1, long1, lat2, long2, height, "CHU ");
237 doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
238 "CHU summer propagation, ");
239 doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
240 "CHU winter propagation, ");
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);
248 lat2 = latlong(goes_sat_lat, 1);
250 long2 = latlong(goes_west_long, 0);
251 satdoit(lat1, long1, lat2, long2, lat3, long3,
252 "GOES Delay via WEST");
254 long2 = latlong(goes_stby_long, 0);
255 satdoit(lat1, long1, lat2, long2, lat3, long3,
256 "GOES Delay via STBY");
258 long2 = latlong(goes_east_long, 0);
259 satdoit(lat1, long1, lat2, long2, lat3, long3,
260 "GOES Delay via EAST");
268 * doit - compute a delay and print it
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);
290 * latlong - decode a latitude/longitude value
308 * Must be north or south
310 if (*str == 'N' || *str == 'n')
312 else if (*str == 'S' || *str == 's')
318 * East is positive, west is negative
320 if (*str == 'E' || *str == 'e')
322 else if (*str == 'W' || *str == 'w')
331 colon = strchr(str, ':');
344 colon = strchr(cp, ':');
351 arg += atof(buf) / div;
355 arg += atof(cp) / div;
364 (void) printf("latitude/longitude %s = %g\n", str, arg);
371 * greatcircle - compute the great circle distance in kilometers
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)));
391 "greatcircle lat1 %g long1 %g lat2 %g long2 %g dist %g\n",
392 lat1, long1, lat2, long2, dg);
398 * waveangle - compute the wave angle for the given distance, virtual
399 * height and number of hops.
411 theta = dg / (EARTHRADIUS * (double)(2 * n));
412 delta = atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2)) - theta;
414 printf("waveangle dist %g height %g hops %d angle %g\n",
415 dg, h, n, delta / RADPERDEG);
421 * propdelay - compute the propagation delay
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));
438 printf("propdelay dist %g height %g hops %d time %g\n",
445 * finddelay - find the propagation delay
457 double dg; /* great circle distance */
458 double delta; /* wave angle */
459 int n; /* number of hops */
461 dg = greatcircle(lat1, long1, lat2, long2);
463 printf("great circle distance %g km %g miles\n", dg, dg/MILE);
466 while ((delta = waveangle(dg, h, n)) < 0.0) {
468 printf("tried %d hop%s, no good\n", n, n>1?"s":"");
472 printf("%d hop%s okay, wave angle is %g\n", n, n>1?"s":"",
475 *delay = propdelay(dg, h, n);
480 * satdoit - compute a delay and print it
493 double up_delay,down_delay;
495 satfinddelay(lat1, long1, lat2, long2, &up_delay);
496 satfinddelay(lat3, long3, lat2, long2, &down_delay);
498 printf("%s, delay %g seconds\n", str, up_delay + down_delay);
502 * satfinddelay - calculate the one-way delay time between a ground station
514 double dg; /* great circle distance */
516 dg = greatcircle(lat1, long1, lat2, long2);
518 *delay = satpropdelay(dg);
522 * satpropdelay - calculate the one-way delay time between a ground station
534 theta = dg / (EARTHRADIUS);
535 k1 = EARTHRADIUS * sin(theta);
536 k2 = SATHEIGHT - (EARTHRADIUS * cos(theta));
538 printf("Theta %g k1 %g k2 %g\n", theta, k1, k2);
539 dist = sqrt(k1*k1 + k2*k2);
540 td = dist / LIGHTSPEED;
542 printf("propdelay dist %g height %g time %g\n", dg, dist, td);