home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume19 / xephem / part21 / eq_ecl.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-15  |  2.4 KB  |  84 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 nutation P_((double Mjd, double *deps, double *dpsi));
  12. extern void obliquity P_((double Mjd, double *eps));
  13. extern void range P_((double *v, double r));
  14.  
  15. void eq_ecl P_((double mjd, double ra, double dec, double *lat, double *lng));
  16. void ecl_eq P_((double mjd, double lat, double lng, double *ra, double *dec));
  17. static void ecleq_aux P_((int sw, double mjd, double x, double y, double *p, double *q));
  18.  
  19. #undef P_
  20.  
  21. #define    EQtoECL    1
  22. #define    ECLtoEQ    (-1)
  23.  
  24. /* given the modified Julian date, mjd, and an equitorial ra and dec, each in
  25.  * radians, find the corresponding geocentric ecliptic latitude, *lat, and
  26.  * longititude, *lng, also each in radians.
  27.  * correction for the effect on the angle of the obliquity due to nutation is
  28.  * included.
  29.  */
  30. void
  31. eq_ecl (mjd, ra, dec, lat, lng)
  32. double mjd, ra, dec;
  33. double *lat, *lng;
  34. {
  35.     ecleq_aux (EQtoECL, mjd, ra, dec, lng, lat);
  36. }
  37.  
  38. /* given the modified Julian date, mjd, and a geocentric ecliptic latitude,
  39.  * *lat, and longititude, *lng, each in radians, find the corresponding
  40.  * equitorial ra and dec, also each in radians.
  41.  * correction for the effect on the angle of the obliquity due to nutation is
  42.  * included.
  43.  */
  44. void
  45. ecl_eq (mjd, lat, lng, ra, dec)
  46. double mjd, lat, lng;
  47. double *ra, *dec;
  48. {
  49.     ecleq_aux (ECLtoEQ, mjd, lng, lat, ra, dec);
  50. }
  51.  
  52. static void
  53. ecleq_aux (sw, mjd, x, y, p, q)
  54. int sw;            /* +1 for eq to ecliptic, -1 for vv. */
  55. double mjd, x, y;    /* sw==1: x==ra, y==dec.  sw==-1: x==lng, y==lat. */
  56. double *p, *q;        /* sw==1: p==lng, q==lat. sw==-1: p==ra, q==dec. */
  57. {
  58.     static double lastmjd = -10000;    /* last mjd calculated */
  59.     static double seps, ceps;    /* sin and cos of mean obliquity */
  60.     double sx, cx, sy, cy, ty;
  61.  
  62.     if (mjd != lastmjd) {
  63.         double eps;
  64.         double deps, dpsi;
  65.         obliquity (mjd, &eps);        /* mean obliquity for date */
  66.         nutation (mjd, &deps, &dpsi);
  67.         eps += deps;
  68.             seps = sin(eps);
  69.         ceps = cos(eps);
  70.         lastmjd = mjd;
  71.     }
  72.  
  73.     sy = sin(y);
  74.     cy = cos(y);                /* always non-negative */
  75.         if (fabs(cy)<1e-20) cy = 1e-20;        /* insure > 0 */
  76.         ty = sy/cy;
  77.     cx = cos(x);
  78.     sx = sin(x);
  79.         *q = asin((sy*ceps)-(cy*seps*sx*sw));
  80.         *p = atan(((sx*ceps)+(ty*seps*sw))/cx);
  81.         if (cx<0) *p += PI;        /* account for atan quad ambiguity */
  82.     range (p, 2*PI);
  83. }
  84.