home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / astrnomy / ephem421.zip / EQ_ECL.C < prev    next >
C/C++ Source or Header  |  1990-09-13  |  2KB  |  66 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "astro.h"
  4.  
  5. #define    EQtoECL    1
  6. #define    ECLtoEQ    (-1)
  7.  
  8. /* given the modified Julian date, mjd, and an equitorial ra and dec, each in
  9.  * radians, find the corresponding geocentric ecliptic latitude, *lat, and
  10.  * longititude, *lng, also each in radians.
  11.  * correction for the effect on the angle of the obliquity due to nutation is
  12.  * included.
  13.  */
  14. eq_ecl (mjd, ra, dec, lat, lng)
  15. double mjd, ra, dec;
  16. double *lat, *lng;
  17. {
  18.     ecleq_aux (EQtoECL, mjd, ra, dec, lng, lat);
  19. }
  20.  
  21. /* given the modified Julian date, mjd, and a geocentric ecliptic latitude,
  22.  * *lat, and longititude, *lng, each in radians, find the corresponding
  23.  * equitorial ra and dec, also each in radians.
  24.  * correction for the effect on the angle of the obliquity due to nutation is
  25.  * included.
  26.  */
  27. ecl_eq (mjd, lat, lng, ra, dec)
  28. double mjd, lat, lng;
  29. double *ra, *dec;
  30. {
  31.     ecleq_aux (ECLtoEQ, mjd, lng, lat, ra, dec);
  32. }
  33.  
  34. static
  35. ecleq_aux (sw, mjd, x, y, p, q)
  36. int sw;            /* +1 for eq to ecliptic, -1 for vv. */
  37. double mjd, x, y;    /* sw==1: x==ra, y==dec.  sw==-1: x==lng, y==lat. */
  38. double *p, *q;        /* sw==1: p==lng, q==lat. sw==-1: p==ra, q==dec. */
  39. {
  40.     static double lastmjd = -10000;    /* last mjd calculated */
  41.     static double seps, ceps;    /* sin and cos of mean obliquity */
  42.     double sx, cx, sy, cy, ty;
  43.  
  44.     if (mjd != lastmjd) {
  45.         double eps;
  46.         double deps, dpsi;
  47.         obliquity (mjd, &eps);        /* mean obliquity for date */
  48.         nutation (mjd, &deps, &dpsi);
  49.         eps += deps;
  50.             seps = sin(eps);
  51.         ceps = cos(eps);
  52.         lastmjd = mjd;
  53.     }
  54.  
  55.     sy = sin(y);
  56.     cy = cos(y);                /* always non-negative */
  57.         if (fabs(cy)<1e-20) cy = 1e-20;        /* insure > 0 */
  58.         ty = sy/cy;
  59.     cx = cos(x);
  60.     sx = sin(x);
  61.         *q = asin((sy*ceps)-(cy*seps*sx*sw));
  62.         *p = atan(((sx*ceps)+(ty*seps*sw))/cx);
  63.         if (cx<0) *p += PI;        /* account for atan quad ambiguity */
  64.     range (p, 2*PI);
  65. }
  66.