File:  [ELWIX - Embedded LightWeight unIX -] / embedaddon / ntp / clockstuff / propdelay.c
Revision 1.1.1.1 (vendor branch): download - view: text, annotated - select for diffs - revision graph
Tue May 29 12:08:38 2012 UTC (12 years, 7 months ago) by misho
Branches: ntp, MAIN
CVS tags: v4_2_6p5p0, v4_2_6p5, HEAD
ntp 4.2.6p5

    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>