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