source: Sophya/trunk/SophyaExt/XephemAstroLib/helio.c@ 2643

#include <stdlib.h>
#include <math.h>
#include <math.h>

#include "astro.h"

/* given geocentric time jd and coords of a distant object at ra/dec (J2000),
 * find the difference in time between light arriving at earth vs the sun.
 * *hcp must be subtracted from geocentric jd to get heliocentric jd.
 * From RLM Oct 12, 1995.
 */
void
heliocorr (double jd, double ra, double dec, double *hcp)
{
        double e;       /* obliquity of ecliptic */
        double n;       /* day number */
        double g;       /* solar mean anomaly */
        double L;       /* solar ecliptic longitude */
        double l;       /* mean solar ecliptic longitude */
        double R;       /* sun distance, AU */
        double X, Y;    /* equatorial rectangular solar coords */
        double cdec, sdec;
        double cra, sra;

        /* following algorithm really wants EOD */
        precess (J2000, jd - MJD0, &ra, &dec);

        cdec = cos(dec);
        sdec = sin(dec);
        cra = cos(ra);
        sra = sin(ra);

        n = jd - 2451545.0;     /* use epoch 2000 */
        e = degrad(23.439 - 0.0000004*n);
        g = degrad(357.528) + degrad(0.9856003)*n;
        L = degrad(280.461) + degrad(0.9856474)*n;
        l = L + degrad(1.915)*sin(g) + degrad(0.02)*sin(2.0*g);
        R = 1.00014 - 0.01671*cos(g) - 0.00014*cos(2.0*g);
        X = R*cos(l);
        Y = R*cos(e)*sin(l);

#if 0
            printf ("n=%g g=%g L=%g l=%g R=%g X=%g Y=%g\n",
                                n, raddeg(g), raddeg(L), raddeg(l), R, X, Y);
#endif

        *hcp = 0.0057755 * (cdec*cra*X + (cdec*sra + tan(e)*sdec)*Y);
}

/* For RCS Only -- Do Not Edit */
static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: helio.c,v $ $Date: 2005-01-17 10:13:05 $ $Revision: 1.4 $ $Name: not supported by cvs2svn $"};