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