propdelay.c revision 289764
192494Ssobomax/* propdelay.c,v 3.1 1993/07/06 01:05:24 jbj Exp 292494Ssobomax * propdelay - compute propagation delays 392494Ssobomax * 492494Ssobomax * cc -o propdelay propdelay.c -lm 592494Ssobomax * 692494Ssobomax * "Time and Frequency Users' Manual", NBS Technical Note 695 (1977). 792494Ssobomax */ 892494Ssobomax 992494Ssobomax/* 1092494Ssobomax * This can be used to get a rough idea of the HF propagation delay 1192494Ssobomax * between two points (usually between you and the radio station). 1292494Ssobomax * The usage is 1392494Ssobomax * 1492494Ssobomax * propdelay latitudeA longitudeA latitudeB longitudeB 1592494Ssobomax * 1692494Ssobomax * where points A and B are the locations in question. You obviously 1792494Ssobomax * need to know the latitude and longitude of each of the places. 1892494Ssobomax * The program expects the latitude to be preceded by an 'n' or 's' 1992494Ssobomax * and the longitude to be preceded by an 'e' or 'w'. It understands 2092494Ssobomax * either decimal degrees or degrees:minutes:seconds. Thus to compute 2192494Ssobomax * the delay between the WWVH (21:59:26N, 159:46:00W) and WWV (40:40:49N, 2292494Ssobomax * 105:02:27W) you could use: 2392494Ssobomax * 2492494Ssobomax * propdelay n21:59:26 w159:46 n40:40:49 w105:02:27 2592494Ssobomax * 2692494Ssobomax * By default it prints out a summer (F2 average virtual height 350 km) and 2792494Ssobomax * winter (F2 average virtual height 250 km) number. The results will be 2892494Ssobomax * quite approximate but are about as good as you can do with HF time anyway. 2992494Ssobomax * You might pick a number between the values to use, or use the summer 3092494Ssobomax * value in the summer and switch to the winter value when the static 3192494Ssobomax * above 10 MHz starts to drop off in the fall. You can also use the 3292494Ssobomax * -h switch if you want to specify your own virtual height. 3392494Ssobomax * 3492494Ssobomax * You can also do a 3592494Ssobomax * 3692494Ssobomax * propdelay -W n45:17:47 w75:45:22 3792494Ssobomax * 3892494Ssobomax * to find the propagation delays to WWV and WWVH (from CHU in this 3992494Ssobomax * case), a 4092494Ssobomax * 4192494Ssobomax * propdelay -C n40:40:49 w105:02:27 4292494Ssobomax * 4392494Ssobomax * to find the delays to CHU, and a 4492494Ssobomax * 4592494Ssobomax * propdelay -G n52:03:17 w98:34:18 4692494Ssobomax * 47124572Sjhb * to find delays to GOES via each of the three satellites. 4892494Ssobomax */ 4992494Ssobomax 5092494Ssobomax#ifdef HAVE_CONFIG_H 5192494Ssobomax# include <config.h> 5292494Ssobomax#endif 5392494Ssobomax#include <stdio.h> 5492494Ssobomax#include <string.h> 5592494Ssobomax 5692494Ssobomax#include "ntp_stdlib.h" 5792494Ssobomax 5892494Ssobomaxextern double sin (double); 5992494Ssobomaxextern double cos (double); 6092494Ssobomaxextern double acos (double); 6192494Ssobomaxextern double tan (double); 6292494Ssobomaxextern double atan (double); 6392494Ssobomaxextern double sqrt (double); 6492494Ssobomax 6592494Ssobomax#define STREQ(a, b) (*(a) == *(b) && strcmp((a), (b)) == 0) 6692494Ssobomax 6792494Ssobomax/* 68124571Sjhb * Program constants 6992494Ssobomax */ 70124571Sjhb#define EARTHRADIUS (6370.0) /* raduis of earth (km) */ 7192494Ssobomax#define LIGHTSPEED (299800.0) /* speed of light, km/s */ 7292494Ssobomax#define PI (3.1415926536) 7392494Ssobomax#define RADPERDEG (PI/180.0) /* radians per degree */ 7492494Ssobomax#define MILE (1.609344) /* km in a mile */ 7592494Ssobomax 7692494Ssobomax#define SUMMERHEIGHT (350.0) /* summer height in km */ 77124571Sjhb#define WINTERHEIGHT (250.0) /* winter height in km */ 78124571Sjhb 7992494Ssobomax#define SATHEIGHT (6.6110 * 6378.0) /* geosync satellite height in km 8092494Ssobomax from centre of earth */ 8192494Ssobomax 82124572Sjhb#define WWVLAT "n40:40:49" 83124572Sjhb#define WWVLONG "w105:02:27" 84124572Sjhb 85124572Sjhb#define WWVHLAT "n21:59:26" 86124572Sjhb#define WWVHLONG "w159:46:00" 87124572Sjhb 88124572Sjhb#define CHULAT "n45:17:47" 89124572Sjhb#define CHULONG "w75:45:22" 90124572Sjhb 91124572Sjhb#define GOES_UP_LAT "n37:52:00" 92124572Sjhb#define GOES_UP_LONG "w75:27:00" 93124572Sjhb#define GOES_EAST_LONG "w75:00:00" 94124572Sjhb#define GOES_STBY_LONG "w105:00:00" 95124572Sjhb#define GOES_WEST_LONG "w135:00:00" 96124572Sjhb#define GOES_SAT_LAT "n00:00:00" 97124572Sjhb 98124572Sjhbchar *wwvlat = WWVLAT; 99124572Sjhbchar *wwvlong = WWVLONG; 100124572Sjhb 101124572Sjhbchar *wwvhlat = WWVHLAT; 102124572Sjhbchar *wwvhlong = WWVHLONG; 103124572Sjhb 10492494Ssobomaxchar *chulat = CHULAT; 10592494Ssobomaxchar *chulong = CHULONG; 10692494Ssobomax 10792494Ssobomaxchar *goes_up_lat = GOES_UP_LAT; 10892494Ssobomaxchar *goes_up_long = GOES_UP_LONG; 10992494Ssobomaxchar *goes_east_long = GOES_EAST_LONG; 11092494Ssobomaxchar *goes_stby_long = GOES_STBY_LONG; 11192494Ssobomaxchar *goes_west_long = GOES_WEST_LONG; 11292494Ssobomaxchar *goes_sat_lat = GOES_SAT_LAT; 11392494Ssobomax 11492494Ssobomaxint hflag = 0; 11592494Ssobomaxint Wflag = 0; 11692494Ssobomaxint Cflag = 0; 11792494Ssobomaxint Gflag = 0; 11892494Ssobomaxint height; 11992494Ssobomax 12092494Ssobomaxchar const *progname; 12192494Ssobomax 12292494Ssobomaxstatic void doit (double, double, double, double, double, char *); 12392494Ssobomaxstatic double latlong (char *, int); 12492494Ssobomaxstatic double greatcircle (double, double, double, double); 12592494Ssobomaxstatic double waveangle (double, double, int); 12692494Ssobomaxstatic double propdelay (double, double, int); 12792494Ssobomaxstatic int finddelay (double, double, double, double, double, double *); 12892494Ssobomaxstatic void satdoit (double, double, double, double, double, double, char *); 12992494Ssobomaxstatic void satfinddelay (double, double, double, double, double *); 13092494Ssobomaxstatic double satpropdelay (double); 13192494Ssobomax 13292494Ssobomax/* 13392494Ssobomax * main - parse arguments and handle options 13492494Ssobomax */ 13592494Ssobomaxint 13692494Ssobomaxmain( 13792494Ssobomax int argc, 13892494Ssobomax char *argv[] 13992494Ssobomax ) 14092494Ssobomax{ 14192494Ssobomax int c; 14292494Ssobomax int errflg = 0; 14392494Ssobomax double lat1, long1; 14492494Ssobomax double lat2, long2; 14592494Ssobomax double lat3, long3; 14692494Ssobomax 14792494Ssobomax init_lib(); 14892494Ssobomax 14992494Ssobomax progname = argv[0]; 15092494Ssobomax while ((c = ntp_getopt(argc, argv, "dh:CWG")) != EOF) 15192494Ssobomax switch (c) { 15292494Ssobomax case 'd': 15392494Ssobomax ++debug; 15492494Ssobomax break; 15592494Ssobomax case 'h': 15692494Ssobomax hflag++; 15792494Ssobomax height = atof(ntp_optarg); 15892494Ssobomax if (height <= 0.0) { 15992494Ssobomax (void) fprintf(stderr, "height %s unlikely\n", 16092494Ssobomax ntp_optarg); 16192494Ssobomax errflg++; 16292494Ssobomax } 16392494Ssobomax break; 16492494Ssobomax case 'C': 165124572Sjhb Cflag++; 16692494Ssobomax break; 16792494Ssobomax case 'W': 16892494Ssobomax Wflag++; 169124572Sjhb break; 170124572Sjhb case 'G': 171124572Sjhb Gflag++; 172124572Sjhb break; 173124572Sjhb default: 17492494Ssobomax errflg++; 17592494Ssobomax break; 17692494Ssobomax } 17792494Ssobomax if (errflg || (!(Cflag || Wflag || Gflag) && ntp_optind+4 != argc) || 17892494Ssobomax ((Cflag || Wflag || Gflag) && ntp_optind+2 != argc)) { 17992494Ssobomax (void) fprintf(stderr, 18092494Ssobomax "usage: %s [-d] [-h height] lat1 long1 lat2 long2\n", 18192494Ssobomax progname); 18292494Ssobomax (void) fprintf(stderr," - or -\n"); 18392494Ssobomax (void) fprintf(stderr, 18492494Ssobomax "usage: %s -CWG [-d] lat long\n", 18592494Ssobomax progname); 18692494Ssobomax exit(2); 18792494Ssobomax } 18892494Ssobomax 18992494Ssobomax 19092494Ssobomax if (!(Cflag || Wflag || Gflag)) { 19192494Ssobomax lat1 = latlong(argv[ntp_optind], 1); 19292494Ssobomax long1 = latlong(argv[ntp_optind + 1], 0); 19392494Ssobomax lat2 = latlong(argv[ntp_optind + 2], 1); 19492494Ssobomax long2 = latlong(argv[ntp_optind + 3], 0); 19592494Ssobomax if (hflag) { 19692494Ssobomax doit(lat1, long1, lat2, long2, height, ""); 19792494Ssobomax } else { 19892494Ssobomax doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 19992494Ssobomax "summer propagation, "); 20092494Ssobomax doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 20192494Ssobomax "winter propagation, "); 20292494Ssobomax } 20392494Ssobomax } else if (Wflag) { 20492494Ssobomax /* 20592494Ssobomax * Compute delay from WWV 20692494Ssobomax */ 20792494Ssobomax lat1 = latlong(argv[ntp_optind], 1); 20892494Ssobomax long1 = latlong(argv[ntp_optind + 1], 0); 20992494Ssobomax lat2 = latlong(wwvlat, 1); 210136093Sstefanf long2 = latlong(wwvlong, 0); 21192494Ssobomax if (hflag) { 21292494Ssobomax doit(lat1, long1, lat2, long2, height, "WWV "); 21392494Ssobomax } else { 21492494Ssobomax doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 21592494Ssobomax "WWV summer propagation, "); 21692494Ssobomax doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 21792494Ssobomax "WWV winter propagation, "); 21892494Ssobomax } 21992494Ssobomax 22092494Ssobomax /* 221124572Sjhb * Compute delay from WWVH 222124572Sjhb */ 22392494Ssobomax lat2 = latlong(wwvhlat, 1); 22492494Ssobomax long2 = latlong(wwvhlong, 0); 22592494Ssobomax if (hflag) { 22692494Ssobomax doit(lat1, long1, lat2, long2, height, "WWVH "); 22792494Ssobomax } else { 22892494Ssobomax doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 22992494Ssobomax "WWVH summer propagation, "); 23092494Ssobomax doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 23192494Ssobomax "WWVH winter propagation, "); 23292494Ssobomax } 23392494Ssobomax } else if (Cflag) { 23492494Ssobomax lat1 = latlong(argv[ntp_optind], 1); 23592494Ssobomax long1 = latlong(argv[ntp_optind + 1], 0); 23692494Ssobomax lat2 = latlong(chulat, 1); 23792494Ssobomax long2 = latlong(chulong, 0); 23892494Ssobomax if (hflag) { 23992494Ssobomax doit(lat1, long1, lat2, long2, height, "CHU "); 24092494Ssobomax } else { 24192494Ssobomax doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 24292494Ssobomax "CHU summer propagation, "); 24392494Ssobomax doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 24492494Ssobomax "CHU winter propagation, "); 24592494Ssobomax } 24692494Ssobomax } else if (Gflag) { 24792494Ssobomax lat1 = latlong(goes_up_lat, 1); 24892494Ssobomax long1 = latlong(goes_up_long, 0); 24992494Ssobomax lat3 = latlong(argv[ntp_optind], 1); 25092494Ssobomax long3 = latlong(argv[ntp_optind + 1], 0); 25192494Ssobomax 25292494Ssobomax lat2 = latlong(goes_sat_lat, 1); 253124811Sjhb 254124811Sjhb long2 = latlong(goes_west_long, 0); 255124811Sjhb satdoit(lat1, long1, lat2, long2, lat3, long3, 25692494Ssobomax "GOES Delay via WEST"); 25792494Ssobomax 25892494Ssobomax long2 = latlong(goes_stby_long, 0); 25992494Ssobomax satdoit(lat1, long1, lat2, long2, lat3, long3, 26092494Ssobomax "GOES Delay via STBY"); 26192494Ssobomax 26292494Ssobomax long2 = latlong(goes_east_long, 0); 26392494Ssobomax satdoit(lat1, long1, lat2, long2, lat3, long3, 26492494Ssobomax "GOES Delay via EAST"); 26592494Ssobomax 26692494Ssobomax } 267124811Sjhb exit(0); 268124811Sjhb} 26992494Ssobomax 270124811Sjhb 27192494Ssobomax/* 27292494Ssobomax * doit - compute a delay and print it 27392494Ssobomax */ 27492494Ssobomaxstatic void 27592494Ssobomaxdoit( 27692494Ssobomax double lat1, 27792494Ssobomax double long1, 27892494Ssobomax double lat2, 27992494Ssobomax double long2, 28092494Ssobomax double h, 28192494Ssobomax char *str 28292494Ssobomax ) 28392494Ssobomax{ 28492494Ssobomax int hops; 28592494Ssobomax double delay; 28692494Ssobomax 28792494Ssobomax hops = finddelay(lat1, long1, lat2, long2, h, &delay); 28892494Ssobomax printf("%sheight %g km, hops %d, delay %g seconds\n", 28992494Ssobomax str, h, hops, delay); 29092494Ssobomax} 291124811Sjhb 292124811Sjhb 29392494Ssobomax/* 294124811Sjhb * latlong - decode a latitude/longitude value 29592494Ssobomax */ 29692494Ssobomaxstatic double 29792494Ssobomaxlatlong( 29892494Ssobomax char *str, 29992494Ssobomax int islat 30092494Ssobomax ) 30192494Ssobomax{ 30292494Ssobomax register char *cp; 30392494Ssobomax register char *bp; 30492494Ssobomax double arg; 30592494Ssobomax double divby; 30692494Ssobomax int isneg; 30792494Ssobomax char buf[32]; 30892494Ssobomax char *colon; 30992494Ssobomax 31092494Ssobomax if (islat) { 31192494Ssobomax /* 31292494Ssobomax * 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 */ 377static double 378greatcircle( 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 */ 405static double 406waveangle( 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 */ 427static double 428propdelay( 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 */ 451static int 452finddelay( 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 */ 486static void 487satdoit( 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 */ 509static void 510satfinddelay( 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 */ 529static double 530satpropdelay( 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