home *** CD-ROM | disk | FTP | other *** search
/ Collection of Education / collectionofeducationcarat1997.iso / SCIENCE / EPHEM421.ZIP / AA_HADEC.C next >
C/C++ Source or Header  |  1990-09-13  |  2KB  |  84 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "astro.h"
  4.  
  5. /* given latitude (n+, radians), lat, altitude (up+, radians), alt, and
  6.  * azimuth (angle round to the east from north+, radians),
  7.  * return hour angle (radians), ha, and declination (radians), dec.
  8.  */
  9. aa_hadec (lat, alt, az, ha, dec)
  10. double lat;
  11. double alt, az;
  12. double *ha, *dec;
  13. {
  14.     aaha_aux (lat, az, alt, ha, dec);
  15. }
  16.  
  17. /* given latitude (n+, radians), lat, hour angle (radians), ha, and declination
  18.  * (radians), dec,
  19.  * return altitude (up+, radians), alt, and
  20.  * azimuth (angle round to the east from north+, radians),
  21.  */
  22. hadec_aa (lat, ha, dec, alt, az)
  23. double lat;
  24. double ha, dec;
  25. double *alt, *az;
  26. {
  27.     aaha_aux (lat, ha, dec, az, alt);
  28. }
  29.  
  30. /* the actual formula is the same for both transformation directions so
  31.  * do it here once for each way.
  32.  * N.B. all arguments are in radians.
  33.  */
  34. static
  35. aaha_aux (lat, x, y, p, q)
  36. double lat;
  37. double x, y;
  38. double *p, *q;
  39. {
  40.     static double lastlat = -1000.;
  41.     static double sinlastlat, coslastlat;
  42.     double sy, cy;
  43.     double sx, cx;
  44.     double sq, cq;
  45.     double a;
  46.     double cp;
  47.  
  48.     /* latitude doesn't change much, so try to reuse the sin and cos evals.
  49.      */
  50.     if (lat != lastlat) {
  51.         sinlastlat = sin (lat);
  52.         coslastlat = cos (lat);
  53.         lastlat = lat;
  54.     }
  55.  
  56.     sy = sin (y);
  57.     cy = cos (y);
  58.     sx = sin (x);
  59.     cx = cos (x);
  60.  
  61. /* define GOODATAN2 if atan2 returns full range -PI through +PI.
  62.  */
  63. #ifdef GOODATAN2
  64.     *q = asin ((sy*sinlastlat) + (cy*coslastlat*cx));
  65.     *p = atan2 (-cy*sx, -cy*cx*sinlastlat + sy*coslastlat);
  66. #else
  67. #define    EPS    (1e-20)
  68.     sq = (sy*sinlastlat) + (cy*coslastlat*cx);
  69.     *q = asin (sq);
  70.     cq = cos (*q);
  71.     a = coslastlat*cq;
  72.     if (a > -EPS && a < EPS)
  73.         a = a < 0 ? -EPS : EPS; /* avoid / 0 */
  74.     cp = (sy - (sinlastlat*sq))/a;
  75.     if (cp >= 1.0)    /* the /a can be slightly > 1 */
  76.         *p = 0.0;
  77.     else if (cp <= -1.0)
  78.         *p = PI;
  79.     else
  80.         *p = acos ((sy - (sinlastlat*sq))/a);
  81.     if (sx>0) *p = 2.0*PI - *p;
  82. #endif
  83. }
  84.