contrib/ntp/clockstuff/propdelay.c

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