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 #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 }