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 5054359Sroberto#include <stdio.h> 5154359Sroberto#include <string.h> 5254359Sroberto 5354359Sroberto#include "ntp_stdlib.h" 5454359Sroberto 5554359Srobertoextern double sin (double); 5654359Srobertoextern double cos (double); 5754359Srobertoextern double acos (double); 5854359Srobertoextern double tan (double); 5954359Srobertoextern double atan (double); 6054359Srobertoextern double sqrt (double); 6154359Sroberto 6254359Sroberto#define STREQ(a, b) (*(a) == *(b) && strcmp((a), (b)) == 0) 6354359Sroberto 6454359Sroberto/* 6554359Sroberto * Program constants 6654359Sroberto */ 6754359Sroberto#define EARTHRADIUS (6370.0) /* raduis of earth (km) */ 6854359Sroberto#define LIGHTSPEED (299800.0) /* speed of light, km/s */ 6954359Sroberto#define PI (3.1415926536) 7054359Sroberto#define RADPERDEG (PI/180.0) /* radians per degree */ 7154359Sroberto#define MILE (1.609344) /* km in a mile */ 7254359Sroberto 7354359Sroberto#define SUMMERHEIGHT (350.0) /* summer height in km */ 7454359Sroberto#define WINTERHEIGHT (250.0) /* winter height in km */ 7554359Sroberto 7654359Sroberto#define SATHEIGHT (6.6110 * 6378.0) /* geosync satellite height in km 7754359Sroberto from centre of earth */ 7854359Sroberto 7954359Sroberto#define WWVLAT "n40:40:49" 8054359Sroberto#define WWVLONG "w105:02:27" 8154359Sroberto 8254359Sroberto#define WWVHLAT "n21:59:26" 8354359Sroberto#define WWVHLONG "w159:46:00" 8454359Sroberto 8554359Sroberto#define CHULAT "n45:17:47" 8654359Sroberto#define CHULONG "w75:45:22" 8754359Sroberto 8854359Sroberto#define GOES_UP_LAT "n37:52:00" 8954359Sroberto#define GOES_UP_LONG "w75:27:00" 9054359Sroberto#define GOES_EAST_LONG "w75:00:00" 9154359Sroberto#define GOES_STBY_LONG "w105:00:00" 9254359Sroberto#define GOES_WEST_LONG "w135:00:00" 9354359Sroberto#define GOES_SAT_LAT "n00:00:00" 9454359Sroberto 9554359Srobertochar *wwvlat = WWVLAT; 9654359Srobertochar *wwvlong = WWVLONG; 9754359Sroberto 9854359Srobertochar *wwvhlat = WWVHLAT; 9954359Srobertochar *wwvhlong = WWVHLONG; 10054359Sroberto 10154359Srobertochar *chulat = CHULAT; 10254359Srobertochar *chulong = CHULONG; 10354359Sroberto 10454359Srobertochar *goes_up_lat = GOES_UP_LAT; 10554359Srobertochar *goes_up_long = GOES_UP_LONG; 10654359Srobertochar *goes_east_long = GOES_EAST_LONG; 10754359Srobertochar *goes_stby_long = GOES_STBY_LONG; 10854359Srobertochar *goes_west_long = GOES_WEST_LONG; 10954359Srobertochar *goes_sat_lat = GOES_SAT_LAT; 11054359Sroberto 11154359Srobertoint hflag = 0; 11254359Srobertoint Wflag = 0; 11354359Srobertoint Cflag = 0; 11454359Srobertoint Gflag = 0; 11554359Srobertoint height; 11654359Sroberto 11754359Srobertochar *progname; 118182007Srobertovolatile int debug; 11954359Sroberto 12054359Srobertostatic void doit (double, double, double, double, double, char *); 12154359Srobertostatic double latlong (char *, int); 12254359Srobertostatic double greatcircle (double, double, double, double); 12354359Srobertostatic double waveangle (double, double, int); 12454359Srobertostatic double propdelay (double, double, int); 12554359Srobertostatic int finddelay (double, double, double, double, double, double *); 12654359Srobertostatic void satdoit (double, double, double, double, double, double, char *); 12754359Srobertostatic void satfinddelay (double, double, double, double, double *); 12854359Srobertostatic double satpropdelay (double); 12954359Sroberto 13054359Sroberto/* 13154359Sroberto * main - parse arguments and handle options 13254359Sroberto */ 13354359Srobertoint 13454359Srobertomain( 13554359Sroberto int argc, 13654359Sroberto char *argv[] 13754359Sroberto ) 13854359Sroberto{ 13954359Sroberto int c; 14054359Sroberto int errflg = 0; 14154359Sroberto double lat1, long1; 14254359Sroberto double lat2, long2; 14354359Sroberto double lat3, long3; 14454359Sroberto 14554359Sroberto progname = argv[0]; 14654359Sroberto while ((c = ntp_getopt(argc, argv, "dh:CWG")) != EOF) 14754359Sroberto switch (c) { 14854359Sroberto case 'd': 14954359Sroberto ++debug; 15054359Sroberto break; 15154359Sroberto case 'h': 15254359Sroberto hflag++; 15354359Sroberto height = atof(ntp_optarg); 15454359Sroberto if (height <= 0.0) { 15554359Sroberto (void) fprintf(stderr, "height %s unlikely\n", 15654359Sroberto ntp_optarg); 15754359Sroberto errflg++; 15854359Sroberto } 15954359Sroberto break; 16054359Sroberto case 'C': 16154359Sroberto Cflag++; 16254359Sroberto break; 16354359Sroberto case 'W': 16454359Sroberto Wflag++; 16554359Sroberto break; 16654359Sroberto case 'G': 16754359Sroberto Gflag++; 16854359Sroberto break; 16954359Sroberto default: 17054359Sroberto errflg++; 17154359Sroberto break; 17254359Sroberto } 17354359Sroberto if (errflg || (!(Cflag || Wflag || Gflag) && ntp_optind+4 != argc) || 17454359Sroberto ((Cflag || Wflag || Gflag) && ntp_optind+2 != argc)) { 17554359Sroberto (void) fprintf(stderr, 17654359Sroberto "usage: %s [-d] [-h height] lat1 long1 lat2 long2\n", 17754359Sroberto progname); 17854359Sroberto (void) fprintf(stderr," - or -\n"); 17954359Sroberto (void) fprintf(stderr, 18054359Sroberto "usage: %s -CWG [-d] lat long\n", 18154359Sroberto progname); 18254359Sroberto exit(2); 18354359Sroberto } 18454359Sroberto 18554359Sroberto 18654359Sroberto if (!(Cflag || Wflag || Gflag)) { 18754359Sroberto lat1 = latlong(argv[ntp_optind], 1); 18854359Sroberto long1 = latlong(argv[ntp_optind + 1], 0); 18954359Sroberto lat2 = latlong(argv[ntp_optind + 2], 1); 19054359Sroberto long2 = latlong(argv[ntp_optind + 3], 0); 19154359Sroberto if (hflag) { 19254359Sroberto doit(lat1, long1, lat2, long2, height, ""); 19354359Sroberto } else { 19454359Sroberto doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 19554359Sroberto "summer propagation, "); 19654359Sroberto doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 19754359Sroberto "winter propagation, "); 19854359Sroberto } 19954359Sroberto } else if (Wflag) { 20054359Sroberto /* 20154359Sroberto * Compute delay from WWV 20254359Sroberto */ 20354359Sroberto lat1 = latlong(argv[ntp_optind], 1); 20454359Sroberto long1 = latlong(argv[ntp_optind + 1], 0); 20554359Sroberto lat2 = latlong(wwvlat, 1); 20654359Sroberto long2 = latlong(wwvlong, 0); 20754359Sroberto if (hflag) { 20854359Sroberto doit(lat1, long1, lat2, long2, height, "WWV "); 20954359Sroberto } else { 21054359Sroberto doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 21154359Sroberto "WWV summer propagation, "); 21254359Sroberto doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 21354359Sroberto "WWV winter propagation, "); 21454359Sroberto } 21554359Sroberto 21654359Sroberto /* 21754359Sroberto * Compute delay from WWVH 21854359Sroberto */ 21954359Sroberto lat2 = latlong(wwvhlat, 1); 22054359Sroberto long2 = latlong(wwvhlong, 0); 22154359Sroberto if (hflag) { 22254359Sroberto doit(lat1, long1, lat2, long2, height, "WWVH "); 22354359Sroberto } else { 22454359Sroberto doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 22554359Sroberto "WWVH summer propagation, "); 22654359Sroberto doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 22754359Sroberto "WWVH winter propagation, "); 22854359Sroberto } 22954359Sroberto } else if (Cflag) { 23054359Sroberto lat1 = latlong(argv[ntp_optind], 1); 23154359Sroberto long1 = latlong(argv[ntp_optind + 1], 0); 23254359Sroberto lat2 = latlong(chulat, 1); 23354359Sroberto long2 = latlong(chulong, 0); 23454359Sroberto if (hflag) { 23554359Sroberto doit(lat1, long1, lat2, long2, height, "CHU "); 23654359Sroberto } else { 23754359Sroberto doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT, 23854359Sroberto "CHU summer propagation, "); 23954359Sroberto doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT, 24054359Sroberto "CHU winter propagation, "); 24154359Sroberto } 24254359Sroberto } else if (Gflag) { 24354359Sroberto lat1 = latlong(goes_up_lat, 1); 24454359Sroberto long1 = latlong(goes_up_long, 0); 24554359Sroberto lat3 = latlong(argv[ntp_optind], 1); 24654359Sroberto long3 = latlong(argv[ntp_optind + 1], 0); 24754359Sroberto 24854359Sroberto lat2 = latlong(goes_sat_lat, 1); 24954359Sroberto 25054359Sroberto long2 = latlong(goes_west_long, 0); 25154359Sroberto satdoit(lat1, long1, lat2, long2, lat3, long3, 25254359Sroberto "GOES Delay via WEST"); 25354359Sroberto 25454359Sroberto long2 = latlong(goes_stby_long, 0); 25554359Sroberto satdoit(lat1, long1, lat2, long2, lat3, long3, 25654359Sroberto "GOES Delay via STBY"); 25754359Sroberto 25854359Sroberto long2 = latlong(goes_east_long, 0); 25954359Sroberto satdoit(lat1, long1, lat2, long2, lat3, long3, 26054359Sroberto "GOES Delay via EAST"); 26154359Sroberto 26254359Sroberto } 26354359Sroberto exit(0); 26454359Sroberto} 26554359Sroberto 26654359Sroberto 26754359Sroberto/* 26854359Sroberto * doit - compute a delay and print it 26954359Sroberto */ 27054359Srobertostatic void 27154359Srobertodoit( 27254359Sroberto double lat1, 27354359Sroberto double long1, 27454359Sroberto double lat2, 27554359Sroberto double long2, 27654359Sroberto double h, 27754359Sroberto char *str 27854359Sroberto ) 27954359Sroberto{ 28054359Sroberto int hops; 28154359Sroberto double delay; 28254359Sroberto 28354359Sroberto hops = finddelay(lat1, long1, lat2, long2, h, &delay); 28454359Sroberto printf("%sheight %g km, hops %d, delay %g seconds\n", 28554359Sroberto str, h, hops, delay); 28654359Sroberto} 28754359Sroberto 28854359Sroberto 28954359Sroberto/* 29054359Sroberto * latlong - decode a latitude/longitude value 29154359Sroberto */ 29254359Srobertostatic double 29354359Srobertolatlong( 29454359Sroberto char *str, 29554359Sroberto int islat 29654359Sroberto ) 29754359Sroberto{ 29854359Sroberto register char *cp; 29954359Sroberto register char *bp; 30054359Sroberto double arg; 301182007Sroberto double divby; 30254359Sroberto int isneg; 30354359Sroberto char buf[32]; 30454359Sroberto char *colon; 30554359Sroberto 30654359Sroberto if (islat) { 30754359Sroberto /* 30854359Sroberto * Must be north or south 30954359Sroberto */ 31054359Sroberto if (*str == 'N' || *str == 'n') 31154359Sroberto isneg = 0; 31254359Sroberto else if (*str == 'S' || *str == 's') 31354359Sroberto isneg = 1; 31454359Sroberto else 31554359Sroberto isneg = -1; 31654359Sroberto } else { 31754359Sroberto /* 31854359Sroberto * East is positive, west is negative 31954359Sroberto */ 32054359Sroberto if (*str == 'E' || *str == 'e') 32154359Sroberto isneg = 0; 32254359Sroberto else if (*str == 'W' || *str == 'w') 32354359Sroberto isneg = 1; 32454359Sroberto else 32554359Sroberto isneg = -1; 32654359Sroberto } 32754359Sroberto 32854359Sroberto if (isneg >= 0) 32954359Sroberto str++; 33054359Sroberto 33154359Sroberto colon = strchr(str, ':'); 33254359Sroberto if (colon != NULL) { 33354359Sroberto /* 33454359Sroberto * in hhh:mm:ss form 33554359Sroberto */ 33654359Sroberto cp = str; 33754359Sroberto bp = buf; 33854359Sroberto while (cp < colon) 33954359Sroberto *bp++ = *cp++; 34054359Sroberto *bp = '\0'; 34154359Sroberto cp++; 34254359Sroberto arg = atof(buf); 343182007Sroberto divby = 60.0; 34454359Sroberto colon = strchr(cp, ':'); 34554359Sroberto if (colon != NULL) { 34654359Sroberto bp = buf; 34754359Sroberto while (cp < colon) 34854359Sroberto *bp++ = *cp++; 34954359Sroberto *bp = '\0'; 35054359Sroberto cp++; 351182007Sroberto arg += atof(buf) / divby; 352182007Sroberto divby = 3600.0; 35354359Sroberto } 35454359Sroberto if (*cp != '\0') 355182007Sroberto arg += atof(cp) / divby; 35654359Sroberto } else { 35754359Sroberto arg = atof(str); 35854359Sroberto } 35954359Sroberto 36054359Sroberto if (isneg == 1) 36154359Sroberto arg = -arg; 36254359Sroberto 36354359Sroberto if (debug > 2) 36454359Sroberto (void) printf("latitude/longitude %s = %g\n", str, arg); 36554359Sroberto 36654359Sroberto return arg; 36754359Sroberto} 36854359Sroberto 36954359Sroberto 37054359Sroberto/* 37154359Sroberto * greatcircle - compute the great circle distance in kilometers 37254359Sroberto */ 37354359Srobertostatic double 37454359Srobertogreatcircle( 37554359Sroberto double lat1, 37654359Sroberto double long1, 37754359Sroberto double lat2, 37854359Sroberto double long2 37954359Sroberto ) 38054359Sroberto{ 38154359Sroberto double dg; 38254359Sroberto double l1r, l2r; 38354359Sroberto 38454359Sroberto l1r = lat1 * RADPERDEG; 38554359Sroberto l2r = lat2 * RADPERDEG; 38654359Sroberto dg = EARTHRADIUS * acos( 38754359Sroberto (cos(l1r) * cos(l2r) * cos((long2-long1)*RADPERDEG)) 38854359Sroberto + (sin(l1r) * sin(l2r))); 38954359Sroberto if (debug >= 2) 39054359Sroberto printf( 39154359Sroberto "greatcircle lat1 %g long1 %g lat2 %g long2 %g dist %g\n", 39254359Sroberto lat1, long1, lat2, long2, dg); 39354359Sroberto return dg; 39454359Sroberto} 39554359Sroberto 39654359Sroberto 39754359Sroberto/* 39854359Sroberto * waveangle - compute the wave angle for the given distance, virtual 39954359Sroberto * height and number of hops. 40054359Sroberto */ 40154359Srobertostatic double 40254359Srobertowaveangle( 40354359Sroberto double dg, 40454359Sroberto double h, 40554359Sroberto int n 40654359Sroberto ) 40754359Sroberto{ 40854359Sroberto double theta; 40954359Sroberto double delta; 41054359Sroberto 41154359Sroberto theta = dg / (EARTHRADIUS * (double)(2 * n)); 41254359Sroberto delta = atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2)) - theta; 41354359Sroberto if (debug >= 2) 41454359Sroberto printf("waveangle dist %g height %g hops %d angle %g\n", 41554359Sroberto dg, h, n, delta / RADPERDEG); 41654359Sroberto return delta; 41754359Sroberto} 41854359Sroberto 41954359Sroberto 42054359Sroberto/* 42154359Sroberto * propdelay - compute the propagation delay 42254359Sroberto */ 42354359Srobertostatic double 42454359Srobertopropdelay( 42554359Sroberto double dg, 42654359Sroberto double h, 42754359Sroberto int n 42854359Sroberto ) 42954359Sroberto{ 43054359Sroberto double phi; 43154359Sroberto double theta; 43254359Sroberto double td; 43354359Sroberto 43454359Sroberto theta = dg / (EARTHRADIUS * (double)(2 * n)); 43554359Sroberto phi = (PI/2.0) - atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2)); 43654359Sroberto td = dg / (LIGHTSPEED * sin(phi)); 43754359Sroberto if (debug >= 2) 43854359Sroberto printf("propdelay dist %g height %g hops %d time %g\n", 43954359Sroberto dg, h, n, td); 44054359Sroberto return td; 44154359Sroberto} 44254359Sroberto 44354359Sroberto 44454359Sroberto/* 44554359Sroberto * finddelay - find the propagation delay 44654359Sroberto */ 44754359Srobertostatic int 44854359Srobertofinddelay( 44954359Sroberto double lat1, 45054359Sroberto double long1, 45154359Sroberto double lat2, 45254359Sroberto double long2, 45354359Sroberto double h, 45454359Sroberto double *delay 45554359Sroberto ) 45654359Sroberto{ 45754359Sroberto double dg; /* great circle distance */ 45854359Sroberto double delta; /* wave angle */ 45954359Sroberto int n; /* number of hops */ 46054359Sroberto 46154359Sroberto dg = greatcircle(lat1, long1, lat2, long2); 46254359Sroberto if (debug) 46354359Sroberto printf("great circle distance %g km %g miles\n", dg, dg/MILE); 46454359Sroberto 46554359Sroberto n = 1; 46654359Sroberto while ((delta = waveangle(dg, h, n)) < 0.0) { 46754359Sroberto if (debug) 46854359Sroberto printf("tried %d hop%s, no good\n", n, n>1?"s":""); 46954359Sroberto n++; 47054359Sroberto } 47154359Sroberto if (debug) 47254359Sroberto printf("%d hop%s okay, wave angle is %g\n", n, n>1?"s":"", 47354359Sroberto delta / RADPERDEG); 47454359Sroberto 47554359Sroberto *delay = propdelay(dg, h, n); 47654359Sroberto return n; 47754359Sroberto} 47854359Sroberto 47954359Sroberto/* 48054359Sroberto * satdoit - compute a delay and print it 48154359Sroberto */ 48254359Srobertostatic void 48354359Srobertosatdoit( 48454359Sroberto double lat1, 48554359Sroberto double long1, 48654359Sroberto double lat2, 48754359Sroberto double long2, 48854359Sroberto double lat3, 48954359Sroberto double long3, 49054359Sroberto char *str 49154359Sroberto ) 49254359Sroberto{ 49354359Sroberto double up_delay,down_delay; 49454359Sroberto 49554359Sroberto satfinddelay(lat1, long1, lat2, long2, &up_delay); 49654359Sroberto satfinddelay(lat3, long3, lat2, long2, &down_delay); 49754359Sroberto 49854359Sroberto printf("%s, delay %g seconds\n", str, up_delay + down_delay); 49954359Sroberto} 50054359Sroberto 50154359Sroberto/* 50254359Sroberto * satfinddelay - calculate the one-way delay time between a ground station 50354359Sroberto * and a satellite 50454359Sroberto */ 50554359Srobertostatic void 50654359Srobertosatfinddelay( 50754359Sroberto double lat1, 50854359Sroberto double long1, 50954359Sroberto double lat2, 51054359Sroberto double long2, 51154359Sroberto double *delay 51254359Sroberto ) 51354359Sroberto{ 51454359Sroberto double dg; /* great circle distance */ 51554359Sroberto 51654359Sroberto dg = greatcircle(lat1, long1, lat2, long2); 51754359Sroberto 51854359Sroberto *delay = satpropdelay(dg); 51954359Sroberto} 52054359Sroberto 52154359Sroberto/* 52254359Sroberto * satpropdelay - calculate the one-way delay time between a ground station 52354359Sroberto * and a satellite 52454359Sroberto */ 52554359Srobertostatic double 52654359Srobertosatpropdelay( 52754359Sroberto double dg 52854359Sroberto ) 52954359Sroberto{ 53054359Sroberto double k1, k2, dist; 53154359Sroberto double theta; 53254359Sroberto double td; 53354359Sroberto 53454359Sroberto theta = dg / (EARTHRADIUS); 53554359Sroberto k1 = EARTHRADIUS * sin(theta); 53654359Sroberto k2 = SATHEIGHT - (EARTHRADIUS * cos(theta)); 53754359Sroberto if (debug >= 2) 53854359Sroberto printf("Theta %g k1 %g k2 %g\n", theta, k1, k2); 53954359Sroberto dist = sqrt(k1*k1 + k2*k2); 54054359Sroberto td = dist / LIGHTSPEED; 54154359Sroberto if (debug >= 2) 54254359Sroberto printf("propdelay dist %g height %g time %g\n", dg, dist, td); 54354359Sroberto return td; 54454359Sroberto} 545