source: Sophya/trunk/SophyaExt/XephemAstroLib/jupmoon.c@ 3111

/* jupiter moon info */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include <math.h>

#include "astro.h"
#include "bdl.h"

static int use_bdl (double jd, char *dir, MoonData md[J_NMOONS]);
static void meeus_jupiter (double d, double *cmlI, double *cmlII,
    MoonData md[J_NMOONS]);
static void moonradec (double jupsize, MoonData md[J_NMOONS]);
static void moonSVis (Obj *sop, Obj *jop, MoonData md[J_NMOONS]);
static void moonEVis (MoonData md[J_NMOONS]);
static void moonPShad (Obj *sop, Obj *jop, MoonData md[J_NMOONS]);
static void moonTrans (MoonData md[J_NMOONS]);

/* moon table and a few other goodies and when it was last computed */
static double mdmjd = -123456;
static MoonData jmd[J_NMOONS] = {
    {"Jupiter", NULL},
    {"Io", "I"},
    {"Europa", "II"},
    {"Ganymede", "III"},
    {"Callisto", "IV"}
};
static double sizemjd;  /* size at last mjd */
static double cmlImjd;  /* central meridian long sys I, at last mjd */
static double cmlIImjd; /*    "                      II      " */

/* These values are from the Explanatory Supplement.
 * Precession degrades them gradually over time.
 */
#define POLE_RA   degrad(268.05)         /* RA of Jupiter's north pole */
#define POLE_DEC  degrad(64.50)         /* Dec of Jupiter's north pole */


/* get jupiter info in md[0], moon info in md[1..J_NMOONS-1].
 * if !dir always use meeus model.
 * if !jop caller just wants md[] for names
 * N.B. we assume sop and jop are updated.
 */
void
jupiter_data (
double Mjd,             /* mjd */
char dir[],             /* dir in which to look for helper files */
Obj *sop,               /* Sun */
Obj *jop,               /* jupiter */
double *sizep,          /* jup angular diam, rads */
double *cmlI, double *cmlII,            /* central meridian longitude, rads */
double *polera, double *poledec,        /* pole location */
MoonData md[J_NMOONS])  /* return info */
{
        double JD;

        /* always copy back at least for name */
        memcpy (md, jmd, sizeof(jmd));

        /* pole */
        if (polera) *polera = POLE_RA;
        if (poledec) *poledec = POLE_DEC;

        /* nothing else if repeat call or just want names */
        if (Mjd == mdmjd || !jop) {
            if (jop) {
                *sizep = sizemjd;
                *cmlI = cmlImjd;
                *cmlII = cmlIImjd;
            }
            return;
        }
        JD = Mjd + MJD0;

        /* planet in [0] */
        md[0].ra = jop->s_ra;
        md[0].dec = jop->s_dec;
        md[0].mag = get_mag(jop);
        md[0].x = 0;
        md[0].y = 0;
        md[0].z = 0;
        md[0].evis = 1;
        md[0].svis = 1;

        /* size is straight from jop */
        *sizep = degrad(jop->s_size/3600.0);

        /* mags from JPL ephemeris */
        md[1].mag = 5.7;
        md[2].mag = 5.8;
        md[3].mag = 5.3;
        md[4].mag = 6.7;

        /* get moon data from BDL if possible, else Meeus' model.
         * always use Meeus for cml
         */
        if (dir && use_bdl (JD, dir, md) == 0)
            meeus_jupiter (Mjd, cmlI, cmlII, NULL);
        else
            meeus_jupiter (Mjd, cmlI, cmlII, md);

        /* set visibilities */
        moonSVis (sop, jop, md);
        moonPShad (sop, jop, md);
        moonEVis (md);
        moonTrans (md);

        /* fill in moon ra and dec */
        moonradec (*sizep, md);

        /* save */
        mdmjd = Mjd;
        sizemjd = *sizep;
        cmlImjd = *cmlI;
        cmlIImjd = *cmlII;
        memcpy (jmd, md, sizeof(jmd));
}

/* hunt for BDL file in dir[] and use if possible
 * return 0 if ok, else -1
 */
