home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume19 / xephem / part21 / aa_hadec.c next >
Encoding:
C/C++ Source or Header  |  1993-05-15  |  2.2 KB  |  98 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. void aa_hadec P_((double lat, double alt, double az, double *ha, double *dec));
  12. void hadec_aa P_((double lat, double ha, double dec, double *alt, double *az));
  13. static void aaha_aux P_((double lat, double x, double y, double *p, double *q));
  14.  
  15. #undef P_
  16.  
  17. /* given latitude (n+, radians), lat, altitude (up+, radians), alt, and
  18.  * azimuth (angle round to the east from north+, radians),
  19.  * return hour angle (radians), ha, and declination (radians), dec.
  20.  */
  21. void
  22. aa_hadec (lat, alt, az, ha, dec)
  23. double lat;
  24. double alt, az;
  25. double *ha, *dec;
  26. {
  27.     aaha_aux (lat, az, alt, ha, dec);
  28. }
  29.  
  30. /* given latitude (n+, radians), lat, hour angle (radians), ha, and declination
  31.  * (radians), dec,
  32.  * return altitude (up+, radians), alt, and
  33.  * azimuth (angle round to the east from north+, radians),
  34.  */
  35. void
  36. hadec_aa (lat, ha, dec, alt, az)
  37. double lat;
  38. double ha, dec;
  39. double *alt, *az;
  40. {
  41.     aaha_aux (lat, ha, dec, az, alt);
  42. }
  43.  
  44. /* the actual formula is the same for both transformation directions so
  45.  * do it here once for each way.
  46.  * N.B. all arguments are in radians.
  47.  */
  48. static void
  49. aaha_aux (lat, x, y, p, q)
  50. double lat;
  51. double x, y;
  52. double *p, *q;
  53. {
  54.     static double lastlat = -1000.;
  55.     static double sinlastlat, coslastlat;
  56.     double sy, cy;
  57.     double sx, cx;
  58.     double sq, cq;
  59.     double a;
  60.     double cp;
  61.  
  62.     /* latitude doesn't change much, so try to reuse the sin and cos evals.
  63.      */
  64.     if (lat != lastlat) {
  65.         sinlastlat = sin (lat);
  66.         coslastlat = cos (lat);
  67.         lastlat = lat;
  68.     }
  69.  
  70.     sy = sin (y);
  71.     cy = cos (y);
  72.     sx = sin (x);
  73.     cx = cos (x);
  74.  
  75. /* define GOODATAN2 if atan2 returns full range -PI through +PI.
  76.  */
  77. #ifdef GOODATAN2
  78.     *q = asin ((sy*sinlastlat) + (cy*coslastlat*cx));
  79.     *p = atan2 (-cy*sx, -cy*cx*sinlastlat + sy*coslastlat);
  80. #else
  81. #define    EPS    (1e-20)
  82.     sq = (sy*sinlastlat) + (cy*coslastlat*cx);
  83.     *q = asin (sq);
  84.     cq = cos (*q);
  85.     a = coslastlat*cq;
  86.     if (a > -EPS && a < EPS)
  87.         a = a < 0 ? -EPS : EPS; /* avoid / 0 */
  88.     cp = (sy - (sinlastlat*sq))/a;
  89.     if (cp >= 1.0)    /* the /a can be slightly > 1 */
  90.         *p = 0.0;
  91.     else if (cp <= -1.0)
  92.         *p = PI;
  93.     else
  94.         *p = acos ((sy - (sinlastlat*sq))/a);
  95.     if (sx>0) *p = 2.0*PI - *p;
  96. #endif
  97. }
  98.