1/* $NetBSD: propdelay.c,v 1.6 2020/05/25 20:47:19 christos 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#ifdef HAVE_CONFIG_H 53# include <config.h> 54#endif 55#include <stdio.h> 56#include <string.h> 57 58#include "ntp_stdlib.h" 59 60extern double sin (double); 61extern double cos (double); 62extern double acos (double); 63extern double tan (double); 64extern double atan (double); 65extern double sqrt (double); 66 67#define STREQ(a, b) (*(a) == *(b) && strcmp((a), (b)) == 0) 68 69/* 70 * Program constants 71 */ 72#define EARTHRADIUS (6370.0) /* raduis of earth (km) */ 73#define LIGHTSPEED (299800.0) /* speed of light, km/s */ 74#define PI (3.1415926536) 75#define RADPERDEG (PI/180.0) /* radians per degree */ 76#define MILE (1.609344) /* km in a mile */ 77 78#define SUMMERHEIGHT (350.0) /* summer height in km */ 79#define WINTERHEIGHT (250.0) /* winter height in km */ 80 81#define SATHEIGHT (6.6110 * 6378.0) /* geosync satellite height in km 82 from centre of earth */ 83 84#define WWVLAT "n40:40:49" 85#define WWVLONG "w105:02:27" 86 87#define WWVHLAT "n21:59:26" 88#define WWVHLONG "w159:46:00" 89 90#define CHULAT "n45:17:47" 91#define CHULONG "w75:45:22" 92 93#define GOES_UP_LAT "n37:52:00" 94#define GOES_UP_LONG "w75:27:00" 95#define GOES_EAST_LONG "w75:00:00" 96#define GOES_STBY_LONG "w105:00:00" 97#define GOES_WEST_LONG "w135:00:00" 98#define GOES_SAT_LAT "n00:00:00" 99 100char *wwvlat = WWVLAT; 101char *wwvlong = WWVLONG; 102 103char *wwvhlat = WWVHLAT; 104char *wwvhlong = WWVHLONG; 105 106char *chulat = CHULAT; 107char *chulong = CHULONG; 108 109char *goes_up_lat = GOES_UP_LAT; 110char *goes_up_long = GOES_UP_LONG; 111char *goes_east_long = GOES_EAST_LONG; 112char *goes_stby_long = GOES_STBY_LONG; 113char *goes_west_long = GOES_WEST_LONG; 114char *goes_sat_lat = GOES_SAT_LAT; 115 116int hflag = 0; 117int Wflag = 0; 118int Cflag = 0; 119int Gflag = 0; 120int height; 121 122char const *progname; 123 124static void doit (double, double, double, double, double, char *); 125static double latlong (char *, int); 126static double greatcircle (double, double, double, double); 127static double waveangle (double, double, int); 128static double propdelay (double, double, int); 129static int finddelay (double, double, double, double, double, double *); 130static void satdoit (double, double, double, double, double, double, char *); 131static void satfinddelay (double, double, double, double, double *); 132static double satpropdelay (double); 133 134/* 135 * main - parse arguments and handle options 136 */ 137int 138main( 139 int argc, 140 char *argv[] 141 ) 142{ 143 int c; 144 int errflg = 0; 145 double lat1, long1; 146 double lat2, long2; 147 double lat3, long3; 148 149 init_lib(); 150 151 progname = argv[0]; 152 while ((c = ntp_getopt(argc, argv, "dh:CWG")) != EOF) 153 switch (c) { 154 case 'd': 155 ++debug; 156 break; 157 case 'h': 158 hflag++; 159 height = atof(ntp_optarg); 160 if (height <= 0.0) { 161 (void) fprintf(stderr, "height %s unlikely\n", 162 ntp_optarg); 163 errflg++; 164 } 165 break; 166 case 'C': 167 Cflag++; 168 break; 169 case 'W': 170 Wflag++; 171 break; 172 case 'G': 173 Gflag++; 174 break; 175 default: 176 errflg++; 177 break; 178 } 179 if (errflg || (!(Cflag || Wflag || Gflag) && ntp_optind+4 != argc) || 180 ((Cflag || Wflag || Gflag) && ntp_optind+2 != argc)) { 181 (void) fprintf(stderr, 182 "usage: %s [-d] [-h height] lat1 long1 lat2 long2\n", 183 progname); 184 (void) fprintf(stderr," - or -\n"); 185 (void) fprintf(stderr, 186 "usage: %s -CWG [-d] lat long\n", 187 progname); 188 exit(2); 189 } 190 191 192 if (!(Cflag || Wflag || Gflag)) { 193 lat1 = latlong(argv[ntp_optind], 1); 194 long1 = latlong(argv[ntp_optind + 1], 0); 195 lat2 = latlong(argv[ntp_optind + 2], 1); 196 long2 = latlong(argv[ntp_optind + 3], 0); 197 if (hflag) { 198 doit(lat1, long1, lat2, long2, height, ""); 199 } else { 200 doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 201 "summer propagation, "); 202 doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 203 "winter propagation, "); 204 } 205 } else if (Wflag) { 206 /* 207 * Compute delay from WWV 208 */ 209 lat1 = latlong(argv[ntp_optind], 1); 210 long1 = latlong(argv[ntp_optind + 1], 0); 211 lat2 = latlong(wwvlat, 1); 212 long2 = latlong(wwvlong, 0); 213 if (hflag) { 214 doit(lat1, long1, lat2, long2, height, "WWV "); 215 } else { 216 doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 217 "WWV summer propagation, "); 218 doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 219 "WWV winter propagation, "); 220 } 221 222 /* 223 * Compute delay from WWVH 224 */ 225 lat2 = latlong(wwvhlat, 1); 226 long2 = latlong(wwvhlong, 0); 227 if (hflag) { 228 doit(lat1, long1, lat2, long2, height, "WWVH "); 229 } else { 230 doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 231 "WWVH summer propagation, "); 232 doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 233 "WWVH winter propagation, "); 234 } 235 } else if (Cflag) { 236 lat1 = latlong(argv[ntp_optind], 1); 237 long1 = latlong(argv[ntp_optind + 1], 0); 238 lat2 = latlong(chulat, 1); 239 long2 = latlong(chulong, 0); 240 if (hflag) { 241 doit(lat1, long1, lat2, long2, height, "CHU "); 242 } else { 243 doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 244 "CHU summer propagation, "); 245 doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 246 "CHU winter propagation, "); 247 } 248 } else if (Gflag) { 249 lat1 = latlong(goes_up_lat, 1); 250 long1 = latlong(goes_up_long, 0); 251 lat3 = latlong(argv[ntp_optind], 1); 252 long3 = latlong(argv[ntp_optind + 1], 0); 253 254 lat2 = latlong(goes_sat_lat, 1); 255 256 long2 = latlong(goes_west_long, 0); 257 satdoit(lat1, long1, lat2, long2, lat3, long3, 258 "GOES Delay via WEST"); 259 260 long2 = latlong(goes_stby_long, 0); 261 satdoit(lat1, long1, lat2, long2, lat3, long3, 262 "GOES Delay via STBY"); 263 264 long2 = latlong(goes_east_long, 0); 265 satdoit(lat1, long1, lat2, long2, lat3, long3, 266 "GOES Delay via EAST"); 267 268 } 269 exit(0); 270} 271 272 273/* 274 * doit - compute a delay and print it 275 */ 276static void 277doit( 278 double lat1, 279 double long1, 280 double lat2, 281 double long2, 282 double h, 283 char *str 284 ) 285{ 286 int hops; 287 double delay; 288 289 hops = finddelay(lat1, long1, lat2, long2, h, &delay); 290 printf("%sheight %g km, hops %d, delay %g seconds\n", 291 str, h, hops, delay); 292} 293 294 295/* 296 * latlong - decode a latitude/longitude value 297 */ 298static double 299latlong( 300 char *str, 301 int islat 302 ) 303{ 304 register char *cp; 305 register char *bp; 306 double arg; 307 double divby; 308 int isneg; 309 char buf[32]; 310 char *colon; 311 312 if (islat) { 313 /* 314 * Must be north or south 315 */ 316 if (*str == 'N' || *str == 'n') 317 isneg = 0; 318 else if (*str == 'S' || *str == 's') 319 isneg = 1; 320 else 321 isneg = -1; 322 } else { 323 /* 324 * East is positive, west is negative 325 */ 326 if (*str == 'E' || *str == 'e') 327 isneg = 0; 328 else if (*str == 'W' || *str == 'w') 329 isneg = 1; 330 else 331 isneg = -1; 332 } 333 334 if (isneg >= 0) 335 str++; 336 337 colon = strchr(str, ':'); 338 if (colon != NULL) { 339 /* 340 * in hhh:mm:ss form 341 */ 342 cp = str; 343 bp = buf; 344 while (cp < colon) 345 *bp++ = *cp++; 346 *bp = '\0'; 347 cp++; 348 arg = atof(buf); 349 divby = 60.0; 350 colon = strchr(cp, ':'); 351 if (colon != NULL) { 352 bp = buf; 353 while (cp < colon) 354 *bp++ = *cp++; 355 *bp = '\0'; 356 cp++; 357 arg += atof(buf) / divby; 358 divby = 3600.0; 359 } 360 if (*cp != '\0') 361 arg += atof(cp) / divby; 362 } else { 363 arg = atof(str); 364 } 365 366 if (isneg == 1) 367 arg = -arg; 368 369 if (debug > 2) 370 (void) printf("latitude/longitude %s = %g\n", str, arg); 371 372 return arg; 373} 374 375 376/* 377 * greatcircle - compute the great circle distance in kilometers 378 */ 379static double 380greatcircle( 381 double lat1, 382 double long1, 383 double lat2, 384 double long2 385 ) 386{ 387 double dg; 388 double l1r, l2r; 389 390 l1r = lat1 * RADPERDEG; 391 l2r = lat2 * RADPERDEG; 392 dg = EARTHRADIUS * acos( 393 (cos(l1r) * cos(l2r) * cos((long2-long1)*RADPERDEG)) 394 + (sin(l1r) * sin(l2r))); 395 if (debug >= 2) 396 printf( 397 "greatcircle lat1 %g long1 %g lat2 %g long2 %g dist %g\n", 398 lat1, long1, lat2, long2, dg); 399 return dg; 400} 401 402 403/* 404 * waveangle - compute the wave angle for the given distance, virtual 405 * height and number of hops. 406 */ 407static double 408waveangle( 409 double dg, 410 double h, 411 int n 412 ) 413{ 414 double theta; 415 double delta; 416 417 theta = dg / (EARTHRADIUS * (double)(2 * n)); 418 delta = atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2)) - theta; 419 if (debug >= 2) 420 printf("waveangle dist %g height %g hops %d angle %g\n", 421 dg, h, n, delta / RADPERDEG); 422 return delta; 423} 424 425 426/* 427 * propdelay - compute the propagation delay 428 */ 429static double 430propdelay( 431 double dg, 432 double h, 433 int n 434 ) 435{ 436 double phi; 437 double theta; 438 double td; 439 440 theta = dg / (EARTHRADIUS * (double)(2 * n)); 441 phi = (PI/2.0) - atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2)); 442 td = dg / (LIGHTSPEED * sin(phi)); 443 if (debug >= 2) 444 printf("propdelay dist %g height %g hops %d time %g\n", 445 dg, h, n, td); 446 return td; 447} 448 449 450/* 451 * finddelay - find the propagation delay 452 */ 453static int 454finddelay( 455 double lat1, 456 double long1, 457 double lat2, 458 double long2, 459 double h, 460 double *delay 461 ) 462{ 463 double dg; /* great circle distance */ 464 double delta; /* wave angle */ 465 int n; /* number of hops */ 466 467 dg = greatcircle(lat1, long1, lat2, long2); 468 if (debug) 469 printf("great circle distance %g km %g miles\n", dg, dg/MILE); 470 471 n = 1; 472 while ((delta = waveangle(dg, h, n)) < 0.0) { 473 if (debug) 474 printf("tried %d hop%s, no good\n", n, n>1?"s":""); 475 n++; 476 } 477 if (debug) 478 printf("%d hop%s okay, wave angle is %g\n", n, n>1?"s":"", 479 delta / RADPERDEG); 480 481 *delay = propdelay(dg, h, n); 482 return n; 483} 484 485/* 486 * satdoit - compute a delay and print it 487 */ 488static void 489satdoit( 490 double lat1, 491 double long1, 492 double lat2, 493 double long2, 494 double lat3, 495 double long3, 496 char *str 497 ) 498{ 499 double up_delay,down_delay; 500 501 satfinddelay(lat1, long1, lat2, long2, &up_delay); 502 satfinddelay(lat3, long3, lat2, long2, &down_delay); 503 504 printf("%s, delay %g seconds\n", str, up_delay + down_delay); 505} 506 507/* 508 * satfinddelay - calculate the one-way delay time between a ground station 509 * and a satellite 510 */ 511static void 512satfinddelay( 513 double lat1, 514 double long1, 515 double lat2, 516 double long2, 517 double *delay 518 ) 519{ 520 double dg; /* great circle distance */ 521 522 dg = greatcircle(lat1, long1, lat2, long2); 523 524 *delay = satpropdelay(dg); 525} 526 527/* 528 * satpropdelay - calculate the one-way delay time between a ground station 529 * and a satellite 530 */ 531static double 532satpropdelay( 533 double dg 534 ) 535{ 536 double k1, k2, dist; 537 double theta; 538 double td; 539 540 theta = dg / (EARTHRADIUS); 541 k1 = EARTHRADIUS * sin(theta); 542 k2 = SATHEIGHT - (EARTHRADIUS * cos(theta)); 543 if (debug >= 2) 544 printf("Theta %g k1 %g k2 %g\n", theta, k1, k2); 545 dist = sqrt(k1*k1 + k2*k2); 546 td = dist / LIGHTSPEED; 547 if (debug >= 2) 548 printf("propdelay dist %g height %g time %g\n", dg, dist, td); 549 return td; 550} 551