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 58extern double sin (double); 59extern double cos (double); 60extern double acos (double); 61extern double tan (double); 62extern double atan (double); 63extern 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 98char *wwvlat = WWVLAT; 99char *wwvlong = WWVLONG; 100 101char *wwvhlat = WWVHLAT; 102char *wwvhlong = WWVHLONG; 103 104char *chulat = CHULAT; 105char *chulong = CHULONG; 106 107char *goes_up_lat = GOES_UP_LAT; 108char *goes_up_long = GOES_UP_LONG; 109char *goes_east_long = GOES_EAST_LONG; 110char *goes_stby_long = GOES_STBY_LONG; 111char *goes_west_long = GOES_WEST_LONG; 112char *goes_sat_lat = GOES_SAT_LAT; 113 114int hflag = 0; 115int Wflag = 0; 116int Cflag = 0; 117int Gflag = 0; 118int height; 119 120char const *progname; 121 122static void doit (double, double, double, double, double, char *); 123static double latlong (char *, int); 124static double greatcircle (double, double, double, double); 125static double waveangle (double, double, int); 126static double propdelay (double, double, int); 127static int finddelay (double, double, double, double, double, double *); 128static void satdoit (double, double, double, double, double, double, char *); 129static void satfinddelay (double, double, double, double, double *); 130static double satpropdelay (double); 131 132/* 133 * main - parse arguments and handle options 134 */ 135int 136main( 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 */ 274static void 275doit( 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 */ 296static double 297latlong( 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 */ 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