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 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 div; 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 div = 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) / div; 352 div = 3600.0; 353 } 354 if (*cp != '\0') 355 arg += atof(cp) / div; 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 } 545