1 /* $NetBSD: propdelay.c,v 1.1.1.1 2009/12/13 16:53:41 kardel 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 } 547