547 lines
12 KiB
C
547 lines
12 KiB
C
/* $NetBSD: propdelay.c,v 1.3 2006/06/11 19:34:09 kardel Exp $ */
|
|
|
|
/* propdelay.c,v 3.1 1993/07/06 01:05:24 jbj Exp
|
|
* propdelay - compute propagation delays
|
|
*
|
|
* cc -o propdelay propdelay.c -lm
|
|
*
|
|
* "Time and Frequency Users' Manual", NBS Technical Note 695 (1977).
|
|
*/
|
|
|
|
/*
|
|
* This can be used to get a rough idea of the HF propagation delay
|
|
* between two points (usually between you and the radio station).
|
|
* The usage is
|
|
*
|
|
* propdelay latitudeA longitudeA latitudeB longitudeB
|
|
*
|
|
* where points A and B are the locations in question. You obviously
|
|
* need to know the latitude and longitude of each of the places.
|
|
* The program expects the latitude to be preceded by an 'n' or 's'
|
|
* and the longitude to be preceded by an 'e' or 'w'. It understands
|
|
* either decimal degrees or degrees:minutes:seconds. Thus to compute
|
|
* the delay between the WWVH (21:59:26N, 159:46:00W) and WWV (40:40:49N,
|
|
* 105:02:27W) you could use:
|
|
*
|
|
* propdelay n21:59:26 w159:46 n40:40:49 w105:02:27
|
|
*
|
|
* By default it prints out a summer (F2 average virtual height 350 km) and
|
|
* winter (F2 average virtual height 250 km) number. The results will be
|
|
* quite approximate but are about as good as you can do with HF time anyway.
|
|
* You might pick a number between the values to use, or use the summer
|
|
* value in the summer and switch to the winter value when the static
|
|
* above 10 MHz starts to drop off in the fall. You can also use the
|
|
* -h switch if you want to specify your own virtual height.
|
|
*
|
|
* You can also do a
|
|
*
|
|
* propdelay -W n45:17:47 w75:45:22
|
|
*
|
|
* to find the propagation delays to WWV and WWVH (from CHU in this
|
|
* case), a
|
|
*
|
|
* propdelay -C n40:40:49 w105:02:27
|
|
*
|
|
* to find the delays to CHU, and a
|
|
*
|
|
* propdelay -G n52:03:17 w98:34:18
|
|
*
|
|
* to find delays to GOES via each of the three satellites.
|
|
*/
|
|
|
|
#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)
|
|
|
|
/*
|
|
* Program constants
|
|
*/
|
|
#define EARTHRADIUS (6370.0) /* raduis of earth (km) */
|
|
#define LIGHTSPEED (299800.0) /* speed of light, km/s */
|
|
#define PI (3.1415926536)
|
|
#define RADPERDEG (PI/180.0) /* radians per degree */
|
|
#define MILE (1.609344) /* km in a mile */
|
|
|
|
#define SUMMERHEIGHT (350.0) /* summer height in km */
|
|
#define WINTERHEIGHT (250.0) /* winter height in km */
|
|
|
|
#define SATHEIGHT (6.6110 * 6378.0) /* geosync satellite height in km
|
|
from centre of earth */
|
|
|
|
#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);
|
|
|
|
/*
|
|
* main - parse arguments and handle options
|
|
*/
|
|
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) {
|
|
/*
|
|
* Compute delay from WWV
|
|
*/
|
|
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, ");
|
|
}
|
|
|
|
/*
|
|
* Compute delay from WWVH
|
|
*/
|
|
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);
|
|
}
|
|
|
|
|
|
/*
|
|
* doit - compute a delay and print it
|
|
*/
|
|
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);
|
|
}
|
|
|
|
|
|
/*
|
|
* latlong - decode a latitude/longitude value
|
|
*/
|
|
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) {
|
|
/*
|
|
* Must be north or south
|
|
*/
|
|
if (*str == 'N' || *str == 'n')
|
|
isneg = 0;
|
|
else if (*str == 'S' || *str == 's')
|
|
isneg = 1;
|
|
else
|
|
isneg = -1;
|
|
} else {
|
|
/*
|
|
* East is positive, west is negative
|
|
*/
|
|
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) {
|
|
/*
|
|
* in hhh:mm:ss form
|
|
*/
|
|
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;
|
|
}
|
|
|
|
|
|
/*
|
|
* greatcircle - compute the great circle distance in kilometers
|
|
*/
|
|
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;
|
|
}
|
|
|
|
|
|
/*
|
|
* waveangle - compute the wave angle for the given distance, virtual
|
|
* height and number of hops.
|
|
*/
|
|
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;
|
|
}
|
|
|
|
|
|
/*
|
|
* propdelay - compute the propagation delay
|
|
*/
|
|
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;
|
|
}
|
|
|
|
|
|
/*
|
|
* finddelay - find the propagation delay
|
|
*/
|
|
static int
|
|
finddelay(
|
|
double lat1,
|
|
double long1,
|
|
double lat2,
|
|
double long2,
|
|
double h,
|
|
double *delay
|
|
)
|
|
{
|
|
double dg; /* great circle distance */
|
|
double delta; /* wave angle */
|
|
int n; /* number of hops */
|
|
|
|
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;
|
|
}
|
|
|
|
/*
|
|
* satdoit - compute a delay and print it
|
|
*/
|
|
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);
|
|
}
|
|
|
|
/*
|
|
* satfinddelay - calculate the one-way delay time between a ground station
|
|
* and a satellite
|
|
*/
|
|
static void
|
|
satfinddelay(
|
|
double lat1,
|
|
double long1,
|
|
double lat2,
|
|
double long2,
|
|
double *delay
|
|
)
|
|
{
|
|
double dg; /* great circle distance */
|
|
|
|
dg = greatcircle(lat1, long1, lat2, long2);
|
|
|
|
*delay = satpropdelay(dg);
|
|
}
|
|
|
|
/*
|
|
* satpropdelay - calculate the one-way delay time between a ground station
|
|
* and a satellite
|
|
*/
|
|
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;
|
|
}
|