propdelay.c revision 280849
118334Speter/* propdelay.c,v 3.1 1993/07/06 01:05:24 jbj Exp 218334Speter * propdelay - compute propagation delays 318334Speter * 490075Sobrien * cc -o propdelay propdelay.c -lm 590075Sobrien * 618334Speter * "Time and Frequency Users' Manual", NBS Technical Note 695 (1977). 790075Sobrien */ 818334Speter 990075Sobrien/* 1090075Sobrien * This can be used to get a rough idea of the HF propagation delay 1190075Sobrien * between two points (usually between you and the radio station). 1290075Sobrien * The usage is 1318334Speter * 1490075Sobrien * propdelay latitudeA longitudeA latitudeB longitudeB 1590075Sobrien * 1690075Sobrien * where points A and B are the locations in question. You obviously 1790075Sobrien * need to know the latitude and longitude of each of the places. 1818334Speter * The program expects the latitude to be preceded by an 'n' or 's' 1918334Speter * and the longitude to be preceded by an 'e' or 'w'. It understands 2090075Sobrien * either decimal degrees or degrees:minutes:seconds. Thus to compute 2190075Sobrien * the delay between the WWVH (21:59:26N, 159:46:00W) and WWV (40:40:49N, 2290075Sobrien * 105:02:27W) you could use: 2318334Speter * 2418334Speter * propdelay n21:59:26 w159:46 n40:40:49 w105:02:27 2518334Speter * 2618334Speter * By default it prints out a summer (F2 average virtual height 350 km) and 2718334Speter * winter (F2 average virtual height 250 km) number. The results will be 2818334Speter * quite approximate but are about as good as you can do with HF time anyway. 2918334Speter * You might pick a number between the values to use, or use the summer 3018334Speter * value in the summer and switch to the winter value when the static 3118334Speter * above 10 MHz starts to drop off in the fall. You can also use the 3218334Speter * -h switch if you want to specify your own virtual height. 3318334Speter * 3418334Speter * You can also do a 3518334Speter * 3618334Speter * propdelay -W n45:17:47 w75:45:22 3718334Speter * 3818334Speter * to find the propagation delays to WWV and WWVH (from CHU in this 3918334Speter * case), a 4018334Speter * 4118334Speter * propdelay -C n40:40:49 w105:02:27 4218334Speter * 4318334Speter * to find the delays to CHU, and a 4418334Speter * 4518334Speter * propdelay -G n52:03:17 w98:34:18 4618334Speter * 4718334Speter * to find delays to GOES via each of the three satellites. 4818334Speter */ 4918334Speter 5018334Speter#ifdef HAVE_CONFIG_H 5118334Speter# include <config.h> 5218334Speter#endif 5318334Speter#include <stdio.h> 5418334Speter#include <string.h> 5518334Speter 5618334Speter#include "ntp_stdlib.h" 5718334Speter 5850397Sobrienextern double sin (double); 5990075Sobrienextern double cos (double); 6018334Speterextern double acos (double); 6118334Speterextern double tan (double); 6218334Speterextern double atan (double); 6318334Speterextern double sqrt (double); 6418334Speter 6518334Speter#define STREQ(a, b) (*(a) == *(b) && strcmp((a), (b)) == 0) 6618334Speter 6718334Speter/* 6818334Speter * Program constants 6918334Speter */ 7018334Speter#define EARTHRADIUS (6370.0) /* raduis of earth (km) */ 7118334Speter#define LIGHTSPEED (299800.0) /* speed of light, km/s */ 7218334Speter#define PI (3.1415926536) 7318334Speter#define RADPERDEG (PI/180.0) /* radians per degree */ 7418334Speter#define MILE (1.609344) /* km in a mile */ 7518334Speter 7690075Sobrien#define SUMMERHEIGHT (350.0) /* summer height in km */ 7790075Sobrien#define WINTERHEIGHT (250.0) /* winter height in km */ 7890075Sobrien 7990075Sobrien#define SATHEIGHT (6.6110 * 6378.0) /* geosync satellite height in km 8090075Sobrien from centre of earth */ 8118334Speter 8218334Speter#define WWVLAT "n40:40:49" 8318334Speter#define WWVLONG "w105:02:27" 8418334Speter 8518334Speter#define WWVHLAT "n21:59:26" 8618334Speter#define WWVHLONG "w159:46:00" 8718334Speter 8818334Speter#define CHULAT "n45:17:47" 8918334Speter#define CHULONG "w75:45:22" 9018334Speter 9118334Speter#define GOES_UP_LAT "n37:52:00" 9218334Speter#define GOES_UP_LONG "w75:27:00" 9318334Speter#define GOES_EAST_LONG "w75:00:00" 9418334Speter#define GOES_STBY_LONG "w105:00:00" 9518334Speter#define GOES_WEST_LONG "w135:00:00" 9618334Speter#define GOES_SAT_LAT "n00:00:00" 9718334Speter 9818334Speterchar *wwvlat = WWVLAT; 9918334Speterchar *wwvlong = WWVLONG; 10018334Speter 10118334Speterchar *wwvhlat = WWVHLAT; 10218334Speterchar *wwvhlong = WWVHLONG; 10318334Speter 10418334Speterchar *chulat = CHULAT; 10518334Speterchar *chulong = CHULONG; 10618334Speter 10718334Speterchar *goes_up_lat = GOES_UP_LAT; 10818334Speterchar *goes_up_long = GOES_UP_LONG; 10918334Speterchar *goes_east_long = GOES_EAST_LONG; 11018334Speterchar *goes_stby_long = GOES_STBY_LONG; 11118334Speterchar *goes_west_long = GOES_WEST_LONG; 11218334Speterchar *goes_sat_lat = GOES_SAT_LAT; 11350397Sobrien 11418334Speterint hflag = 0; 11518334Speterint Wflag = 0; 11618334Speterint Cflag = 0; 11718334Speterint Gflag = 0; 11818334Speterint height; 11918334Speter 12018334Speterchar *progname; 12118334Speter 12218334Speterstatic void doit (double, double, double, double, double, char *); 12318334Speterstatic double latlong (char *, int); 12418334Speterstatic double greatcircle (double, double, double, double); 12518334Speterstatic double waveangle (double, double, int); 12618334Speterstatic double propdelay (double, double, int); 12718334Speterstatic int finddelay (double, double, double, double, double, double *); 12818334Speterstatic void satdoit (double, double, double, double, double, double, char *); 12918334Speterstatic void satfinddelay (double, double, double, double, double *); 13018334Speterstatic double satpropdelay (double); 13118334Speter 13218334Speter/* 13318334Speter * main - parse arguments and handle options 13418334Speter */ 13518334Speterint 13618334Spetermain( 13718334Speter int argc, 13818334Speter char *argv[] 13918334Speter ) 14018334Speter{ 14118334Speter int c; 14218334Speter int errflg = 0; 14318334Speter double lat1, long1; 14418334Speter double lat2, long2; 14518334Speter double lat3, long3; 14618334Speter 14718334Speter init_lib(); 14818334Speter 14918334Speter progname = argv[0]; 15018334Speter while ((c = ntp_getopt(argc, argv, "dh:CWG")) != EOF) 15118334Speter switch (c) { 15218334Speter case 'd': 15318334Speter ++debug; 15418334Speter break; 15518334Speter case 'h': 15618334Speter hflag++; 15718334Speter height = atof(ntp_optarg); 15818334Speter if (height <= 0.0) { 15918334Speter (void) fprintf(stderr, "height %s unlikely\n", 16018334Speter ntp_optarg); 16118334Speter errflg++; 16218334Speter } 16318334Speter break; 16418334Speter case 'C': 16518334Speter Cflag++; 16618334Speter break; 16718334Speter case 'W': 16818334Speter Wflag++; 16918334Speter break; 17018334Speter case 'G': 17150397Sobrien Gflag++; 17250397Sobrien break; 17390075Sobrien default: 17450397Sobrien errflg++; 17550397Sobrien break; 17618334Speter } 17718334Speter if (errflg || (!(Cflag || Wflag || Gflag) && ntp_optind+4 != argc) || 17818334Speter ((Cflag || Wflag || Gflag) && ntp_optind+2 != argc)) { 17918334Speter (void) fprintf(stderr, 18018334Speter "usage: %s [-d] [-h height] lat1 long1 lat2 long2\n", 18118334Speter progname); 18218334Speter (void) fprintf(stderr," - or -\n"); 18318334Speter (void) fprintf(stderr, 18418334Speter "usage: %s -CWG [-d] lat long\n", 18518334Speter progname); 18618334Speter exit(2); 18718334Speter } 18818334Speter 18918334Speter 19018334Speter if (!(Cflag || Wflag || Gflag)) { 19118334Speter lat1 = latlong(argv[ntp_optind], 1); 19290075Sobrien long1 = latlong(argv[ntp_optind + 1], 0); 19318334Speter lat2 = latlong(argv[ntp_optind + 2], 1); 19418334Speter long2 = latlong(argv[ntp_optind + 3], 0); 19518334Speter if (hflag) { 19618334Speter doit(lat1, long1, lat2, long2, height, ""); 19718334Speter } else { 19818334Speter doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 19918334Speter "summer propagation, "); 20090075Sobrien doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 20118334Speter "winter propagation, "); 20218334Speter } 20318334Speter } else if (Wflag) { 20418334Speter /* 20518334Speter * Compute delay from WWV 20618334Speter */ 20718334Speter lat1 = latlong(argv[ntp_optind], 1); 20890075Sobrien long1 = latlong(argv[ntp_optind + 1], 0); 20990075Sobrien lat2 = latlong(wwvlat, 1); 21090075Sobrien long2 = latlong(wwvlong, 0); 21118334Speter if (hflag) { 21218334Speter doit(lat1, long1, lat2, long2, height, "WWV "); 21390075Sobrien } else { 21490075Sobrien doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 21590075Sobrien "WWV summer propagation, "); 21690075Sobrien doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 21790075Sobrien "WWV winter propagation, "); 21890075Sobrien } 21990075Sobrien 22090075Sobrien /* 22190075Sobrien * Compute delay from WWVH 22290075Sobrien */ 22390075Sobrien lat2 = latlong(wwvhlat, 1); 22490075Sobrien long2 = latlong(wwvhlong, 0); 22590075Sobrien if (hflag) { 22690075Sobrien doit(lat1, long1, lat2, long2, height, "WWVH "); 22790075Sobrien } else { 22890075Sobrien doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 22990075Sobrien "WWVH summer propagation, "); 23090075Sobrien doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 23190075Sobrien "WWVH winter propagation, "); 23290075Sobrien } 23390075Sobrien } else if (Cflag) { 23490075Sobrien lat1 = latlong(argv[ntp_optind], 1); 23590075Sobrien long1 = latlong(argv[ntp_optind + 1], 0); 23690075Sobrien lat2 = latlong(chulat, 1); 23790075Sobrien long2 = latlong(chulong, 0); 23890075Sobrien if (hflag) { 23990075Sobrien doit(lat1, long1, lat2, long2, height, "CHU "); 24090075Sobrien } else { 24190075Sobrien doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 24290075Sobrien "CHU summer propagation, "); 24390075Sobrien doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 24490075Sobrien "CHU winter propagation, "); 24518334Speter } 24618334Speter } else if (Gflag) { 24718334Speter lat1 = latlong(goes_up_lat, 1); 24818334Speter long1 = latlong(goes_up_long, 0); 24918334Speter lat3 = latlong(argv[ntp_optind], 1); 25018334Speter long3 = latlong(argv[ntp_optind + 1], 0); 25118334Speter 25290075Sobrien lat2 = latlong(goes_sat_lat, 1); 25318334Speter 25418334Speter long2 = latlong(goes_west_long, 0); 25518334Speter satdoit(lat1, long1, lat2, long2, lat3, long3, 25618334Speter "GOES Delay via WEST"); 25790075Sobrien 25890075Sobrien long2 = latlong(goes_stby_long, 0); 25990075Sobrien satdoit(lat1, long1, lat2, long2, lat3, long3, 26018334Speter "GOES Delay via STBY"); 26118334Speter 26218334Speter long2 = latlong(goes_east_long, 0); 26318334Speter satdoit(lat1, long1, lat2, long2, lat3, long3, 26418334Speter "GOES Delay via EAST"); 26518334Speter 26618334Speter } 26718334Speter exit(0); 26818334Speter} 26918334Speter 27018334Speter 27118334Speter/* 27218334Speter * doit - compute a delay and print it 27318334Speter */ 27418334Speterstatic void 27518334Speterdoit( 27618334Speter double lat1, 27718334Speter double long1, 27818334Speter double lat2, 27918334Speter double long2, 28018334Speter double h, 28118334Speter char *str 28218334Speter ) 28318334Speter{ 28418334Speter int hops; 28518334Speter double delay; 28618334Speter 28718334Speter hops = finddelay(lat1, long1, lat2, long2, h, &delay); 28818334Speter printf("%sheight %g km, hops %d, delay %g seconds\n", 28918334Speter str, h, hops, delay); 29018334Speter} 29118334Speter 29218334Speter 29318334Speter/* 29418334Speter * latlong - decode a latitude/longitude value 29518334Speter */ 29618334Speterstatic double 29718334Speterlatlong( 29818334Speter char *str, 29918334Speter int islat 30018334Speter ) 30118334Speter{ 30218334Speter register char *cp; 30318334Speter register char *bp; 30418334Speter double arg; 30518334Speter double divby; 30618334Speter int isneg; 30750397Sobrien char buf[32]; 30818334Speter char *colon; 30918334Speter 31018334Speter if (islat) { 31118334Speter /* 31218334Speter * Must be north or south 31318334Speter */ 31418334Speter if (*str == 'N' || *str == 'n') 31518334Speter isneg = 0; 31618334Speter else if (*str == 'S' || *str == 's') 31718334Speter isneg = 1; 31818334Speter else 31918334Speter isneg = -1; 32018334Speter } else { 32118334Speter /* 32290075Sobrien * East is positive, west is negative 32390075Sobrien */ 32490075Sobrien if (*str == 'E' || *str == 'e') 32590075Sobrien isneg = 0; 32690075Sobrien else if (*str == 'W' || *str == 'w') 32790075Sobrien isneg = 1; 32890075Sobrien else 32990075Sobrien isneg = -1; 33090075Sobrien } 33190075Sobrien 33290075Sobrien if (isneg >= 0) 33390075Sobrien str++; 33490075Sobrien 33518334Speter colon = strchr(str, ':'); 33618334Speter if (colon != NULL) { 33718334Speter /* 33818334Speter * in hhh:mm:ss form 33918334Speter */ 34018334Speter cp = str; 34190075Sobrien bp = buf; 34218334Speter while (cp < colon) 34318334Speter *bp++ = *cp++; 34418334Speter *bp = '\0'; 34518334Speter cp++; 34618334Speter arg = atof(buf); 34718334Speter divby = 60.0; 34818334Speter colon = strchr(cp, ':'); 34918334Speter if (colon != NULL) { 35018334Speter bp = buf; 35118334Speter while (cp < colon) 35218334Speter *bp++ = *cp++; 35318334Speter *bp = '\0'; 35418334Speter cp++; 35518334Speter arg += atof(buf) / divby; 35618334Speter divby = 3600.0; 35718334Speter } 35818334Speter if (*cp != '\0') 35918334Speter arg += atof(cp) / divby; 36018334Speter } else { 36118334Speter arg = atof(str); 36218334Speter } 36318334Speter 36418334Speter if (isneg == 1) 36518334Speter arg = -arg; 36618334Speter 36718334Speter if (debug > 2) 36818334Speter (void) printf("latitude/longitude %s = %g\n", str, arg); 36918334Speter 37018334Speter return arg; 37118334Speter} 37218334Speter 37318334Speter 37418334Speter/* 37518334Speter * greatcircle - compute the great circle distance in kilometers 37618334Speter */ 37718334Speterstatic double 37818334Spetergreatcircle( 37918334Speter double lat1, 38018334Speter double long1, 38118334Speter double lat2, 38218334Speter double long2 38318334Speter ) 38418334Speter{ 38518334Speter double dg; 38618334Speter double l1r, l2r; 38718334Speter 38818334Speter l1r = lat1 * RADPERDEG; 38918334Speter l2r = lat2 * RADPERDEG; 39018334Speter dg = EARTHRADIUS * acos( 39118334Speter (cos(l1r) * cos(l2r) * cos((long2-long1)*RADPERDEG)) 39218334Speter + (sin(l1r) * sin(l2r))); 39318334Speter if (debug >= 2) 39418334Speter printf( 39518334Speter "greatcircle lat1 %g long1 %g lat2 %g long2 %g dist %g\n", 39618334Speter lat1, long1, lat2, long2, dg); 39718334Speter return dg; 39818334Speter} 39918334Speter 40018334Speter 40118334Speter/* 40218334Speter * waveangle - compute the wave angle for the given distance, virtual 40318334Speter * height and number of hops. 40418334Speter */ 40518334Speterstatic double 40618334Speterwaveangle( 40718334Speter double dg, 40818334Speter double h, 40918334Speter int n 41018334Speter ) 41118334Speter{ 41218334Speter double theta; 41318334Speter double delta; 41418334Speter 41518334Speter theta = dg / (EARTHRADIUS * (double)(2 * n)); 41618334Speter delta = atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2)) - theta; 41718334Speter if (debug >= 2) 41818334Speter printf("waveangle dist %g height %g hops %d angle %g\n", 41918334Speter dg, h, n, delta / RADPERDEG); 42018334Speter return delta; 42118334Speter} 42218334Speter 42318334Speter 42418334Speter/* 42590075Sobrien * propdelay - compute the propagation delay 42690075Sobrien */ 42790075Sobrienstatic double 42890075Sobrienpropdelay( 42990075Sobrien double dg, 43090075Sobrien double h, 43118334Speter int n 43218334Speter ) 43390075Sobrien{ 43490075Sobrien double phi; 43590075Sobrien double theta; 43690075Sobrien double td; 43790075Sobrien 43818334Speter theta = dg / (EARTHRADIUS * (double)(2 * n)); 43918334Speter phi = (PI/2.0) - atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2)); 44018334Speter td = dg / (LIGHTSPEED * sin(phi)); 44118334Speter if (debug >= 2) 44218334Speter printf("propdelay dist %g height %g hops %d time %g\n", 44390075Sobrien dg, h, n, td); 44490075Sobrien return td; 44590075Sobrien} 44690075Sobrien 44790075Sobrien 44890075Sobrien/* 44990075Sobrien * finddelay - find the propagation delay 45090075Sobrien */ 45190075Sobrienstatic int 45290075Sobrienfinddelay( 45318334Speter double lat1, 45418334Speter double long1, 45518334Speter double lat2, 45618334Speter double long2, 45718334Speter double h, 45818334Speter double *delay 45918334Speter ) 46018334Speter{ 46118334Speter double dg; /* great circle distance */ 46218334Speter double delta; /* wave angle */ 46318334Speter int n; /* number of hops */ 46418334Speter 46518334Speter dg = greatcircle(lat1, long1, lat2, long2); 46618334Speter if (debug) 46718334Speter printf("great circle distance %g km %g miles\n", dg, dg/MILE); 46818334Speter 46918334Speter n = 1; 47018334Speter while ((delta = waveangle(dg, h, n)) < 0.0) { 47118334Speter if (debug) 47218334Speter printf("tried %d hop%s, no good\n", n, n>1?"s":""); 47318334Speter n++; 47418334Speter } 47518334Speter if (debug) 47618334Speter printf("%d hop%s okay, wave angle is %g\n", n, n>1?"s":"", 47718334Speter delta / RADPERDEG); 47818334Speter 47918334Speter *delay = propdelay(dg, h, n); 48018334Speter return n; 48118334Speter} 48218334Speter 48318334Speter/* 48418334Speter * satdoit - compute a delay and print it 48518334Speter */ 48618334Speterstatic void 48718334Spetersatdoit( 48818334Speter double lat1, 48918334Speter double long1, 49018334Speter double lat2, 49118334Speter double long2, 49218334Speter double lat3, 49318334Speter double long3, 49418334Speter char *str 49518334Speter ) 49618334Speter{ 49718334Speter double up_delay,down_delay; 49850397Sobrien 49950397Sobrien satfinddelay(lat1, long1, lat2, long2, &up_delay); 50050397Sobrien satfinddelay(lat3, long3, lat2, long2, &down_delay); 50150397Sobrien 50250397Sobrien printf("%s, delay %g seconds\n", str, up_delay + down_delay); 50350397Sobrien} 50450397Sobrien 50550397Sobrien/* 50650397Sobrien * satfinddelay - calculate the one-way delay time between a ground station 50750397Sobrien * and a satellite 50850397Sobrien */ 50950397Sobrienstatic void 51050397Sobriensatfinddelay( 51150397Sobrien double lat1, 51250397Sobrien double long1, 51318334Speter double lat2, 51450397Sobrien double long2, 51550397Sobrien double *delay 51650397Sobrien ) 51750397Sobrien{ 51890075Sobrien double dg; /* great circle distance */ 51950397Sobrien 52090075Sobrien dg = greatcircle(lat1, long1, lat2, long2); 52190075Sobrien 52290075Sobrien *delay = satpropdelay(dg); 52390075Sobrien} 52490075Sobrien 52590075Sobrien/* 52690075Sobrien * satpropdelay - calculate the one-way delay time between a ground station 52790075Sobrien * and a satellite 52890075Sobrien */ 52990075Sobrienstatic double 53090075Sobriensatpropdelay( 53118334Speter double dg 53218334Speter ) 53318334Speter{ 53418334Speter double k1, k2, dist; 53518334Speter double theta; 53618334Speter double td; 53718334Speter 53818334Speter theta = dg / (EARTHRADIUS); 53918334Speter k1 = EARTHRADIUS * sin(theta); 54018334Speter k2 = SATHEIGHT - (EARTHRADIUS * cos(theta)); 54118334Speter if (debug >= 2) 54218334Speter printf("Theta %g k1 %g k2 %g\n", theta, k1, k2); 54318334Speter dist = sqrt(k1*k1 + k2*k2); 54418334Speter td = dist / LIGHTSPEED; 54518334Speter if (debug >= 2) 54618334Speter printf("propdelay dist %g height %g time %g\n", dg, dist, td); 54718334Speter return td; 54818334Speter} 54918334Speter