home *** CD-ROM | disk | FTP | other *** search
/ rtsi.com / 2014.01.www.rtsi.com.tar / www.rtsi.com / OS9 / OSK / EFFO / pd8.lzh / SRC / riset.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-04-13  |  2.3 KB  |  85 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, in radians, the observer's
  6.  *   latitude in radians, and a horizon displacement correction, dis, in radians
  7.  *   find the local sidereal times and azimuths of rising and setting, lstr/s
  8.  *   and azr/s, 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.  * status: 0: normal; 1: never rises; -1: circumpolar; 2: trouble.
  19.  */
  20. riset (ra, dec, lat, dis, lstr, lsts, azr, azs, status)
  21. double ra, dec;
  22. double lat, dis;
  23. double *lstr, *lsts;
  24. double *azr, *azs;
  25. int *status;
  26. {
  27.     static double lastlat = 0, slat = 0, clat = 1.0;
  28.     double dec1, sdec, cdec, tdec;
  29.     double psi, spsi, cpsi;
  30.     double h, dh, ch;    /* hr angle, delta and cos */
  31.  
  32.     /* avoid needless sin/cos since latitude doesn't change often */
  33.         if (lat != lastlat) {
  34.         clat = cos(lat);
  35.         slat = sin(lat);
  36.         lastlat = lat;
  37.     }
  38.  
  39.     /* can't cope with objects very near the celestial pole nor if we 
  40.      * are located very near the earth's poles.
  41.      */
  42.     cdec = cos(dec);
  43.         if (fabs(cdec*clat) < 1e-10) {
  44.         /* trouble */
  45.         *status = 2;
  46.         return;
  47.     }
  48.  
  49.         cpsi = slat/cdec;
  50.         if (cpsi>1.0) cpsi = 1.0;
  51.         else if (cpsi<-1.0) cpsi = -1.0;
  52.         psi = acos(cpsi);
  53.     spsi = sin(psi);
  54.  
  55.         dh = dis*spsi;
  56.     dec1 = dec + (dis*cpsi);
  57.         sdec = sin(dec1);
  58.     tdec = tan(dec1);
  59.  
  60.         ch = (-1*slat*tdec)/clat;
  61.  
  62.         if (ch < -1) {
  63.         /* circumpolar */
  64.         *status = -1;
  65.         return;
  66.     }
  67.         if (ch > 1) {
  68.         /* never rises */
  69.         *status = 1;
  70.         return;
  71.     }
  72.  
  73.         *status = 0;
  74.     h = acos(ch)+dh;
  75.  
  76.         *lstr = 24+radhr(ra-h);
  77.     *lsts = radhr(ra+h);
  78.     range (lstr, 24.0);
  79.     range (lsts, 24.0);
  80.  
  81.     *azr = acos(sdec/clat);
  82.     range (azr, 2*PI);
  83.         *azs = 2*PI - *azr;
  84. }
  85.