source: trunk/third/xntp/clockstuff/propdelay.c @ 10832

Revision 10832, 12.0 KB checked in by brlewis, 27 years ago (diff)
This commit was generated by cvs2svn to compensate for changes in r10831, which included commits to RCS files with non-trunk default branches.
RevLine 
[10831]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
55extern  double  sin     P((double));
56extern  double  cos     P((double));
57extern  double  acos    P((double));
58extern  double  tan     P((double));
59extern  double  atan    P((double));
60extern  double  sqrt    P((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
95char *wwvlat = WWVLAT;
96char *wwvlong = WWVLONG;
97
98char *wwvhlat = WWVHLAT;
99char *wwvhlong = WWVHLONG;
100
101char *chulat = CHULAT;
102char *chulong = CHULONG;
103
104char *goes_up_lat = GOES_UP_LAT;
105char *goes_up_long = GOES_UP_LONG;
106char *goes_east_long = GOES_EAST_LONG;
107char *goes_stby_long = GOES_STBY_LONG;
108char *goes_west_long = GOES_WEST_LONG;
109char *goes_sat_lat = GOES_SAT_LAT;
110
111int hflag = 0;
112int Wflag = 0;
113int Cflag = 0;
114int Gflag = 0;
115int height;
116
117char *progname;
118int debug;
119
120static  void    doit            P((double, double, double, double, double, char *));
121static  double  latlong         P((char *, int));
122static  double  greatcircle     P((double, double, double, double));
123static  double  waveangle       P((double, double, int));
124static  double  propdelay       P((double, double, int));
125static  int     finddelay       P((double, double, double, double, double, double *));
126static  void    satdoit         P((double, double, double, double, double, double, char *));
127static  void    satfinddelay    P((double, double, double, double, double *));
128static  double  satpropdelay    P((double));
129
130/*
131 * main - parse arguments and handle options
132 */
133void
134main(argc, argv)
135int argc;
136char *argv[];
137{
138        int c;
139        int errflg = 0;
140        double lat1, long1;
141        double lat2, long2;
142        double lat3, long3;
143        extern int ntp_optind;
144        extern char *ntp_optarg;
145
146        progname = argv[0];
147        while ((c = ntp_getopt(argc, argv, "dh:CWG")) != EOF)
148                switch (c) {
149                case 'd':
150                        ++debug;
151                        break;
152                case 'h':
153                        hflag++;
154                        height = atof(ntp_optarg);
155                        if (height <= 0.0) {
156                                (void) fprintf(stderr, "height %s unlikely\n",
157                                               ntp_optarg);
158                                errflg++;
159                        }
160                        break;
161                case 'C':
162                        Cflag++;
163                        break;
164                case 'W':
165                        Wflag++;
166                        break;
167                case 'G':
168                        Gflag++;
169                        break;
170                default:
171                        errflg++;
172                        break;
173                }
174        if (errflg || (!(Cflag || Wflag || Gflag) && ntp_optind+4 != argc) ||
175            ((Cflag || Wflag || Gflag) && ntp_optind+2 != argc)) {
176                (void) fprintf(stderr,
177                    "usage: %s [-d] [-h height] lat1 long1 lat2 long2\n",
178                    progname);
179                (void) fprintf(stderr," - or -\n");
180                (void) fprintf(stderr,
181                    "usage: %s -CWG [-d] lat long\n",
182                    progname);
183                exit(2);
184        }
185
186                   
187        if (!(Cflag || Wflag || Gflag)) {
188                lat1 = latlong(argv[ntp_optind], 1);
189                long1 = latlong(argv[ntp_optind + 1], 0);
190                lat2 = latlong(argv[ntp_optind + 2], 1);
191                long2 = latlong(argv[ntp_optind + 3], 0);
192                if (hflag) {
193                        doit(lat1, long1, lat2, long2, height, "");
194                } else {
195                        doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
196                            "summer propagation, ");
197                        doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
198                            "winter propagation, ");
199                }
200        } else if (Wflag) {
201                /*
202                 * Compute delay from WWV
203                 */
204                lat1 = latlong(argv[ntp_optind], 1);
205                long1 = latlong(argv[ntp_optind + 1], 0);
206                lat2 = latlong(wwvlat, 1);
207                long2 = latlong(wwvlong, 0);
208                if (hflag) {
209                        doit(lat1, long1, lat2, long2, height, "WWV  ");
210                } else {
211                        doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
212                            "WWV  summer propagation, ");
213                        doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
214                            "WWV  winter propagation, ");
215                }
216
217                /*
218                 * Compute delay from WWVH
219                 */
220                lat2 = latlong(wwvhlat, 1);
221                long2 = latlong(wwvhlong, 0);
222                if (hflag) {
223                        doit(lat1, long1, lat2, long2, height, "WWVH ");
224                } else {
225                        doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
226                            "WWVH summer propagation, ");
227                        doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
228                            "WWVH winter propagation, ");
229                }
230        } else if (Cflag) {
231                lat1 = latlong(argv[ntp_optind], 1);
232                long1 = latlong(argv[ntp_optind + 1], 0);
233                lat2 = latlong(chulat, 1);
234                long2 = latlong(chulong, 0);
235                if (hflag) {
236                        doit(lat1, long1, lat2, long2, height, "CHU ");
237                } else {
238                        doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
239                            "CHU summer propagation, ");
240                        doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
241                            "CHU winter propagation, ");
242                }
243        } else if (Gflag) {
244                lat1 = latlong(goes_up_lat, 1);
245                long1 = latlong(goes_up_long, 0);
246                lat3 = latlong(argv[ntp_optind], 1);
247                long3 = latlong(argv[ntp_optind + 1], 0);
248
249                lat2 = latlong(goes_sat_lat, 1);
250
251                long2 = latlong(goes_west_long, 0);
252                satdoit(lat1, long1, lat2, long2, lat3, long3,
253                        "GOES Delay via WEST");
254
255                long2 = latlong(goes_stby_long, 0);
256                satdoit(lat1, long1, lat2, long2, lat3, long3,
257                        "GOES Delay via STBY");
258
259                long2 = latlong(goes_east_long, 0);
260                satdoit(lat1, long1, lat2, long2, lat3, long3,
261                        "GOES Delay via EAST");
262
263        }
264        exit(0);
265}
266
267
268/*
269 * doit - compute a delay and print it
270 */
271static void
272doit(lat1, long1, lat2, long2, h, str)
273        double lat1;
274        double long1;
275        double lat2;
276        double long2;
277        double h;
278        char *str;
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 */
292static double
293latlong(str, islat)
294        char *str;
295        int islat;
296{
297        register char *cp;
298        register char *bp;
299        double arg;
300        double div;
301        int isneg;
302        char buf[32];
303        char *colon;
304
305        if (islat) {
306                /*
307                 * Must be north or south
308                 */
309                if (*str == 'N' || *str == 'n')
310                        isneg = 0;
311                else if (*str == 'S' || *str == 's')
312                        isneg = 1;
313                else
314                        isneg = -1;
315        } else {
316                /*
317                 * East is positive, west is negative
318                 */
319                if (*str == 'E' || *str == 'e')
320                        isneg = 0;
321                else if (*str == 'W' || *str == 'w')
322                        isneg = 1;
323                else
324                        isneg = -1;
325        }
326
327        if (isneg >= 0)
328                str++;
329
330        colon = strchr(str, ':');
331        if (colon != NULL) {
332                /*
333                 * in hhh:mm:ss form
334                 */
335                cp = str;
336                bp = buf;
337                while (cp < colon)
338                        *bp++ = *cp++;
339                *bp = '\0';
340                cp++;
341                arg = atof(buf);
342                div = 60.0;
343                colon = strchr(cp, ':');
344                if (colon != NULL) {
345                        bp = buf;
346                        while (cp < colon)
347                                *bp++ = *cp++;
348                        *bp = '\0';
349                        cp++;
350                        arg += atof(buf) / div;
351                        div = 3600.0;
352                }
353                if (*cp != '\0')
354                        arg += atof(cp) / div;
355        } else {
356                arg = atof(str);
357        }
358
359        if (isneg == 1)
360                arg = -arg;
361
362        if (debug > 2)
363                (void) printf("latitude/longitude %s = %g\n", str, arg);
364
365        return arg;
366}
367
368
369/*
370 * greatcircle - compute the great circle distance in kilometers
371 */
372static double
373greatcircle(lat1, long1, lat2, long2)
374        double lat1;
375        double long1;
376        double lat2;
377        double long2;
378{
379        double dg;
380        double l1r, l2r;
381
382        l1r = lat1 * RADPERDEG;
383        l2r = lat2 * RADPERDEG;
384        dg = EARTHRADIUS * acos(
385            (cos(l1r) * cos(l2r) * cos((long2-long1)*RADPERDEG))
386            + (sin(l1r) * sin(l2r)));
387        if (debug >= 2)
388                printf(
389                    "greatcircle lat1 %g long1 %g lat2 %g long2 %g dist %g\n",
390                    lat1, long1, lat2, long2, dg);
391        return dg;
392}
393
394
395/*
396 * waveangle - compute the wave angle for the given distance, virtual
397 *             height and number of hops.
398 */
399static double
400waveangle(dg, h, n)
401        double dg;
402        double h;
403        int n;
404{
405        double theta;
406        double delta;
407
408        theta = dg / (EARTHRADIUS * (double)(2 * n));
409        delta = atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2)) - theta;
410        if (debug >= 2)
411                printf("waveangle dist %g height %g hops %d angle %g\n",
412                    dg, h, n, delta / RADPERDEG);
413        return delta;
414}
415
416
417/*
418 * propdelay - compute the propagation delay
419 */
420static double
421propdelay(dg, h, n)
422        double dg;
423        double h;
424        int n;
425{
426        double phi;
427        double theta;
428        double td;
429
430        theta = dg / (EARTHRADIUS * (double)(2 * n));
431        phi = (PI/2.0) - atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2));
432        td = dg / (LIGHTSPEED * sin(phi));
433        if (debug >= 2)
434                printf("propdelay dist %g height %g hops %d time %g\n",
435                    dg, h, n, td);
436        return td;
437}
438
439
440/*
441 * finddelay - find the propagation delay
442 */
443static int
444finddelay(lat1, long1, lat2, long2, h, delay)
445        double lat1;
446        double long1;
447        double lat2;
448        double long2;
449        double h;
450        double *delay;
451{
452        double dg;      /* great circle distance */
453        double delta;   /* wave angle */
454        int n;          /* number of hops */
455
456        dg = greatcircle(lat1, long1, lat2, long2);
457        if (debug)
458                printf("great circle distance %g km %g miles\n", dg, dg/MILE);
459       
460        n = 1;
461        while ((delta = waveangle(dg, h, n)) < 0.0) {
462                if (debug)
463                        printf("tried %d hop%s, no good\n", n, n>1?"s":"");
464                n++;
465        }
466        if (debug)
467                printf("%d hop%s okay, wave angle is %g\n", n, n>1?"s":"",
468                    delta / RADPERDEG);
469
470        *delay = propdelay(dg, h, n);
471        return n;
472}
473
474/*
475 * satdoit - compute a delay and print it
476 */
477static void
478satdoit(lat1, long1, lat2, long2, lat3, long3, str)
479        double lat1;
480        double long1;
481        double lat2;
482        double long2;
483        double lat3;
484        double long3;
485        char *str;
486{
487        double up_delay,down_delay;
488
489        satfinddelay(lat1, long1, lat2, long2, &up_delay);
490        satfinddelay(lat3, long3, lat2, long2, &down_delay);
491
492        printf("%s, delay %g seconds\n", str, up_delay + down_delay);
493}
494
495/*
496 * satfinddelay - calculate the one-way delay time between a ground station
497 * and a satellite
498 */
499static void
500satfinddelay(lat1, long1, lat2, long2, delay)
501        double lat1;
502        double long1;
503        double lat2;
504        double long2;
505        double *delay;
506{
507        double dg;      /* great circle distance */
508
509        dg = greatcircle(lat1, long1, lat2, long2);
510
511        *delay = satpropdelay(dg);
512}
513
514/*
515 * satpropdelay - calculate the one-way delay time between a ground station
516 * and a satellite
517 */
518static double
519satpropdelay(dg)
520        double dg;
521{
522        double k1, k2, dist;
523        double theta;
524        double td;
525
526        theta = dg / (EARTHRADIUS);
527        k1 = EARTHRADIUS * sin(theta);
528        k2 = SATHEIGHT - (EARTHRADIUS * cos(theta));
529        if (debug >= 2)
530                printf("Theta %g k1 %g k2 %g\n", theta, k1, k2);
531        dist = sqrt(k1*k1 + k2*k2);
532        td = dist / LIGHTSPEED;
533        if (debug >= 2)
534                printf("propdelay dist %g height %g time %g\n", dg, dist, td);
535        return td;
536}
Note: See TracBrowser for help on using the repository browser.