home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume19 / xephem / part21 / riset.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-15  |  4.8 KB  |  139 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "astro.h"
  4.  
  5. #if defined(__STDC__) || defined(__cplusplus)
  6. #define P_(s) s
  7. #else
  8. #define P_(s) ()
  9. #endif
  10.  
  11. extern void range P_((double *v, double r));
  12.  
  13. void riset P_((double ra, double dec, double lat, double dis, double *lstr, double *lsts, double *azr, double *azs, int *status));
  14.  
  15. #undef P_
  16.  
  17.  
  18. /* given the true geocentric ra and dec of an object, the observer's latitude,
  19.  *   lat, and a horizon displacement correction, dis, all in radians, find the
  20.  *   local sidereal times and azimuths of rising and setting, lstr/s
  21.  *   and azr/s, also all in radians, respectively.
  22.  * dis is the vertical displacement from the true position of the horizon. it
  23.  *   is positive if the apparent position is lower than the true position.
  24.  *   said another way, it is positive if the shift causes the object to spend
  25.  *   longer above the horizon. for example, atmospheric refraction is typically
  26.  *   assumed to produce a vertical shift of 34 arc minutes at the horizon; dis
  27.  *   would then take on the value +9.89e-3 (radians). On the other hand, if
  28.  *   your horizon has hills such that your apparent horizon is, say, 1 degree
  29.  *   above sea level, you would allow for this by setting dis to -1.75e-2
  30.  *   (radians).
  31.  *
  32.  * algorithm:
  33.  *   the situation is described by two spherical triangles with two equal angles
  34.  *    (the right horizon intercepts, and the common horizon transverse):
  35.  *   given lat, d(=d1+d2), and dis find z(=z1+z2) and rho, where      /| eq pole
  36.  *     lat = latitude,                                              /  |
  37.  *     dis = horizon displacement (>0 is below ideal)             / rho|
  38.  *     d = angle from pole = PI/2 - declination                /       |
  39.  *     z = azimuth east of north                            /          |
  40.  *     rho = polar rotation from down = PI - hour angle    /           | 
  41.  *   solve simultaneous equations for d1 and d2:         /             |
  42.  *     1) cos(d) = cos(d1+d2)                           / d2           | lat
  43.  *            = cos(d1)cos(d2) - sin(d1)sin(d2)        /               |
  44.  *     2) sin(d2) = sin(lat)sin(d1)/sin(dis)          /                |
  45.  *   then can solve for z1, z2 and rho, taking       /                 |
  46.  *     care to preserve quadrant information.       /                 -|
  47.  *                                              z1 /        z2       | |
  48.  *                      ideal horizon ------------/--------------------| 
  49.  *                                         | |   /                     N
  50.  *                                          -|  / d1
  51.  *                                       dis | /
  52.  *                                           |/
  53.  *                  apparent horizon  ---------------------------------
  54.  *
  55.  * note that when lat=0 this all breaks down (because d2 and z2 degenerate to 0)
  56.  *   but fortunately then we can solve for z and rho directly.
  57.  *
  58.  * status: 0: normal; 1: never rises; -1: circumpolar; 2: trouble.
  59.  */
  60. void
  61. riset (ra, dec, lat, dis, lstr, lsts, azr, azs, status)
  62. double ra, dec;
  63. double lat, dis;
  64. double *lstr, *lsts;
  65. double *azr, *azs;
  66. int *status;
  67. {
  68. #define    EPS    (1e-6)    /* math rounding fudge - always the way, eh? */
  69.     double d;    /* angle from pole */
  70.     double h;    /* hour angle */
  71.     double crho;    /* cos hour-angle complement */
  72.     int shemi;    /* flag for southern hemisphere reflection */
  73.  
  74.     d = PI/2 - dec;
  75.  
  76.     /* reflect if in southern hemisphere.
  77.      * (then reflect azimuth back after computation.)
  78.      */
  79.     if ((shemi = (lat < 0)) != 0) {
  80.         lat = -lat;
  81.         d = PI - d;
  82.     }
  83.  
  84.     /* do the easy ones (and avoid violated assumptions) if d arc never
  85.      * meets horizon. 
  86.      */
  87.     if (d <= lat + dis + EPS) {
  88.         *status = -1; /* never sets */
  89.         return;
  90.     }
  91.     if (d >= PI - lat + dis - EPS) {
  92.         *status = 1; /* never rises */
  93.         return;
  94.     }
  95.  
  96.     /* find rising azimuth and cosine of hour-angle complement */
  97.     if (lat > EPS) {
  98.         double d2, d1; /* polr arc to ideal hzn, and corrctn for apparent */
  99.         double z2, z1; /* azimuth to ideal horizon, and " */
  100.         double a;       /* intermediate temp */
  101.         double sdis, slat, clat, cz2, cd2;    /* trig temps */
  102.         sdis = sin(dis);
  103.         slat = sin(lat);
  104.         a = sdis*sdis + slat*slat + 2*cos(d)*sdis*slat;
  105.         if (a <= 0) {
  106.         *status = 2; /* can't happen - hah! */
  107.         return;
  108.         }
  109.         d1 = asin (sin(d) * sdis / sqrt(a));
  110.         d2 = d - d1;
  111.         cd2 = cos(d2);
  112.         clat = cos(lat);
  113.         cz2 = cd2/clat;
  114.         z2 = acos (cz2);
  115.         z1 = acos (cos(d1)/cos(dis));
  116.         if (dis < 0)
  117.         z1 = -z1;
  118.         *azr = z1 + z2;
  119.         range (azr, PI);
  120.         crho = (cz2 - cd2*clat)/(sin(d2)*slat);
  121.     } else {
  122.         *azr = acos (cos(d)/cos(dis));
  123.         crho = sin(dis)/sin(d);
  124.     }
  125.  
  126.     if (shemi)
  127.         *azr = PI - *azr;
  128.         *azs = 2*PI - *azr;
  129.     
  130.     /* find hour angle */
  131.     h = PI - acos (crho);
  132.         *lstr = radhr(ra-h);
  133.     *lsts = radhr(ra+h);
  134.     range (lstr, 24.0);
  135.     range (lsts, 24.0);
  136.  
  137.     *status = 0;
  138. }
  139.