static int
use_bdl (
double JD,              /* julian date */
char dir[],             /* directory */
MoonData md[J_NMOONS])  /* fill md[1..NM-1].x/y/z for each moon */
{
#define JUPRAU  .0004769108     /* jupiter radius, AU */
        double x[J_NMOONS], y[J_NMOONS], z[J_NMOONS];
        char buf[1024];
        FILE *fp;
        char *fn;
        int i;

        /* check ranges and appropriate data file */
        if (JD < 2451179.50000)         /* Jan 1 1999 UTC */
            return (-1);
        if (JD < 2455562.5)             /* Jan 1 2011 UTC */
            fn = "jupiter.9910";
        else if (JD < 2459215.5)        /* Jan 1 2021 UTC */
            fn = "jupiter.1020";
        else
            return (-1);

        /* open */
        (void) sprintf (buf, "%s/%s", dir, fn);
        fp = fopen (buf, "r");
        if (!fp) {
            fprintf (stderr, "%s: %s\n", fn, strerror(errno));
            return (-1);
        }

        /* use it */
        if ((i = read_bdl (fp, JD, x, y, z, buf)) < 0) {
            fprintf (stderr, "%s: %s\n", fn, buf);
            fclose (fp);
            return (-1);
        }
        if (i != J_NMOONS-1) {
            fprintf (stderr, "%s: BDL says %d moons, code expects %d", fn, 
                                                                i, J_NMOONS-1);
            fclose (fp);
            return (-1);
        }

        /* copy into md[1..NM-1] with our scale and sign conventions */
        for (i = 1; i < J_NMOONS; i++) {
            md[i].x =  x[i-1]/JUPRAU;   /* we want jup radii +E */
            md[i].y = -y[i-1]/JUPRAU;   /* we want jup radii +S */
            md[i].z = -z[i-1]/JUPRAU;   /* we want jup radii +front */
        }

        /* ok */
        fclose (fp);
        return (0);
}

/* compute location of GRS and Galilean moons.
 * if md == NULL, just to cml.
 * from "Astronomical Formulae for Calculators", 2nd ed, by Jean Meeus,
 *   Willmann-Bell, Richmond, Va., U.S.A. (c) 1982, chapters 35 and 36.
 */
static void
meeus_jupiter(
double d,
double *cmlI, double *cmlII,    /* central meridian longitude, rads */
MoonData md[J_NMOONS])  /* fill in md[1..NM-1].x/y/z for each moon.
                         * N.B. md[0].ra/dec must already be set
                         */
{
#define dsin(x) sin(degrad(x))
#define dcos(x) cos(degrad(x))
        double A, B, Del, J, K, M, N, R, V;
        double cor_u1, cor_u2, cor_u3, cor_u4;
        double solc, tmp, G, H, psi, r, r1, r2, r3, r4;
        double u1, u2, u3, u4;
        double lam, Ds;
        double z1, z2, z3,  z4;
        double De, dsinDe;
        double theta, phi;
        double tvc, pvc;
        double salpha, calpha;
        int i;

        V = 134.63 + 0.00111587 * d;

        M = (358.47583 + 0.98560003*d);
        N = (225.32833 + 0.0830853*d) + 0.33 * dsin (V);

        J = 221.647 + 0.9025179*d - 0.33 * dsin(V);

        A = 1.916*dsin(M) + 0.02*dsin(2*M);
        B = 5.552*dsin(N) + 0.167*dsin(2*N);
        K = (J+A-B);
        R = 1.00014 - 0.01672 * dcos(M) - 0.00014 * dcos(2*M);
        r = 5.20867 - 0.25192 * dcos(N) - 0.00610 * dcos(2*N);
        Del = sqrt (R*R + r*r - 2*R*r*dcos(K));
        psi = raddeg (asin (R/Del*dsin(K)));

        *cmlI  = degrad(268.28 + 877.8169088*(d - Del/173) + psi - B);
        range (cmlI, 2*PI);
        *cmlII = degrad(290.28 + 870.1869088*(d - Del/173) + psi - B);
        range (cmlII, 2*PI);

        /* that's it if don't want moon info too */
        if (!md)
            return;

        solc = (d - Del/173.);  /* speed of light correction */
        tmp = psi - B;

        u1 = 84.5506 + 203.4058630 * solc + tmp;
        u2 = 41.5015 + 101.2916323 * solc + tmp;
        u3 = 109.9770 + 50.2345169 * solc + tmp;
        u4 = 176.3586 + 21.4879802 * solc + tmp;

        G = 187.3 + 50.310674 * solc;
        H = 311.1 + 21.569229 * solc;
      
        cor_u1 =  0.472 * dsin (2*(u1-u2));
        cor_u2 =  1.073 * dsin (2*(u2-u3));
        cor_u3 =  0.174 * dsin (G);
        cor_u4 =  0.845 * dsin (H);
      
        r1 = 5.9061 - 0.0244 * dcos (2*(u1-u2));
        r2 = 9.3972 - 0.0889 * dcos (2*(u2-u3));
        r3 = 14.9894 - 0.0227 * dcos (G);
        r4 = 26.3649 - 0.1944 * dcos (H);

        md[1].x = -r1 * dsin (u1+cor_u1);
        md[2].x = -r2 * dsin (u2+cor_u2);
        md[3].x = -r3 * dsin (u3+cor_u3);
        md[4].x = -r4 * dsin (u4+cor_u4);

        lam = 238.05 + 0.083091*d + 0.33*dsin(V) + B;
        Ds = 3.07*dsin(lam + 44.5);
        De = Ds - 2.15*dsin(psi)*dcos(lam+24.)
                - 1.31*(r-Del)/Del*dsin(lam-99.4);
        dsinDe = dsin(De);

        z1 = r1 * dcos(u1+cor_u1);
        z2 = r2 * dcos(u2+cor_u2);
        z3 = r3 * dcos(u3+cor_u3);
        z4 = r4 * dcos(u4+cor_u4);

        md[1].y = z1*dsinDe;
        md[2].y = z2*dsinDe;
        md[3].y = z3*dsinDe;
        md[4].y = z4*dsinDe;

        /* compute sky transformation angle as triple vector product */
        tvc = PI/2.0 - md[0].dec;
        pvc = md[0].ra;
        theta = PI/2.0 - POLE_DEC;
        phi = POLE_RA;
        salpha = -sin(tvc)*sin(theta)*(cos(pvc)*sin(phi) - sin(pvc)*cos(phi));
        calpha = sqrt (1.0 - salpha*salpha);

        for (i = 0; i < J_NMOONS; i++) {
            double tx =  md[i].x*calpha + md[i].y*salpha;
            double ty = -md[i].x*salpha + md[i].y*calpha;
            md[i].x = tx;
            md[i].y = ty;
        }

        md[1].z = z1;
        md[2].z = z2;
        md[3].z = z3;
        md[4].z = z4;
}


