/contrib/ntp/clockstuff/propdelay.c

https://bitbucket.org/freebsd/freebsd-head/ · C · 544 lines · 383 code · 62 blank · 99 comment · 77 complexity · 76816f3a682d39d5828364328a85d8bf MD5 · raw file

  1. /* propdelay.c,v 3.1 1993/07/06 01:05:24 jbj Exp
  2. * propdelay - compute propagation delays
  3. *
  4. * cc -o propdelay propdelay.c -lm
  5. *
  6. * "Time and Frequency Users' Manual", NBS Technical Note 695 (1977).
  7. */
  8. /*
  9. * This can be used to get a rough idea of the HF propagation delay
  10. * between two points (usually between you and the radio station).
  11. * The usage is
  12. *
  13. * propdelay latitudeA longitudeA latitudeB longitudeB
  14. *
  15. * where points A and B are the locations in question. You obviously
  16. * need to know the latitude and longitude of each of the places.
  17. * The program expects the latitude to be preceded by an 'n' or 's'
  18. * and the longitude to be preceded by an 'e' or 'w'. It understands
  19. * either decimal degrees or degrees:minutes:seconds. Thus to compute
  20. * the delay between the WWVH (21:59:26N, 159:46:00W) and WWV (40:40:49N,
  21. * 105:02:27W) you could use:
  22. *
  23. * propdelay n21:59:26 w159:46 n40:40:49 w105:02:27
  24. *
  25. * By default it prints out a summer (F2 average virtual height 350 km) and
  26. * winter (F2 average virtual height 250 km) number. The results will be
  27. * quite approximate but are about as good as you can do with HF time anyway.
  28. * You might pick a number between the values to use, or use the summer
  29. * value in the summer and switch to the winter value when the static
  30. * above 10 MHz starts to drop off in the fall. You can also use the
  31. * -h switch if you want to specify your own virtual height.
  32. *
  33. * You can also do a
  34. *
  35. * propdelay -W n45:17:47 w75:45:22
  36. *
  37. * to find the propagation delays to WWV and WWVH (from CHU in this
  38. * case), a
  39. *
  40. * propdelay -C n40:40:49 w105:02:27
  41. *
  42. * to find the delays to CHU, and a
  43. *
  44. * propdelay -G n52:03:17 w98:34:18
  45. *
  46. * to find delays to GOES via each of the three satellites.
  47. */
  48. #include <stdio.h>
  49. #include <string.h>
  50. #include "ntp_stdlib.h"
  51. extern double sin (double);
  52. extern double cos (double);
  53. extern double acos (double);
  54. extern double tan (double);
  55. extern double atan (double);
  56. extern double sqrt (double);
  57. #define STREQ(a, b) (*(a) == *(b) && strcmp((a), (b)) == 0)
  58. /*
  59. * Program constants
  60. */
  61. #define EARTHRADIUS (6370.0) /* raduis of earth (km) */
  62. #define LIGHTSPEED (299800.0) /* speed of light, km/s */
  63. #define PI (3.1415926536)
  64. #define RADPERDEG (PI/180.0) /* radians per degree */
  65. #define MILE (1.609344) /* km in a mile */
  66. #define SUMMERHEIGHT (350.0) /* summer height in km */
  67. #define WINTERHEIGHT (250.0) /* winter height in km */
  68. #define SATHEIGHT (6.6110 * 6378.0) /* geosync satellite height in km
  69. from centre of earth */
  70. #define WWVLAT "n40:40:49"
  71. #define WWVLONG "w105:02:27"
  72. #define WWVHLAT "n21:59:26"
  73. #define WWVHLONG "w159:46:00"
  74. #define CHULAT "n45:17:47"
  75. #define CHULONG "w75:45:22"
  76. #define GOES_UP_LAT "n37:52:00"
  77. #define GOES_UP_LONG "w75:27:00"
  78. #define GOES_EAST_LONG "w75:00:00"
  79. #define GOES_STBY_LONG "w105:00:00"
  80. #define GOES_WEST_LONG "w135:00:00"
  81. #define GOES_SAT_LAT "n00:00:00"
  82. char *wwvlat = WWVLAT;
  83. char *wwvlong = WWVLONG;
  84. char *wwvhlat = WWVHLAT;
  85. char *wwvhlong = WWVHLONG;
  86. char *chulat = CHULAT;
  87. char *chulong = CHULONG;
  88. char *goes_up_lat = GOES_UP_LAT;
  89. char *goes_up_long = GOES_UP_LONG;
  90. char *goes_east_long = GOES_EAST_LONG;
  91. char *goes_stby_long = GOES_STBY_LONG;
  92. char *goes_west_long = GOES_WEST_LONG;
  93. char *goes_sat_lat = GOES_SAT_LAT;
  94. int hflag = 0;
  95. int Wflag = 0;
  96. int Cflag = 0;
  97. int Gflag = 0;
  98. int height;
  99. char *progname;
  100. volatile int debug;
  101. static void doit (double, double, double, double, double, char *);
  102. static double latlong (char *, int);
  103. static double greatcircle (double, double, double, double);
  104. static double waveangle (double, double, int);
  105. static double propdelay (double, double, int);
  106. static int finddelay (double, double, double, double, double, double *);
  107. static void satdoit (double, double, double, double, double, double, char *);
  108. static void satfinddelay (double, double, double, double, double *);
  109. static double satpropdelay (double);
  110. /*
  111. * main - parse arguments and handle options
  112. */
  113. int
  114. main(
  115. int argc,
  116. char *argv[]
  117. )
  118. {
  119. int c;
  120. int errflg = 0;
  121. double lat1, long1;
  122. double lat2, long2;
  123. double lat3, long3;
  124. progname = argv[0];
  125. while ((c = ntp_getopt(argc, argv, "dh:CWG")) != EOF)
  126. switch (c) {
  127. case 'd':
  128. ++debug;
  129. break;
  130. case 'h':
  131. hflag++;
  132. height = atof(ntp_optarg);
  133. if (height <= 0.0) {
  134. (void) fprintf(stderr, "height %s unlikely\n",
  135. ntp_optarg);
  136. errflg++;
  137. }
  138. break;
  139. case 'C':
  140. Cflag++;
  141. break;
  142. case 'W':
  143. Wflag++;
  144. break;
  145. case 'G':
  146. Gflag++;
  147. break;
  148. default:
  149. errflg++;
  150. break;
  151. }
  152. if (errflg || (!(Cflag || Wflag || Gflag) && ntp_optind+4 != argc) ||
  153. ((Cflag || Wflag || Gflag) && ntp_optind+2 != argc)) {
  154. (void) fprintf(stderr,
  155. "usage: %s [-d] [-h height] lat1 long1 lat2 long2\n",
  156. progname);
  157. (void) fprintf(stderr," - or -\n");
  158. (void) fprintf(stderr,
  159. "usage: %s -CWG [-d] lat long\n",
  160. progname);
  161. exit(2);
  162. }
  163. if (!(Cflag || Wflag || Gflag)) {
  164. lat1 = latlong(argv[ntp_optind], 1);
  165. long1 = latlong(argv[ntp_optind + 1], 0);
  166. lat2 = latlong(argv[ntp_optind + 2], 1);
  167. long2 = latlong(argv[ntp_optind + 3], 0);
  168. if (hflag) {
  169. doit(lat1, long1, lat2, long2, height, "");
  170. } else {
  171. doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
  172. "summer propagation, ");
  173. doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
  174. "winter propagation, ");
  175. }
  176. } else if (Wflag) {
  177. /*
  178. * Compute delay from WWV
  179. */
  180. lat1 = latlong(argv[ntp_optind], 1);
  181. long1 = latlong(argv[ntp_optind + 1], 0);
  182. lat2 = latlong(wwvlat, 1);
  183. long2 = latlong(wwvlong, 0);
  184. if (hflag) {
  185. doit(lat1, long1, lat2, long2, height, "WWV ");
  186. } else {
  187. doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
  188. "WWV summer propagation, ");
  189. doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
  190. "WWV winter propagation, ");
  191. }
  192. /*
  193. * Compute delay from WWVH
  194. */
  195. lat2 = latlong(wwvhlat, 1);
  196. long2 = latlong(wwvhlong, 0);
  197. if (hflag) {
  198. doit(lat1, long1, lat2, long2, height, "WWVH ");
  199. } else {
  200. doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
  201. "WWVH summer propagation, ");
  202. doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
  203. "WWVH winter propagation, ");
  204. }
  205. } else if (Cflag) {
  206. lat1 = latlong(argv[ntp_optind], 1);
  207. long1 = latlong(argv[ntp_optind + 1], 0);
  208. lat2 = latlong(chulat, 1);
  209. long2 = latlong(chulong, 0);
  210. if (hflag) {
  211. doit(lat1, long1, lat2, long2, height, "CHU ");
  212. } else {
  213. doit(lat1, long1, lat2, long2, (double)SUMMERHEIGHT,
  214. "CHU summer propagation, ");
  215. doit(lat1, long1, lat2, long2, (double)WINTERHEIGHT,
  216. "CHU winter propagation, ");
  217. }
  218. } else if (Gflag) {
  219. lat1 = latlong(goes_up_lat, 1);
  220. long1 = latlong(goes_up_long, 0);
  221. lat3 = latlong(argv[ntp_optind], 1);
  222. long3 = latlong(argv[ntp_optind + 1], 0);
  223. lat2 = latlong(goes_sat_lat, 1);
  224. long2 = latlong(goes_west_long, 0);
  225. satdoit(lat1, long1, lat2, long2, lat3, long3,
  226. "GOES Delay via WEST");
  227. long2 = latlong(goes_stby_long, 0);
  228. satdoit(lat1, long1, lat2, long2, lat3, long3,
  229. "GOES Delay via STBY");
  230. long2 = latlong(goes_east_long, 0);
  231. satdoit(lat1, long1, lat2, long2, lat3, long3,
  232. "GOES Delay via EAST");
  233. }
  234. exit(0);
  235. }
  236. /*
  237. * doit - compute a delay and print it
  238. */
  239. static void
  240. doit(
  241. double lat1,
  242. double long1,
  243. double lat2,
  244. double long2,
  245. double h,
  246. char *str
  247. )
  248. {
  249. int hops;
  250. double delay;
  251. hops = finddelay(lat1, long1, lat2, long2, h, &delay);
  252. printf("%sheight %g km, hops %d, delay %g seconds\n",
  253. str, h, hops, delay);
  254. }
  255. /*
  256. * latlong - decode a latitude/longitude value
  257. */
  258. static double
  259. latlong(
  260. char *str,
  261. int islat
  262. )
  263. {
  264. register char *cp;
  265. register char *bp;
  266. double arg;
  267. double divby;
  268. int isneg;
  269. char buf[32];
  270. char *colon;
  271. if (islat) {
  272. /*
  273. * Must be north or south
  274. */
  275. if (*str == 'N' || *str == 'n')
  276. isneg = 0;
  277. else if (*str == 'S' || *str == 's')
  278. isneg = 1;
  279. else
  280. isneg = -1;
  281. } else {
  282. /*
  283. * East is positive, west is negative
  284. */
  285. if (*str == 'E' || *str == 'e')
  286. isneg = 0;
  287. else if (*str == 'W' || *str == 'w')
  288. isneg = 1;
  289. else
  290. isneg = -1;
  291. }
  292. if (isneg >= 0)
  293. str++;
  294. colon = strchr(str, ':');
  295. if (colon != NULL) {
  296. /*
  297. * in hhh:mm:ss form
  298. */
  299. cp = str;
  300. bp = buf;
  301. while (cp < colon)
  302. *bp++ = *cp++;
  303. *bp = '\0';
  304. cp++;
  305. arg = atof(buf);
  306. divby = 60.0;
  307. colon = strchr(cp, ':');
  308. if (colon != NULL) {
  309. bp = buf;
  310. while (cp < colon)
  311. *bp++ = *cp++;
  312. *bp = '\0';
  313. cp++;
  314. arg += atof(buf) / divby;
  315. divby = 3600.0;
  316. }
  317. if (*cp != '\0')
  318. arg += atof(cp) / divby;
  319. } else {
  320. arg = atof(str);
  321. }
  322. if (isneg == 1)
  323. arg = -arg;
  324. if (debug > 2)
  325. (void) printf("latitude/longitude %s = %g\n", str, arg);
  326. return arg;
  327. }
  328. /*
  329. * greatcircle - compute the great circle distance in kilometers
  330. */
  331. static double
  332. greatcircle(
  333. double lat1,
  334. double long1,
  335. double lat2,
  336. double long2
  337. )
  338. {
  339. double dg;
  340. double l1r, l2r;
  341. l1r = lat1 * RADPERDEG;
  342. l2r = lat2 * RADPERDEG;
  343. dg = EARTHRADIUS * acos(
  344. (cos(l1r) * cos(l2r) * cos((long2-long1)*RADPERDEG))
  345. + (sin(l1r) * sin(l2r)));
  346. if (debug >= 2)
  347. printf(
  348. "greatcircle lat1 %g long1 %g lat2 %g long2 %g dist %g\n",
  349. lat1, long1, lat2, long2, dg);
  350. return dg;
  351. }
  352. /*
  353. * waveangle - compute the wave angle for the given distance, virtual
  354. * height and number of hops.
  355. */
  356. static double
  357. waveangle(
  358. double dg,
  359. double h,
  360. int n
  361. )
  362. {
  363. double theta;
  364. double delta;
  365. theta = dg / (EARTHRADIUS * (double)(2 * n));
  366. delta = atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2)) - theta;
  367. if (debug >= 2)
  368. printf("waveangle dist %g height %g hops %d angle %g\n",
  369. dg, h, n, delta / RADPERDEG);
  370. return delta;
  371. }
  372. /*
  373. * propdelay - compute the propagation delay
  374. */
  375. static double
  376. propdelay(
  377. double dg,
  378. double h,
  379. int n
  380. )
  381. {
  382. double phi;
  383. double theta;
  384. double td;
  385. theta = dg / (EARTHRADIUS * (double)(2 * n));
  386. phi = (PI/2.0) - atan((h / (EARTHRADIUS * sin(theta))) + tan(theta/2));
  387. td = dg / (LIGHTSPEED * sin(phi));
  388. if (debug >= 2)
  389. printf("propdelay dist %g height %g hops %d time %g\n",
  390. dg, h, n, td);
  391. return td;
  392. }
  393. /*
  394. * finddelay - find the propagation delay
  395. */
  396. static int
  397. finddelay(
  398. double lat1,
  399. double long1,
  400. double lat2,
  401. double long2,
  402. double h,
  403. double *delay
  404. )
  405. {
  406. double dg; /* great circle distance */
  407. double delta; /* wave angle */
  408. int n; /* number of hops */
  409. dg = greatcircle(lat1, long1, lat2, long2);
  410. if (debug)
  411. printf("great circle distance %g km %g miles\n", dg, dg/MILE);
  412. n = 1;
  413. while ((delta = waveangle(dg, h, n)) < 0.0) {
  414. if (debug)
  415. printf("tried %d hop%s, no good\n", n, n>1?"s":"");
  416. n++;
  417. }
  418. if (debug)
  419. printf("%d hop%s okay, wave angle is %g\n", n, n>1?"s":"",
  420. delta / RADPERDEG);
  421. *delay = propdelay(dg, h, n);
  422. return n;
  423. }
  424. /*
  425. * satdoit - compute a delay and print it
  426. */
  427. static void
  428. satdoit(
  429. double lat1,
  430. double long1,
  431. double lat2,
  432. double long2,
  433. double lat3,
  434. double long3,
  435. char *str
  436. )
  437. {
  438. double up_delay,down_delay;
  439. satfinddelay(lat1, long1, lat2, long2, &up_delay);
  440. satfinddelay(lat3, long3, lat2, long2, &down_delay);
  441. printf("%s, delay %g seconds\n", str, up_delay + down_delay);
  442. }
  443. /*
  444. * satfinddelay - calculate the one-way delay time between a ground station
  445. * and a satellite
  446. */
  447. static void
  448. satfinddelay(
  449. double lat1,
  450. double long1,
  451. double lat2,
  452. double long2,
  453. double *delay
  454. )
  455. {
  456. double dg; /* great circle distance */
  457. dg = greatcircle(lat1, long1, lat2, long2);
  458. *delay = satpropdelay(dg);
  459. }
  460. /*
  461. * satpropdelay - calculate the one-way delay time between a ground station
  462. * and a satellite
  463. */
  464. static double
  465. satpropdelay(
  466. double dg
  467. )
  468. {
  469. double k1, k2, dist;
  470. double theta;
  471. double td;
  472. theta = dg / (EARTHRADIUS);
  473. k1 = EARTHRADIUS * sin(theta);
  474. k2 = SATHEIGHT - (EARTHRADIUS * cos(theta));
  475. if (debug >= 2)
  476. printf("Theta %g k1 %g k2 %g\n", theta, k1, k2);
  477. dist = sqrt(k1*k1 + k2*k2);
  478. td = dist / LIGHTSPEED;
  479. if (debug >= 2)
  480. printf("propdelay dist %g height %g time %g\n", dg, dist, td);
  481. return td;
  482. }