Annotation of embedaddon/ntp/clockstuff/propdelay.c, revision 1.1
1.1 ! misho 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: volatile 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 divby;
! 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: divby = 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) / divby;
! 352: divby = 3600.0;
! 353: }
! 354: if (*cp != '\0')
! 355: arg += atof(cp) / divby;
! 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: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to propdelay.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [ELWIX - Embedded LightWeight unIX -] / embedaddon / ntp / clockstuff