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