154359Sroberto/* propdelay.c,v 3.1 1993/07/06 01:05:24 jbj Exp 254359Sroberto * propdelay - compute propagation delays 354359Sroberto * 454359Sroberto * cc -o propdelay propdelay.c -lm 554359Sroberto * 654359Sroberto * "Time and Frequency Users' Manual", NBS Technical Note 695 (1977). 754359Sroberto */ 854359Sroberto 954359Sroberto/* 1054359Sroberto * This can be used to get a rough idea of the HF propagation delay 1154359Sroberto * between two points (usually between you and the radio station). 1254359Sroberto * The usage is 1354359Sroberto * 1454359Sroberto * propdelay latitudeA longitudeA latitudeB longitudeB 1554359Sroberto * 1654359Sroberto * where points A and B are the locations in question. You obviously 1754359Sroberto * need to know the latitude and longitude of each of the places. 1854359Sroberto * The program expects the latitude to be preceded by an 'n' or 's' 1954359Sroberto * and the longitude to be preceded by an 'e' or 'w'. It understands 2054359Sroberto * either decimal degrees or degrees:minutes:seconds. Thus to compute 2154359Sroberto * the delay between the WWVH (21:59:26N, 159:46:00W) and WWV (40:40:49N, 2254359Sroberto * 105:02:27W) you could use: 2354359Sroberto * 2454359Sroberto * propdelay n21:59:26 w159:46 n40:40:49 w105:02:27 2554359Sroberto * 2654359Sroberto * By default it prints out a summer (F2 average virtual height 350 km) and 2754359Sroberto * winter (F2 average virtual height 250 km) number. The results will be 2854359Sroberto * quite approximate but are about as good as you can do with HF time anyway. 2954359Sroberto * You might pick a number between the values to use, or use the summer 3054359Sroberto * value in the summer and switch to the winter value when the static 3154359Sroberto * above 10 MHz starts to drop off in the fall. You can also use the 3254359Sroberto * -h switch if you want to specify your own virtual height. 3354359Sroberto * 3454359Sroberto * You can also do a 3554359Sroberto * 3654359Sroberto * propdelay -W n45:17:47 w75:45:22 3754359Sroberto * 3854359Sroberto * to find the propagation delays to WWV and WWVH (from CHU in this 3954359Sroberto * case), a 4054359Sroberto * 4154359Sroberto * propdelay -C n40:40:49 w105:02:27 4254359Sroberto * 4354359Sroberto * to find the delays to CHU, and a 4454359Sroberto * 4554359Sroberto * propdelay -G n52:03:17 w98:34:18 4654359Sroberto * 4754359Sroberto * to find delays to GOES via each of the three satellites. 4854359Sroberto */ 4954359Sroberto 50285612Sdelphij#ifdef HAVE_CONFIG_H 51285612Sdelphij# include <config.h> 52285612Sdelphij#endif 5354359Sroberto#include <stdio.h> 5454359Sroberto#include <string.h> 5554359Sroberto 5654359Sroberto#include "ntp_stdlib.h" 5754359Sroberto 5854359Srobertoextern double sin (double); 5954359Srobertoextern double cos (double); 6054359Srobertoextern double acos (double); 6154359Srobertoextern double tan (double); 6254359Srobertoextern double atan (double); 6354359Srobertoextern double sqrt (double); 6454359Sroberto 6554359Sroberto#define STREQ(a, b) (*(a) == *(b) && strcmp((a), (b)) == 0) 6654359Sroberto 6754359Sroberto/* 6854359Sroberto * Program constants 6954359Sroberto */ 7054359Sroberto#define EARTHRADIUS (6370.0) /* raduis of earth (km) */ 7154359Sroberto#define LIGHTSPEED (299800.0) /* speed of light, km/s */ 7254359Sroberto#define PI (3.1415926536) 7354359Sroberto#define RADPERDEG (PI/180.0) /* radians per degree */ 7454359Sroberto#define MILE (1.609344) /* km in a mile */ 7554359Sroberto 7654359Sroberto#define SUMMERHEIGHT (350.0) /* summer height in km */ 7754359Sroberto#define WINTERHEIGHT (250.0) /* winter height in km */ 7854359Sroberto 7954359Sroberto#define SATHEIGHT (6.6110 * 6378.0) /* geosync satellite height in km 8054359Sroberto from centre of earth */ 8154359Sroberto 8254359Sroberto#define WWVLAT "n40:40:49" 8354359Sroberto#define WWVLONG "w105:02:27" 8454359Sroberto 8554359Sroberto#define WWVHLAT "n21:59:26" 8654359Sroberto#define WWVHLONG "w159:46:00" 8754359Sroberto 8854359Sroberto#define CHULAT "n45:17:47" 8954359Sroberto#define CHULONG "w75:45:22" 9054359Sroberto 9154359Sroberto#define GOES_UP_LAT "n37:52:00" 9254359Sroberto#define GOES_UP_LONG "w75:27:00" 9354359Sroberto#define GOES_EAST_LONG "w75:00:00" 9454359Sroberto#define GOES_STBY_LONG "w105:00:00" 9554359Sroberto#define GOES_WEST_LONG "w135:00:00" 9654359Sroberto#define GOES_SAT_LAT "n00:00:00" 9754359Sroberto 9854359Srobertochar *wwvlat = WWVLAT; 9954359Srobertochar *wwvlong = WWVLONG; 10054359Sroberto 10154359Srobertochar *wwvhlat = WWVHLAT; 10254359Srobertochar *wwvhlong = WWVHLONG; 10354359Sroberto 10454359Srobertochar *chulat = CHULAT; 10554359Srobertochar *chulong = CHULONG; 10654359Sroberto 10754359Srobertochar *goes_up_lat = GOES_UP_LAT; 10854359Srobertochar *goes_up_long = GOES_UP_LONG; 10954359Srobertochar *goes_east_long = GOES_EAST_LONG; 11054359Srobertochar *goes_stby_long = GOES_STBY_LONG; 11154359Srobertochar *goes_west_long = GOES_WEST_LONG; 11254359Srobertochar *goes_sat_lat = GOES_SAT_LAT; 11354359Sroberto 11454359Srobertoint hflag = 0; 11554359Srobertoint Wflag = 0; 11654359Srobertoint Cflag = 0; 11754359Srobertoint Gflag = 0; 11854359Srobertoint height; 11954359Sroberto 120289997Sglebiuschar const *progname; 12154359Sroberto 12254359Srobertostatic void doit (double, double, double, double, double, char *); 12354359Srobertostatic double latlong (char *, int); 12454359Srobertostatic double greatcircle (double, double, double, double); 12554359Srobertostatic double waveangle (double, double, int); 12654359Srobertostatic double propdelay (double, double, int); 12754359Srobertostatic int finddelay (double, double, double, double, double, double *); 12854359Srobertostatic void satdoit (double, double, double, double, double, double, char *); 12954359Srobertostatic void satfinddelay (double, double, double, double, double *); 13054359Srobertostatic double satpropdelay (double); 13154359Sroberto 13254359Sroberto/* 13354359Sroberto * main - parse arguments and handle options 13454359Sroberto */ 13554359Srobertoint 13654359Srobertomain( 13754359Sroberto int argc, 13854359Sroberto char *argv[] 13954359Sroberto ) 14054359Sroberto{ 14154359Sroberto int c; 14254359Sroberto int errflg = 0; 14354359Sroberto double lat1, long1; 14454359Sroberto double lat2, long2; 14554359Sroberto double lat3, long3; 14654359Sroberto 147285612Sdelphij init_lib(); 148285612Sdelphij 14954359Sroberto progname = argv[0]; 15054359Sroberto while ((c = ntp_getopt(argc, argv, "dh:CWG")) != EOF) 15154359Sroberto switch (c) { 15254359Sroberto case 'd': 15354359Sroberto ++debug; 15454359Sroberto break; 15554359Sroberto case 'h': 15654359Sroberto hflag++; 15754359Sroberto height = atof(ntp_optarg); 15854359Sroberto if (height <= 0.0) { 15954359Sroberto (void) fprintf(stderr, "height %s unlikely\n", 16054359Sroberto ntp_optarg); 16154359Sroberto errflg++; 16254359Sroberto } 16354359Sroberto break; 16454359Sroberto case 'C': 16554359Sroberto Cflag++; 16654359Sroberto break; 16754359Sroberto case 'W': 16854359Sroberto Wflag++; 16954359Sroberto break; 17054359Sroberto case 'G': 17154359Sroberto Gflag++; 17254359Sroberto break; 17354359Sroberto default: 17454359Sroberto errflg++; 17554359Sroberto break; 17654359Sroberto } 17754359Sroberto if (errflg || (!(Cflag || Wflag || Gflag) && ntp_optind+4 != argc) || 17854359Sroberto ((Cflag || Wflag || Gflag) && ntp_optind+2 != argc)) { 17954359Sroberto (void) fprintf(stderr, 18054359Sroberto "usage: %s [-d] [-h height] lat1 long1 lat2 long2\n", 18154359Sroberto progname); 18254359Sroberto (void) fprintf(stderr," - or -\n"); 18354359Sroberto (void) fprintf(stderr, 18454359Sroberto "usage: %s -CWG [-d] lat long\n", 18554359Sroberto progname); 18654359Sroberto exit(2); 18754359Sroberto } 18854359Sroberto 18954359Sroberto 19054359Sroberto if (!(Cflag || Wflag || Gflag)) { 19154359Sroberto lat1 = latlong(argv[ntp_optind], 1); 19254359Sroberto long1 = latlong(argv[ntp_optind + 1], 0); 19354359Sroberto lat2 = latlong(argv[ntp_optind + 2], 1); 19454359Sroberto long2 = latlong(argv[ntp_optind + 3], 0); 19554359Sroberto if (hflag) { 19654359Sroberto doit(lat1, long1, lat2, long2, height, ""); 19754359Sroberto } else { 19854359Sroberto doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 19954359Sroberto "summer propagation, "); 20054359Sroberto doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 20154359Sroberto "winter propagation, "); 20254359Sroberto } 20354359Sroberto } else if (Wflag) { 20454359Sroberto /* 20554359Sroberto * Compute delay from WWV 20654359Sroberto */ 20754359Sroberto lat1 = latlong(argv[ntp_optind], 1); 20854359Sroberto long1 = latlong(argv[ntp_optind + 1], 0); 20954359Sroberto lat2 = latlong(wwvlat, 1); 21054359Sroberto long2 = latlong(wwvlong, 0); 21154359Sroberto if (hflag) { 21254359Sroberto doit(lat1, long1, lat2, long2, height, "WWV "); 21354359Sroberto } else { 21454359Sroberto doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 21554359Sroberto "WWV summer propagation, "); 21654359Sroberto doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 21754359Sroberto "WWV winter propagation, "); 21854359Sroberto } 21954359Sroberto 22054359Sroberto /* 22154359Sroberto * Compute delay from WWVH 22254359Sroberto */ 22354359Sroberto lat2 = latlong(wwvhlat, 1); 22454359Sroberto long2 = latlong(wwvhlong, 0); 22554359Sroberto if (hflag) { 22654359Sroberto doit(lat1, long1, lat2, long2, height, "WWVH "); 22754359Sroberto } else { 22854359Sroberto doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 22954359Sroberto "WWVH summer propagation, "); 23054359Sroberto doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 23154359Sroberto "WWVH winter propagation, "); 23254359Sroberto } 23354359Sroberto } else if (Cflag) { 23454359Sroberto lat1 = latlong(argv[ntp_optind], 1); 23554359Sroberto long1 = latlong(argv[ntp_optind + 1], 0); 23654359Sroberto lat2 = latlong(chulat, 1); 23754359Sroberto long2 = latlong(chulong, 0); 23854359Sroberto if (hflag) { 23954359Sroberto doit(lat1, long1, lat2, long2, height, "CHU "); 24054359Sroberto } else { 24154359Sroberto doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 24254359Sroberto "CHU summer propagation, "); 24354359Sroberto doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 24454359Sroberto "CHU winter propagation, "); 24554359Sroberto } 24654359Sroberto } else if (Gflag) { 24754359Sroberto lat1 = latlong(goes_up_lat, 1); 24854359Sroberto long1 = latlong(goes_up_long, 0); 24954359Sroberto lat3 = latlong(argv[ntp_optind], 1); 25054359Sroberto long3 = latlong(argv[ntp_optind + 1], 0); 25154359Sroberto 25254359Sroberto lat2 = latlong(goes_sat_lat, 1); 25354359Sroberto 25454359Sroberto long2 = latlong(goes_west_long, 0); 25554359Sroberto satdoit(lat1, long1, lat2, long2, lat3, long3, 25654359Sroberto "GOES Delay via WEST"); 25754359Sroberto 25854359Sroberto long2 = latlong(goes_stby_long, 0); 25954359Sroberto satdoit(lat1, long1, lat2, long2, lat3, long3, 26054359Sroberto "GOES Delay via STBY"); 26154359Sroberto 26254359Sroberto long2 = latlong(goes_east_long, 0); 26354359Sroberto satdoit(lat1, long1, lat2, long2, lat3, long3, 26454359Sroberto "GOES Delay via EAST"); 26554359Sroberto 26654359Sroberto } 26754359Sroberto exit(0); 26854359Sroberto} 26954359Sroberto 27054359Sroberto 27154359Sroberto/* 27254359Sroberto * doit - compute a delay and print it 27354359Sroberto */ 27454359Srobertostatic void 27554359Srobertodoit( 27654359Sroberto double lat1, 27754359Sroberto double long1, 27854359Sroberto double lat2, 27954359Sroberto double long2, 28054359Sroberto double h, 28154359Sroberto char *str 28254359Sroberto ) 28354359Sroberto{ 28454359Sroberto int hops; 28554359Sroberto double delay; 28654359Sroberto 28754359Sroberto hops = finddelay(lat1, long1, lat2, long2, h, &delay); 28854359Sroberto printf("%sheight %g km, hops %d, delay %g seconds\n", 28954359Sroberto str, h, hops, delay); 29054359Sroberto} 29154359Sroberto 29254359Sroberto 29354359Sroberto/* 29454359Sroberto * latlong - decode a latitude/longitude value 29554359Sroberto */ 29654359Srobertostatic double 29754359Srobertolatlong( 29854359Sroberto char *str, 29954359Sroberto int islat 30054359Sroberto ) 30154359Sroberto{ 30254359Sroberto register char *cp; 30354359Sroberto register char *bp; 30454359Sroberto double arg; 305182007Sroberto double divby; 30654359Sroberto int isneg; 30754359Sroberto char buf[32]; 30854359Sroberto char *colon; 30954359Sroberto 31054359Sroberto if (islat) { 31154359Sroberto /* 31254359Sroberto * Must be north or south 31354359Sroberto */ 31454359Sroberto if (*str == 'N' || *str == 'n') 31554359Sroberto isneg = 0; 31654359Sroberto else if (*str == 'S' || *str == 's') 31754359Sroberto isneg = 1; 31854359Sroberto else 31954359Sroberto isneg = -1; 32054359Sroberto } else { 32154359Sroberto /* 32254359Sroberto * East is positive, west is negative 32354359Sroberto */ 32454359Sroberto if (*str == 'E' || *str == 'e') 32554359Sroberto isneg = 0; 32654359Sroberto else if (*str == 'W' || *str == 'w') 32754359Sroberto isneg = 1; 32854359Sroberto else 32954359Sroberto isneg = -1; 33054359Sroberto } 33154359Sroberto 33254359Sroberto if (isneg >= 0) 33354359Sroberto str++; 33454359Sroberto 33554359Sroberto colon = strchr(str, ':'); 33654359Sroberto if (colon != NULL) { 33754359Sroberto /* 33854359Sroberto * in hhh:mm:ss form 33954359Sroberto */ 34054359Sroberto cp = str; 34154359Sroberto bp = buf; 34254359Sroberto while (cp < colon) 34354359Sroberto *bp++ = *cp++; 34454359Sroberto *bp = '\0'; 34554359Sroberto cp++; 34654359Sroberto arg = atof(buf); 347182007Sroberto divby = 60.0; 34854359Sroberto colon = strchr(cp, ':'); 34954359Sroberto if (colon != NULL) { 35054359Sroberto bp = buf; 35154359Sroberto while (cp < colon) 35254359Sroberto *bp++ = *cp++; 35354359Sroberto *bp = '\0'; 35454359Sroberto cp++; 355182007Sroberto arg += atof(buf) / divby; 356182007Sroberto divby = 3600.0; 35754359Sroberto } 35854359Sroberto if (*cp != '\0') 359182007Sroberto arg += atof(cp) / divby; 36054359Sroberto } else { 36154359Sroberto arg = atof(str); 36254359Sroberto } 36354359Sroberto 36454359Sroberto if (isneg == 1) 36554359Sroberto arg = -arg; 36654359Sroberto 36754359Sroberto if (debug > 2) 36854359Sroberto (void) printf("latitude/longitude %s = %g\n", str, arg); 36954359Sroberto 37054359Sroberto return arg; 37154359Sroberto} 37254359Sroberto 37354359Sroberto 37454359Sroberto/* 37554359Sroberto * greatcircle - compute the great circle distance in kilometers 37654359Sroberto */ 37754359Srobertostatic double 37854359Srobertogreatcircle( 37954359Sroberto double lat1, 38054359Sroberto double long1, 38154359Sroberto double lat2, 38254359Sroberto double long2 38354359Sroberto ) 38454359Sroberto{ 38554359Sroberto double dg; 38654359Sroberto double l1r, l2r; 38754359Sroberto 38854359Sroberto l1r = lat1 * RADPERDEG; 38954359Sroberto l2r = lat2 * RADPERDEG; 39054359Sroberto dg = EARTHRADIUS * acos( 39154359Sroberto (cos(l1r) * cos(l2r) * cos((long2-long1)*RADPERDEG)) 39254359Sroberto + (sin(l1r) * sin(l2r))); 39354359Sroberto if (debug >= 2) 39454359Sroberto printf( 39554359Sroberto "greatcircle lat1 %g long1 %g lat2 %g long2 %g dist %g\n", 39654359Sroberto lat1, long1, lat2, long2, dg); 39754359Sroberto return dg; 39854359Sroberto} 39954359Sroberto 40054359Sroberto 40154359Sroberto/* 40254359Sroberto * waveangle - compute the wave angle for the given distance, virtual 40354359Sroberto * height and number of hops. 40454359Sroberto */ 40554359Srobertostatic double 40654359Srobertowaveangle( 40754359Sroberto double dg, 40854359Sroberto double h, 40954359Sroberto int n 41054359Sroberto ) 41154359Sroberto{ 41254359Sroberto double theta; 41354359Sroberto double delta; 41454359Sroberto 41554359Sroberto theta = dg / (EARTHRADIUS * (double)(2 * n)); 41654359Sroberto delta = atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2)) - theta; 41754359Sroberto if (debug >= 2) 41854359Sroberto printf("waveangle dist %g height %g hops %d angle %g\n", 41954359Sroberto dg, h, n, delta / RADPERDEG); 42054359Sroberto return delta; 42154359Sroberto} 42254359Sroberto 42354359Sroberto 42454359Sroberto/* 42554359Sroberto * propdelay - compute the propagation delay 42654359Sroberto */ 42754359Srobertostatic double 42854359Srobertopropdelay( 42954359Sroberto double dg, 43054359Sroberto double h, 43154359Sroberto int n 43254359Sroberto ) 43354359Sroberto{ 43454359Sroberto double phi; 43554359Sroberto double theta; 43654359Sroberto double td; 43754359Sroberto 43854359Sroberto theta = dg / (EARTHRADIUS * (double)(2 * n)); 43954359Sroberto phi = (PI/2.0) - atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2)); 44054359Sroberto td = dg / (LIGHTSPEED * sin(phi)); 44154359Sroberto if (debug >= 2) 44254359Sroberto printf("propdelay dist %g height %g hops %d time %g\n", 44354359Sroberto dg, h, n, td); 44454359Sroberto return td; 44554359Sroberto} 44654359Sroberto 44754359Sroberto 44854359Sroberto/* 44954359Sroberto * finddelay - find the propagation delay 45054359Sroberto */ 45154359Srobertostatic int 45254359Srobertofinddelay( 45354359Sroberto double lat1, 45454359Sroberto double long1, 45554359Sroberto double lat2, 45654359Sroberto double long2, 45754359Sroberto double h, 45854359Sroberto double *delay 45954359Sroberto ) 46054359Sroberto{ 46154359Sroberto double dg; /* great circle distance */ 46254359Sroberto double delta; /* wave angle */ 46354359Sroberto int n; /* number of hops */ 46454359Sroberto 46554359Sroberto dg = greatcircle(lat1, long1, lat2, long2); 46654359Sroberto if (debug) 46754359Sroberto printf("great circle distance %g km %g miles\n", dg, dg/MILE); 46854359Sroberto 46954359Sroberto n = 1; 47054359Sroberto while ((delta = waveangle(dg, h, n)) < 0.0) { 47154359Sroberto if (debug) 47254359Sroberto printf("tried %d hop%s, no good\n", n, n>1?"s":""); 47354359Sroberto n++; 47454359Sroberto } 47554359Sroberto if (debug) 47654359Sroberto printf("%d hop%s okay, wave angle is %g\n", n, n>1?"s":"", 47754359Sroberto delta / RADPERDEG); 47854359Sroberto 47954359Sroberto *delay = propdelay(dg, h, n); 48054359Sroberto return n; 48154359Sroberto} 48254359Sroberto 48354359Sroberto/* 48454359Sroberto * satdoit - compute a delay and print it 48554359Sroberto */ 48654359Srobertostatic void 48754359Srobertosatdoit( 48854359Sroberto double lat1, 48954359Sroberto double long1, 49054359Sroberto double lat2, 49154359Sroberto double long2, 49254359Sroberto double lat3, 49354359Sroberto double long3, 49454359Sroberto char *str 49554359Sroberto ) 49654359Sroberto{ 49754359Sroberto double up_delay,down_delay; 49854359Sroberto 49954359Sroberto satfinddelay(lat1, long1, lat2, long2, &up_delay); 50054359Sroberto satfinddelay(lat3, long3, lat2, long2, &down_delay); 50154359Sroberto 50254359Sroberto printf("%s, delay %g seconds\n", str, up_delay + down_delay); 50354359Sroberto} 50454359Sroberto 50554359Sroberto/* 50654359Sroberto * satfinddelay - calculate the one-way delay time between a ground station 50754359Sroberto * and a satellite 50854359Sroberto */ 50954359Srobertostatic void 51054359Srobertosatfinddelay( 51154359Sroberto double lat1, 51254359Sroberto double long1, 51354359Sroberto double lat2, 51454359Sroberto double long2, 51554359Sroberto double *delay 51654359Sroberto ) 51754359Sroberto{ 51854359Sroberto double dg; /* great circle distance */ 51954359Sroberto 52054359Sroberto dg = greatcircle(lat1, long1, lat2, long2); 52154359Sroberto 52254359Sroberto *delay = satpropdelay(dg); 52354359Sroberto} 52454359Sroberto 52554359Sroberto/* 52654359Sroberto * satpropdelay - calculate the one-way delay time between a ground station 52754359Sroberto * and a satellite 52854359Sroberto */ 52954359Srobertostatic double 53054359Srobertosatpropdelay( 53154359Sroberto double dg 53254359Sroberto ) 53354359Sroberto{ 53454359Sroberto double k1, k2, dist; 53554359Sroberto double theta; 53654359Sroberto double td; 53754359Sroberto 53854359Sroberto theta = dg / (EARTHRADIUS); 53954359Sroberto k1 = EARTHRADIUS * sin(theta); 54054359Sroberto k2 = SATHEIGHT - (EARTHRADIUS * cos(theta)); 54154359Sroberto if (debug >= 2) 54254359Sroberto printf("Theta %g k1 %g k2 %g\n", theta, k1, k2); 54354359Sroberto dist = sqrt(k1*k1 + k2*k2); 54454359Sroberto td = dist / LIGHTSPEED; 54554359Sroberto if (debug >= 2) 54654359Sroberto printf("propdelay dist %g height %g time %g\n", dg, dist, td); 54754359Sroberto return td; 54854359Sroberto} 549