dist/ntp/clockstuff/propdelay.c

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