#include <stdio.h>
#include <string.h>
#include "ntp_stdlib.h"
extern double sin (double);
extern double cos (double);
extern double acos (double);
extern double tan (double);
extern double atan (double);
extern double sqrt (double);
#define STREQ(a, b) (*(a) == *(b) && strcmp((a), (b)) == 0)
#define EARTHRADIUS (6370.0)
#define LIGHTSPEED (299800.0)
#define PI (3.1415926536)
#define RADPERDEG (PI/180.0)
#define MILE (1.609344)
#define SUMMERHEIGHT (350.0)
#define WINTERHEIGHT (250.0)
#define SATHEIGHT (6.6110 * 6378.0)
#define WWVLAT "n40:40:49"
#define WWVLONG "w105:02:27"
#define WWVHLAT "n21:59:26"
#define WWVHLONG "w159:46:00"
#define CHULAT "n45:17:47"
#define CHULONG "w75:45:22"
#define GOES_UP_LAT "n37:52:00"
#define GOES_UP_LONG "w75:27:00"
#define GOES_EAST_LONG "w75:00:00"
#define GOES_STBY_LONG "w105:00:00"
#define GOES_WEST_LONG "w135:00:00"
#define GOES_SAT_LAT "n00:00:00"
char *wwvlat = WWVLAT;
char *wwvlong = WWVLONG;
char *wwvhlat = WWVHLAT;
char *wwvhlong = WWVHLONG;
char *chulat = CHULAT;
char *chulong = CHULONG;
char *goes_up_lat = GOES_UP_LAT;
char *goes_up_long = GOES_UP_LONG;
char *goes_east_long = GOES_EAST_LONG;
char *goes_stby_long = GOES_STBY_LONG;
char *goes_west_long = GOES_WEST_LONG;
char *goes_sat_lat = GOES_SAT_LAT;
int hflag = 0;
int Wflag = 0;
int Cflag = 0;
int Gflag = 0;
int height;
char *progname;
volatile int debug;
static void doit (double, double, double, double, double, char *);
static double latlong (char *, int);
static double greatcircle (double, double, double, double);
static double waveangle (double, double, int);
static double propdelay (double, double, int);
static int finddelay (double, double, double, double, double, double *);
static void satdoit (double, double, double, double, double, double, char *);
static void satfinddelay (double, double, double, double, double *);
static double satpropdelay (double);
int
main(
int argc,
char *argv[]
)
{
int c;
int errflg = 0;
double lat1, long1;
double lat2, long2;
double lat3, long3;
progname = argv[0];
while ((c = ntp_getopt(argc, argv, "dh:CWG")) != EOF)
switch (c) {
case 'd':
++debug;
break;
case 'h':
hflag++;
height = atof(ntp_optarg);
if (height <= 0.0) {
(void) fprintf(stderr, "height %s unlikely\n",
ntp_optarg);
errflg++;
}
break;
case 'C':
Cflag++;
break;
case 'W':
Wflag++;
break;
case 'G':
Gflag++;
break;
default:
errflg++;
break;
}
if (errflg || (!(Cflag || Wflag || Gflag) && ntp_optind+4 != argc) ||
((Cflag || Wflag || Gflag) && ntp_optind+2 != argc)) {
(void) fprintf(stderr,
"usage: %s [-d] [-h height] lat1 long1 lat2 long2\n",
progname);
(void) fprintf(stderr," - or -\n");
(void) fprintf(stderr,
"usage: %s -CWG [-d] lat long\n",
progname);
exit(2);
}
if (!(Cflag || Wflag || Gflag)) {
lat1 = latlong(argv[ntp_optind], 1);
long1 = latlong(argv[ntp_optind + 1], 0);
lat2 = latlong(argv[ntp_optind + 2], 1);
long2 = latlong(argv[ntp_optind + 3], 0);
if (hflag) {
doit(lat1, long1, lat2, long2, height, "");
} else {
doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
"summer propagation, ");
doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
"winter propagation, ");
}
} else if (Wflag) {
lat1 = latlong(argv[ntp_optind], 1);
long1 = latlong(argv[ntp_optind + 1], 0);
lat2 = latlong(wwvlat, 1);
long2 = latlong(wwvlong, 0);
if (hflag) {
doit(lat1, long1, lat2, long2, height, "WWV ");
} else {
doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
"WWV summer propagation, ");
doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
"WWV winter propagation, ");
}
lat2 = latlong(wwvhlat, 1);
long2 = latlong(wwvhlong, 0);
if (hflag) {
doit(lat1, long1, lat2, long2, height, "WWVH ");
} else {
doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
"WWVH summer propagation, ");
doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
"WWVH winter propagation, ");
}
} else if (Cflag) {
lat1 = latlong(argv[ntp_optind], 1);
long1 = latlong(argv[ntp_optind + 1], 0);
lat2 = latlong(chulat, 1);
long2 = latlong(chulong, 0);
if (hflag) {
doit(lat1, long1, lat2, long2, height, "CHU ");
} else {
doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
"CHU summer propagation, ");
doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
"CHU winter propagation, ");
}
} else if (Gflag) {
lat1 = latlong(goes_up_lat, 1);
long1 = latlong(goes_up_long, 0);
lat3 = latlong(argv[ntp_optind], 1);
long3 = latlong(argv[ntp_optind + 1], 0);
lat2 = latlong(goes_sat_lat, 1);
long2 = latlong(goes_west_long, 0);
satdoit(lat1, long1, lat2, long2, lat3, long3,
"GOES Delay via WEST");
long2 = latlong(goes_stby_long, 0);
satdoit(lat1, long1, lat2, long2, lat3, long3,
"GOES Delay via STBY");
long2 = latlong(goes_east_long, 0);
satdoit(lat1, long1, lat2, long2, lat3, long3,
"GOES Delay via EAST");
}
exit(0);
}
static void
doit(
double lat1,
double long1,
double lat2,
double long2,
double h,
char *str
)
{
int hops;
double delay;
hops = finddelay(lat1, long1, lat2, long2, h, &delay);
printf("%sheight %g km, hops %d, delay %g seconds\n",
str, h, hops, delay);
}
static double
latlong(
char *str,
int islat
)
{
register char *cp;
register char *bp;
double arg;
double divby;
int isneg;
char buf[32];
char *colon;
if (islat) {
if (*str == 'N' || *str == 'n')
isneg = 0;
else if (*str == 'S' || *str == 's')
isneg = 1;
else
isneg = -1;
} else {
if (*str == 'E' || *str == 'e')
isneg = 0;
else if (*str == 'W' || *str == 'w')
isneg = 1;
else
isneg = -1;
}
if (isneg >= 0)
str++;
colon = strchr(str, ':');
if (colon != NULL) {
cp = str;
bp = buf;
while (cp < colon)
*bp++ = *cp++;
*bp = '\0';
cp++;
arg = atof(buf);
divby = 60.0;
colon = strchr(cp, ':');
if (colon != NULL) {
bp = buf;
while (cp < colon)
*bp++ = *cp++;
*bp = '\0';
cp++;
arg += atof(buf) / divby;
divby = 3600.0;
}
if (*cp != '\0')
arg += atof(cp) / divby;
} else {
arg = atof(str);
}
if (isneg == 1)
arg = -arg;
if (debug > 2)
(void) printf("latitude/longitude %s = %g\n", str, arg);
return arg;
}
static double
greatcircle(
double lat1,
double long1,
double lat2,
double long2
)
{
double dg;
double l1r, l2r;
l1r = lat1 * RADPERDEG;
l2r = lat2 * RADPERDEG;
dg = EARTHRADIUS * acos(
(cos(l1r) * cos(l2r) * cos((long2-long1)*RADPERDEG))
+ (sin(l1r) * sin(l2r)));
if (debug >= 2)
printf(
"greatcircle lat1 %g long1 %g lat2 %g long2 %g dist %g\n",
lat1, long1, lat2, long2, dg);
return dg;
}
static double
waveangle(
double dg,
double h,
int n
)
{
double theta;
double delta;
theta = dg / (EARTHRADIUS * (double)(2 * n));
delta = atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2)) - theta;
if (debug >= 2)
printf("waveangle dist %g height %g hops %d angle %g\n",
dg, h, n, delta / RADPERDEG);
return delta;
}
static double
propdelay(
double dg,
double h,
int n
)
{
double phi;
double theta;
double td;
theta = dg / (EARTHRADIUS * (double)(2 * n));
phi = (PI/2.0) - atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2));
td = dg / (LIGHTSPEED * sin(phi));
if (debug >= 2)
printf("propdelay dist %g height %g hops %d time %g\n",
dg, h, n, td);
return td;
}
static int
finddelay(
double lat1,
double long1,
double lat2,
double long2,
double h,
double *delay
)
{
double dg;
double delta;
int n;
dg = greatcircle(lat1, long1, lat2, long2);
if (debug)
printf("great circle distance %g km %g miles\n", dg, dg/MILE);
n = 1;
while ((delta = waveangle(dg, h, n)) < 0.0) {
if (debug)
printf("tried %d hop%s, no good\n", n, n>1?"s":"");
n++;
}
if (debug)
printf("%d hop%s okay, wave angle is %g\n", n, n>1?"s":"",
delta / RADPERDEG);
*delay = propdelay(dg, h, n);
return n;
}
static void
satdoit(
double lat1,
double long1,
double lat2,
double long2,
double lat3,
double long3,
char *str
)
{
double up_delay,down_delay;
satfinddelay(lat1, long1, lat2, long2, &up_delay);
satfinddelay(lat3, long3, lat2, long2, &down_delay);
printf("%s, delay %g seconds\n", str, up_delay + down_delay);
}
static void
satfinddelay(
double lat1,
double long1,
double lat2,
double long2,
double *delay
)
{
double dg;
dg = greatcircle(lat1, long1, lat2, long2);
*delay = satpropdelay(dg);
}
static double
satpropdelay(
double dg
)
{
double k1, k2, dist;
double theta;
double td;
theta = dg / (EARTHRADIUS);
k1 = EARTHRADIUS * sin(theta);
k2 = SATHEIGHT - (EARTHRADIUS * cos(theta));
if (debug >= 2)
printf("Theta %g k1 %g k2 %g\n", theta, k1, k2);
dist = sqrt(k1*k1 + k2*k2);
td = dist / LIGHTSPEED;
if (debug >= 2)
printf("propdelay dist %g height %g time %g\n", dg, dist, td);
return td;
}