Home | History | Annotate | Line # | Download | only in clockstuff
propdelay.c revision 1.1.1.2 
      1 /*	$NetBSD: propdelay.c,v 1.1.1.2 2013/12/27 23:30:33 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 
     60 extern	double	sin	(double);
     61 extern	double	cos	(double);
     62 extern	double	acos	(double);
     63 extern	double	tan	(double);
     64 extern	double	atan	(double);
     65 extern	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 
    100 char *wwvlat = WWVLAT;
    101 char *wwvlong = WWVLONG;
    102 
    103 char *wwvhlat = WWVHLAT;
    104 char *wwvhlong = WWVHLONG;
    105 
    106 char *chulat = CHULAT;
    107 char *chulong = CHULONG;
    108 
    109 char *goes_up_lat = GOES_UP_LAT;
    110 char *goes_up_long = GOES_UP_LONG;
    111 char *goes_east_long = GOES_EAST_LONG;
    112 char *goes_stby_long = GOES_STBY_LONG;
    113 char *goes_west_long = GOES_WEST_LONG;
    114 char *goes_sat_lat = GOES_SAT_LAT;
    115 
    116 int hflag = 0;
    117 int Wflag = 0;
    118 int Cflag = 0;
    119 int Gflag = 0;
    120 int height;
    121 
    122 char *progname;
    123 
    124 static	void	doit		(double, double, double, double, double, char *);
    125 static	double	latlong		(char *, int);
    126 static	double	greatcircle	(double, double, double, double);
    127 static	double	waveangle	(double, double, int);
    128 static	double	propdelay	(double, double, int);
    129 static	int	finddelay	(double, double, double, double, double, double *);
    130 static	void	satdoit		(double, double, double, double, double, double, char *);
    131 static	void	satfinddelay	(double, double, double, double, double *);
    132 static	double	satpropdelay	(double);
    133 
    134 /*
    135  * main - parse arguments and handle options
    136  */
    137 int
    138 main(
    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  */
    276 static void
    277 doit(
    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  */
    298 static double
    299 latlong(
    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  */
    379 static double
    380 greatcircle(
    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  */
    407 static double
    408 waveangle(
    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  */
    429 static double
    430 propdelay(
    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  */
    453 static int
    454 finddelay(
    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  */
    488 static void
    489 satdoit(
    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  */
    511 static void
    512 satfinddelay(
    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  */
    531 static double
    532 satpropdelay(
    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