source: Sophya/trunk/SophyaExt/XephemAstroLib/eq_ecl.c@ 3237

#include <stdio.h>
#include <math.h>

#include "astro.h"

static void ecleq_aux (int sw, double mj, double x, double y,
    double *p, double *q);

#define EQtoECL 1
#define ECLtoEQ (-1)


/* given the modified Julian date, mj, and an equitorial ra and dec, each in
 * radians, find the corresponding geocentric ecliptic latitude, *lt, and
 * longititude, *lg, also each in radians.
 * correction for the effect on the angle of the obliquity due to nutation is
 * not included.
 */
void
eq_ecl (double mj, double ra, double dec, double *lt, double *lg)
{
        ecleq_aux (EQtoECL, mj, ra, dec, lg, lt);
}

/* given the modified Julian date, mj, and a geocentric ecliptic latitude,
 * *lt, and longititude, *lg, each in radians, find the corresponding
 * equitorial ra and dec, also each in radians.
 * correction for the effect on the angle of the obliquity due to nutation is
 * not included.
 */
void
ecl_eq (double mj, double lt, double lg, double *ra, double *dec)
{
        ecleq_aux (ECLtoEQ, mj, lg, lt, ra, dec);
}

static void
ecleq_aux (
int sw,                 /* +1 for eq to ecliptic, -1 for vv. */
double mj,
double x, double y,     /* sw==1: x==ra, y==dec.  sw==-1: x==lg, y==lt. */
double *p, double *q)   /* sw==1: p==lg, q==lt. sw==-1: p==ra, q==dec. */
{
        static double lastmj = -10000;  /* last mj calculated */
        static double seps, ceps;       /* sin and cos of mean obliquity */
        double sx, cx, sy, cy, ty, sq;

        if (mj != lastmj) {
            double eps;
            obliquity (mj, &eps);               /* mean obliquity for date */
            seps = sin(eps);
            ceps = cos(eps);
            lastmj = mj;
        }

        sy = sin(y);
        cy = cos(y);                            /* always non-negative */
        if (fabs(cy)<1e-20) cy = 1e-20;         /* insure > 0 */
        ty = sy/cy;
        cx = cos(x);
        sx = sin(x);
        sq = (sy*ceps)-(cy*seps*sx*sw);
        if (sq < -1) sq = -1;
        if (sq >  1) sq =  1;
        *q = asin(sq);
        *p = atan(((sx*ceps)+(ty*seps*sw))/cx);
        if (cx<0) *p += PI;             /* account for atan quad ambiguity */
        range (p, 2*PI);
}

/* For RCS Only -- Do Not Edit */
static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: eq_ecl.c,v $ $Date: 2006-11-22 13:53:29 $ $Revision: 1.6 $ $Name: not supported by cvs2svn $"};