/* given jupiter loc in md[0].ra/dec and size, and location of each moon in 
 * md[1..NM-1].x/y in jup radii, find ra/dec of each moon in md[1..NM-1].ra/dec.
 */
static void
moonradec (
double jupsize,         /* jup diameter, rads */
MoonData md[J_NMOONS])  /* fill in RA and Dec */
{
        double juprad = jupsize/2;
        double jupra = md[0].ra;
        double jupdec = md[0].dec;
        int i;

        for (i = 1; i < J_NMOONS; i++) {
            double dra  = juprad * md[i].x;
            double ddec = juprad * md[i].y;
            md[i].ra  = jupra + dra;
            md[i].dec = jupdec - ddec;
        }
}

/* set svis according to whether moon is in sun light */
static void
moonSVis(
Obj *sop,               /* SUN */
Obj *jop,               /* jupiter */
MoonData md[J_NMOONS])
{
        double esd = sop->s_edist;
        double eod = jop->s_edist;
        double sod = jop->s_sdist;
        double soa = degrad(jop->s_elong);
        double esa = asin(esd*sin(soa)/sod);
        double   h = sod*jop->s_hlat;
        double nod = h*(1./eod - 1./sod);
        double sca = cos(esa), ssa = sin(esa);
        int i;

        for (i = 1; i < J_NMOONS; i++) {
            MoonData *mdp = &md[i];
            double xp =  sca*mdp->x + ssa*mdp->z;
            double yp =  mdp->y;
            double zp = -ssa*mdp->x + sca*mdp->z;
            double ca = cos(nod), sa = sin(nod);
            double xpp = xp;
            double ypp = ca*yp - sa*zp;
            double zpp = sa*yp + ca*zp;
            int outside = xpp*xpp + ypp*ypp > 1.0;
            int infront = zpp > 0.0;
            mdp->svis = outside || infront;
        }
}

/* set evis according to whether moon is geometrically visible from earth */
static void
moonEVis (MoonData md[J_NMOONS])
{
        int i;

        for (i = 1; i < J_NMOONS; i++) {
            MoonData *mdp = &md[i];
            int outside = mdp->x*mdp->x + mdp->y*mdp->y > 1.0;
            int infront = mdp->z > 0.0;
            mdp->evis = outside || infront;
        }
}

/* set pshad and sx,sy shadow info */
static void
moonPShad(
Obj *sop,               /* SUN */
Obj *jop,               /* jupiter */
MoonData md[J_NMOONS])
{
        int i;

        for (i = 1; i < J_NMOONS; i++) {
            MoonData *mdp = &md[i];
            mdp->pshad = !plshadow (jop, sop, POLE_RA, POLE_DEC, mdp->x,
                                            mdp->y, mdp->z, &mdp->sx, &mdp->sy);
        }
}

/* set whether moons are transiting */
static void
moonTrans (MoonData md[J_NMOONS])
{
        int i;

        for (i = 1; i < J_NMOONS; i++) {
            MoonData *mdp = &md[i];
            mdp->trans = mdp->z > 0 && mdp->x*mdp->x + mdp->y*mdp->y < 1;
        }
}

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