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

Revision 16192, 11.8 KB checked in by ghudson, 24 years ago (diff)
This commit was generated by cvs2svn to compensate for changes in r16191, which included commits to RCS files with non-trunk default branches.
Line 
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     (double);
56extern  double  cos     (double);
57extern  double  acos    (double);
58extern  double  tan     (double);
59extern  double  atan    (double);
60extern  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
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            (double, double, double, double, double, char *);
121static  double  latlong         (char *, int);
122static  double  greatcircle     (double, double, double, double);
123static  double  waveangle       (double, double, int);
124static  double  propdelay       (double, double, int);
125static  int     finddelay       (double, double, double, double, double, double *);
126static  void    satdoit         (double, double, double, double, double, double, char *);
127static  void    satfinddelay    (double, double, double, double, double *);
128static  double  satpropdelay    (double);
129
130/*
131 * main - parse arguments and handle options
132 */
133int
134main(
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 */
270static void
271doit(
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 */
292static double
293latlong(
294        char *str,
295        int islat
296        )
297{
298        register char *cp;
299        register char *bp;
300        double arg;
301        double div;
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                div = 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) / div;
352                        div = 3600.0;
353                }
354                if (*cp != '\0')
355                    arg += atof(cp) / div;
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 */
373static double
374greatcircle(
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 */
401static double
402waveangle(
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 */
423static double
424propdelay(
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 */
447static int
448finddelay(
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 */
482static void
483satdoit(
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 */
505static void
506satfinddelay(
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 */
525static double
526satpropdelay(
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}
Note: See TracBrowser for help on using the repository browser.