source: Sophya/trunk/SophyaExt/XephemAstroLib/circum.c@ 3885

/* given a Now and an Obj with the object definition portion filled in,
 * fill in the sky position (s_*) portions.
 * calculation of positional coordinates reworked by
 * Michael Sternberg <sternberg@physik.tu-chemnitz.de>
 *  3/11/98: deflect was using op->s_hlong before being set in cir_pos().
 *  4/19/98: just edit a comment
 */

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

#include "astro.h"
#include "preferences.h"


static int obj_planet (Now *np, Obj *op);
static int obj_binary (Now *np, Obj *op);
static int obj_2binary (Now *np, Obj *op);
static int obj_fixed (Now *np, Obj *op);
static int obj_elliptical (Now *np, Obj *op);
static int obj_hyperbolic (Now *np, Obj *op);
static int obj_parabolic (Now *np, Obj *op);
static int sun_cir (Now *np, Obj *op);
static int moon_cir (Now *np, Obj *op);
static double solveKepler (double M, double e);
static void binaryStarOrbit (double t, double T, double e, double o, double O,
    double i, double a, double P, double *thetap, double *rhop);
static void cir_sky (Now *np, double lpd, double psi, double rp, double *rho,
    double lam, double bet, double lsn, double rsn, Obj *op);
static void cir_pos (Now *np, double bet, double lam, double *rho, Obj *op);
static void elongation (double lam, double bet, double lsn, double *el);
static void deflect (double mjd1, double lpd, double psi, double rsn,
    double lsn, double rho, double *ra, double *dec);
static double h_albsize (double H);

/* given a Now and an Obj, fill in the approprirate s_* fields within Obj.
 * return 0 if all ok, else -1.
 */
int
obj_cir (Now *np, Obj *op)
{
        op->o_flags &= ~NOCIRCUM;
        switch (op->o_type) {
        case BINARYSTAR: return (obj_binary (np, op));
        case FIXED:      return (obj_fixed (np, op));
        case ELLIPTICAL: return (obj_elliptical (np, op));
        case HYPERBOLIC: return (obj_hyperbolic (np, op));
        case PARABOLIC:  return (obj_parabolic (np, op));
        case EARTHSAT:   return (obj_earthsat (np, op));
        case PLANET:     return (obj_planet (np, op));
        default:
            printf ("obj_cir() called with type %d %s\n", op->o_type, op->o_name);
            abort();
            return (-1);        /* just for lint */
        }
}

static int
obj_planet (Now *np, Obj *op)
{
        double lsn, rsn;        /* true geoc lng of sun; dist from sn to earth*/
        double lpd, psi;        /* heliocentric ecliptic long and lat */
        double rp;              /* dist from sun */
        double rho;             /* dist from earth */
        double lam, bet;        /* geocentric ecliptic long and lat */
        double dia, mag;        /* angular diameter at 1 AU and magnitude */
        PLCode p;

        /* validate code and check for a few special cases */
        p = op->pl_code;
        if (p == SUN)
            return (sun_cir (np, op));
        if (p == MOON)
            return (moon_cir (np, op));
        if (op->pl_moon != X_PLANET)
            return (plmoon_cir (np, op));
        if (p < 0 || p > MOON) {
            printf ("unknown planet code: %d\n", p);
            abort();
        }

        /* planet itself */

        /* find solar ecliptical longitude and distance to sun from earth */
        sunpos (mjed, &lsn, &rsn, 0);

        /* find helio long/lat; sun/planet and earth/planet dist; ecliptic
         * long/lat; diameter and mag.
         */
        plans(mjed, p, &lpd, &psi, &rp, &rho, &lam, &bet, &dia, &mag);

        /* fill in all of op->s_* stuff except s_size and s_mag */
        cir_sky (np, lpd, psi, rp, &rho, lam, bet, lsn, rsn, op);

        /* set magnitude and angular size */
        set_smag (op, mag);
        op->s_size = (float)(dia/rho);

        return (0);
}

static int
obj_binary (Now *np, Obj *op)
{
        /* always compute circumstances of primary */
        if (obj_fixed (np, op) < 0)
            return (0);

        /* compute secondary only if requested, and always reset request flag */
        if (!op->b_2compute)
            return (0);
        op->b_2compute = 0;
        return (obj_2binary (np, op));
}

/* compute position of secondary component of a BINARYSTAR */
static int
obj_2binary (Now *np, Obj *op)
{
        if (op->b_nbp > 0) {
            /* we just have discrete pa/sep, project each from primary */
            int i;
            for (i = 0; i < op->b_nbp; i++) {
                BinPos *bp = &op->b_bp[i];
                bp->bp_dec = op->s_dec + bp->bp_sep*cos(bp->bp_pa);
                bp->bp_ra = op->s_ra + bp->bp_sep*sin(bp->bp_pa)/cos(op->s_dec);
            }
        } else {
            BinOrbit *bp = &op->b_bo;
            double t, theta, rho;

            mjd_year (mjd, &t);
            binaryStarOrbit (t, bp->bo_T, bp->bo_e, bp->bo_o, bp->bo_O,
                                bp->bo_i, bp->bo_a, bp->bo_P, &theta, &rho);
            bp->bo_pa = (float)theta;
            bp->bo_sep = (float)rho;
            rho = degrad(rho/3600.);    /* arc secs to rads */
            bp->bo_dec = op->s_dec + rho*cos(theta);
            bp->bo_ra =  op->s_ra  + rho*sin(theta)/cos(op->s_dec);
        }

        return (0);
}

/* from W. M. Smart */
static void
binaryStarOrbit (
double t,               /* desired ephemeris epoch, year */
double T,               /* epoch of periastron, year */
double e,               /* eccentricity */
double o,               /* argument of periastron, degrees */
double O,               /* ascending node, degrees */
double i,               /* inclination, degrees */
double a,               /* semi major axis, arcsecs */
double P,               /* period, years */
double *thetap,         /* position angle, rads E of N */
double *rhop)           /* separation, arcsecs */
{
        double M, E, cosE, nu, cosnu, r, rho, theta;

        /* find mean anomaly, insure 0..2*PI */
        M = 2*PI/P*(t-T);
        range (&M, 2*PI);

        /* solve for eccentric anomaly */
        E = solveKepler (M, e);
        cosE = cos(E);

        /* find true anomaly and separation */
        cosnu = (cosE - e)/(1.0 - e*cosE);
        r = a*(1.0 - e*e)/(1.0 + e*cosnu);
        nu = acos(cosnu);
        if (E > PI)
            nu = -nu;

        /* project onto sky */
        theta = atan(tan(nu+degrad(o))*cos(degrad(i))) + degrad(O);
        rho = r*cos(nu+degrad(o))/cos(theta-degrad(O));
        if (rho < 0) {
            theta += PI;
            rho = -rho;
        }
        range (&theta, 2*PI);

        *thetap = theta;
        *rhop = rho;
}

/* solve kepler equation using Newton-Raphson search.
 * Charles and Tatum have shown it always converges starting with PI.
 */
static double
solveKepler (double M, double e)
{
        double E, Eprime = PI;

        do {
            double cosE = cos(Eprime);
            E = Eprime;
            Eprime = (M - e*(E*cosE - sin(E)))/(1.0 - e*cosE);
        } while (fabs(E-Eprime) > 1e-7);

        return (Eprime);
}

static int
obj_fixed (Now *np, Obj *op)
{
        double lsn, rsn;        /* true geoc lng of sun, dist from sn to earth*/
        double lam, bet;        /* geocentric ecliptic long and lat */
        double ha;              /* local hour angle */
        double el;              /* elongation */
        double alt, az;         /* current alt, az */
        double ra, dec;         /* ra and dec at equinox of date */
        double rpm, dpm;        /* astrometric ra and dec with PM to now */
        double lst;

        /* on the assumption that the user will stick with their chosen display
         * epoch for a while, we move the defining values to match and avoid
         * precession for every call until it is changed again.
         * N.B. only compare and store jd's to lowest precission (f_epoch).
         * N.B. maintaining J2k ref (which is arbitrary) helps avoid accum err
         */
        if (epoch != EOD && (float)epoch != (float)op->f_epoch) {
            double pr = op->f_RA, pd = op->f_dec, fe = (float)epoch;
            /* first bring back to 2k */
            precess (op->f_epoch, J2000, &pr, &pd);
            pr += op->f_pmRA*(J2000-op->f_epoch);
            pd += op->f_pmdec*(J2000-op->f_epoch);
            /* then to epoch */
            pr += op->f_pmRA*(fe-J2000);
            pd += op->f_pmdec*(fe-J2000);
            precess (J2000, fe, &pr, &pd);
            op->f_RA = (float)pr;
            op->f_dec = (float)pd;
            op->f_epoch = (float)fe;
        }

        /* apply proper motion .. assume pm epoch reference equals equinox */
        rpm = op->f_RA + op->f_pmRA*(mjd-op->f_epoch);
        dpm = op->f_dec + op->f_pmdec*(mjd-op->f_epoch);

        /* set ra/dec to astrometric @ equinox of date */
        ra = rpm;
        dec = dpm;
        precess (op->f_epoch, mjed, &ra, &dec);

        /* convert equatoreal ra/dec to mean geocentric ecliptic lat/long */
        eq_ecl (mjed, ra, dec, &bet, &lam);

        /* find solar ecliptical long.(mean equinox) and distance from earth */
        sunpos (mjed, &lsn, &rsn, NULL);

        /* allow for relativistic light bending near the sun */
        deflect (mjed, lam, bet, lsn, rsn, 1e10, &ra, &dec);

        /* TODO: correction for annual parallax would go here */

        /* correct EOD equatoreal for nutation/aberation to form apparent 
         * geocentric
         */
        nut_eq(mjed, &ra, &dec);
        ab_eq(mjed, lsn, &ra, &dec);
        op->s_gaera = (float)ra;
        op->s_gaedec = (float)dec;

        /* set s_ra/dec -- apparent if EOD else astrometric */
        if (epoch == EOD) {
            op->s_ra = (float)ra;
            op->s_dec = (float)dec;
        } else {
            /* annual parallax at time mjd is to be added here, too, but
             * technically in the frame of equinox (usually different from mjd)
             */
            op->s_ra = rpm;
            op->s_dec = dpm;
        }

        /* compute elongation from ecliptic long/lat and sun geocentric long */
        elongation (lam, bet, lsn, &el);
        el = raddeg(el);
        op->s_elong = (float)el;

        /* these are really the same fields ...
        op->s_mag = op->f_mag;
        op->s_size = op->f_size;
        */

        /* alt, az: correct for refraction; use eod ra/dec. */
        now_lst (np, &lst);
        ha = hrrad(lst) - ra;
        hadec_aa (lat, ha, dec, &alt, &az);
        refract (pressure, temp, alt, &alt);
        op->s_alt = alt;
        op->s_az = az;

        return (0);
}

/* compute sky circumstances of an object in heliocentric elliptic orbit at *np.
 */
static int
obj_elliptical (Now *np, Obj *op)
{
        double lsn, rsn;        /* true geoc lng of sun; dist from sn to earth*/
        double dt;              /* light travel time to object */
        double lg;              /* helio long of earth */
        double nu;              /* true anomaly */
        double rp=0;            /* distance from the sun */
        double lo, slo, clo;    /* angle from ascending node */
        double inc;             /* inclination */
        double psi=0;           /* heliocentric latitude */
        double spsi=0, cpsi=0;  /* trig of heliocentric latitude */
        double lpd;             /* heliocentric longitude */
        double rho=0;           /* distance from the Earth */
        double om;              /* arg of perihelion */
        double Om;              /* long of ascending node. */
        double lam;             /* geocentric ecliptic longitude */
        double bet;             /* geocentric ecliptic latitude */
        double ll=0, sll, cll;  /* helio angle between object and earth */
        double mag;             /* magnitude */
        double e_n;             /* mean daily motion */
        double tp;              /* time from perihelion (days) */
        double rpd=0;
        double y;
        int pass;

        /* find location of earth from sun now */
        sunpos (mjed, &lsn, &rsn, 0);
        lg = lsn + PI;

        /* mean daily motion is derived fro mean distance */
        e_n = 0.9856076686/pow((double)op->e_a, 1.5);

        /* correct for light time by computing position at time mjd, then
         *   again at mjd-dt, where
         *   dt = time it takes light to travel earth-object distance.
         */
        dt = 0;
        for (pass = 0; pass < 2; pass++) {

            reduce_elements (op->e_epoch, mjd-dt, degrad(op->e_inc),
                                        degrad (op->e_om), degrad (op->e_Om),
                                        &inc, &om, &Om);

            tp = mjed - dt - (op->e_cepoch - op->e_M/e_n);
            if (vrc (&nu, &rp, tp, op->e_e, op->e_a*(1-op->e_e)) < 0)
                op->o_flags |= NOCIRCUM;
            nu = degrad(nu);
            lo = nu + om;
            slo = sin(lo);
            clo = cos(lo);
            spsi = slo*sin(inc);
            y = slo*cos(inc);
            psi = asin(spsi);
            lpd = atan(y/clo)+Om;
            if (clo<0) lpd += PI;
            range (&lpd, 2*PI);
            cpsi = cos(psi);
            rpd = rp*cpsi;
            ll = lpd-lg;
            rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll));

            dt = rho*LTAU/3600.0/24.0;  /* light travel time, in days / AU */
        }

        /* compute sin and cos of ll */
        sll = sin(ll);
        cll = cos(ll);

        /* find geocentric ecliptic longitude and latitude */
        if (rpd < rsn)
            lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI;
        else
            lam = atan(rsn*sll/(rpd-rsn*cll))+lpd;
        range (&lam, 2*PI);
        bet = atan(rpd*spsi*sin(lam-lpd)/(cpsi*rsn*sll));

        /* fill in all of op->s_* stuff except s_size and s_mag */
        cir_sky (np, lpd, psi, rp, &rho, lam, bet, lsn, rsn, op);

        /* compute magnitude and size */
        if (op->e_mag.whichm == MAG_HG) {
            /* the H and G parameters from the Astro. Almanac.
             */
            if (op->e_size)
                op->s_size = (float)(op->e_size / rho);
            else {
                hg_mag (op->e_mag.m1, op->e_mag.m2, rp, rho, rsn, &mag);
                op->s_size = (float)(h_albsize (op->e_mag.m1)/rho);

            }
        } else {
            /* the g/k model of comets */
            gk_mag (op->e_mag.m1, op->e_mag.m2, rp, rho, &mag);
            op->s_size = (float)(op->e_size / rho);
        }
        set_smag (op, mag);

        return (0);
}

/* compute sky circumstances of an object in heliocentric hyperbolic orbit.
 */
static int
obj_hyperbolic (Now *np, Obj *op)
{
        double lsn, rsn;        /* true geoc lng of sun; dist from sn to earth*/
        double dt;              /* light travel time to object */
        double lg;              /* helio long of earth */
        double nu;              /* true anomaly and eccentric anomaly */
        double rp=0;            /* distance from the sun */
        double lo, slo, clo;    /* angle from ascending node */
        double inc;             /* inclination */
        double psi=0;           /* heliocentric latitude */
        double spsi=0, cpsi=0;  /* trig of heliocentric latitude */
        double lpd;             /* heliocentric longitude */
        double rho=0;           /* distance from the Earth */
        double om;              /* arg of perihelion */
        double Om;              /* long of ascending node. */
        double lam;             /* geocentric ecliptic longitude */
        double bet;             /* geocentric ecliptic latitude */
        double e;               /* fast eccentricity */
        double ll=0, sll, cll;  /* helio angle between object and earth */
        double mag;             /* magnitude */
        double a;               /* mean distance */
        double tp;              /* time from perihelion (days) */
        double rpd=0;
        double y;
        int pass;

        /* find solar ecliptical longitude and distance to sun from earth */
        sunpos (mjed, &lsn, &rsn, 0);

        lg = lsn + PI;
        e = op->h_e;
        a = op->h_qp/(e - 1.0);

        /* correct for light time by computing position at time mjd, then
         *   again at mjd-dt, where
         *   dt = time it takes light to travel earth-object distance.
         */
        dt = 0;
        for (pass = 0; pass < 2; pass++) {

            reduce_elements (op->h_epoch, mjd-dt, degrad(op->h_inc),
                            degrad (op->h_om), degrad (op->h_Om),
                            &inc, &om, &Om);

            tp = mjed - dt - op->h_ep;
            if (vrc (&nu, &rp, tp, op->h_e, op->h_qp) < 0)
                op->o_flags |= NOCIRCUM;
            nu = degrad(nu);
            lo = nu + om;
            slo = sin(lo);
            clo = cos(lo);
            spsi = slo*sin(inc);
            y = slo*cos(inc);
            psi = asin(spsi);
            lpd = atan(y/clo)+Om;
            if (clo<0) lpd += PI;
            range (&lpd, 2*PI);
            cpsi = cos(psi);
            rpd = rp*cpsi;
            ll = lpd-lg;
            rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll));

            dt = rho*5.775518e-3;       /* light travel time, in days */
        }

        /* compute sin and cos of ll */
        sll = sin(ll);
        cll = cos(ll);

        /* find geocentric ecliptic longitude and latitude */
        if (rpd < rsn)
            lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI;
        else
            lam = atan(rsn*sll/(rpd-rsn*cll))+lpd;
        range (&lam, 2*PI);
        bet = atan(rpd*spsi*sin(lam-lpd)/(cpsi*rsn*sll));

        /* fill in all of op->s_* stuff except s_size and s_mag */
        cir_sky (np, lpd, psi, rp, &rho, lam, bet, lsn, rsn, op);

        /* compute magnitude and size */
        gk_mag (op->h_g, op->h_k, rp, rho, &mag);
        set_smag (op, mag);
        op->s_size = (float)(op->h_size / rho);

        return (0);
}

/* compute sky circumstances of an object in heliocentric hyperbolic orbit.
 */
static int
obj_parabolic (Now *np, Obj *op)
{
        double lsn, rsn;        /* true geoc lng of sun; dist from sn to earth*/
        double lam;             /* geocentric ecliptic longitude */
        double bet;             /* geocentric ecliptic latitude */
        double mag;             /* magnitude */
        double inc, om, Om;
        double lpd, psi, rp, rho;
        double dt;
        int pass;

        /* find solar ecliptical longitude and distance to sun from earth */
        sunpos (mjed, &lsn, &rsn, 0);

        /* two passes to correct lam and bet for light travel time. */
        dt = 0.0;
        for (pass = 0; pass < 2; pass++) {
            reduce_elements (op->p_epoch, mjd-dt, degrad(op->p_inc),
                degrad(op->p_om), degrad(op->p_Om), &inc, &om, &Om);
            comet (mjed-dt, op->p_ep, inc, om, op->p_qp, Om,
                                    &lpd, &psi, &rp, &rho, &lam, &bet);
            dt = rho*LTAU/3600.0/24.0;  /* light travel time, in days / AU */
        }

        /* fill in all of op->s_* stuff except s_size and s_mag */
        cir_sky (np, lpd, psi, rp, &rho, lam, bet, lsn, rsn, op);

        /* compute magnitude and size */
        gk_mag (op->p_g, op->p_k, rp, rho, &mag);
        set_smag (op, mag);
        op->s_size = (float)(op->p_size / rho);

        return (0);
}

/* find sun's circumstances now.
 */
static int
sun_cir (Now *np, Obj *op)
{
        double lsn, rsn;        /* true geoc lng of sun; dist from sn to earth*/
        double bsn;             /* true latitude beta of sun */
        double dhlong;

        sunpos (mjed, &lsn, &rsn, &bsn);/* sun's true coordinates; mean ecl. */

        op->s_sdist = 0.0;
        op->s_elong = 0.0;
        op->s_phase = 100.0;
        set_smag (op, -26.8);   /* TODO */
        dhlong = lsn-PI;        /* geo- to helio- centric */
        range (&dhlong, 2*PI);
        op->s_hlong = (float)dhlong;
        op->s_hlat = (float)(-bsn);

        /* fill sun's ra/dec, alt/az in op */
        cir_pos (np, bsn, lsn, &rsn, op);
        op->s_edist = (float)rsn;
        op->s_size = (float)(raddeg(4.65242e-3/rsn)*3600*2);

        return (0);
}

/* find moon's circumstances now.
 */
static int
moon_cir (Now *np, Obj *op)
{
        double lsn, rsn;        /* true geoc lng of sun; dist from sn to earth*/
        double lam;             /* geocentric ecliptic longitude */
        double bet;             /* geocentric ecliptic latitude */
        double edistau;         /* earth-moon dist, in au */
        double el;              /* elongation, rads east */
        double ms;              /* sun's mean anomaly */
        double md;              /* moon's mean anomaly */
        double i;

        moon (mjed, &lam, &bet, &edistau, &ms, &md);    /* mean ecliptic & EOD*/
        sunpos (mjed, &lsn, &rsn, NULL);                /* mean ecliptic & EOD*/

        op->s_hlong = (float)lam;               /* save geo in helio fields */
        op->s_hlat = (float)bet;

        /* find angular separation from sun */
        elongation (lam, bet, lsn, &el);
        op->s_elong = (float)raddeg(el);                /* want degrees */

        /* solve triangle of earth, sun, and elongation for moon-sun dist */
        op->s_sdist = (float) sqrt (edistau*edistau + rsn*rsn
                                                    - 2.0*edistau*rsn*cos(el));

        /* TODO: improve mag; this is based on a flat moon model. */
        i = -12.7 + 2.5*(log10(PI) - log10(PI/2*(1+1.e-6-cos(el)))) 
                                        + 5*log10(edistau/.0025) /* dist */;
        set_smag (op, i);

        /* find phase -- allow for projection effects */
        i = 0.1468*sin(el)*(1 - 0.0549*sin(md))/(1 - 0.0167*sin(ms));
        op->s_phase = (float)((1+cos(PI-el-degrad(i)))/2*100);

        /* fill moon's ra/dec, alt/az in op and update for topo dist */
        cir_pos (np, bet, lam, &edistau, op);

        op->s_edist = (float)edistau;
        op->s_size = (float)(3600*2.0*raddeg(asin(MRAD/MAU/edistau)));
                                                /* moon angular dia, seconds */

        return (0);
}

/* fill in all of op->s_* stuff except s_size and s_mag.
 * this is used for sol system objects (except sun and moon); never FIXED.
 */
static void
cir_sky (
Now *np,
double lpd,             /* heliocentric ecliptic longitude */
double psi,             /* heliocentric ecliptic lat */
double rp,              /* dist from sun */
double *rho,            /* dist from earth: in as geo, back as geo or topo */
double lam,             /* true geocentric ecliptic long */
double bet,             /* true geocentric ecliptic lat */
double lsn,             /* true geoc lng of sun */
double rsn,             /* dist from sn to earth*/
Obj *op)
{
        double el;              /* elongation */
        double f;               /* fractional phase from earth */

        /* compute elongation and phase */
        elongation (lam, bet, lsn, &el);
        el = raddeg(el);
        op->s_elong = (float)el;
        f = 0.25 * ((rp+ *rho)*(rp+ *rho) - rsn*rsn)/(rp* *rho);
        op->s_phase = (float)(f*100.0); /* percent */

        /* set heliocentric long/lat; mean ecliptic and EOD */
        op->s_hlong = (float)lpd;
        op->s_hlat = (float)psi;

        /* fill solar sys body's ra/dec, alt/az in op */
        cir_pos (np, bet, lam, rho, op);        /* updates rho */

        /* set earth/planet and sun/planet distance */
        op->s_edist = (float)(*rho);
        op->s_sdist = (float)rp;
}

/* fill equatoreal and horizontal op-> fields; stern
 *
 *    input:          lam/bet/rho geocentric mean ecliptic and equinox of day
 * 
 * algorithm at EOD:
 *   ecl_eq     --> ra/dec      geocentric mean equatoreal EOD (via mean obliq)
 *   deflect    --> ra/dec        relativistic deflection
 *   nut_eq     --> ra/dec      geocentric true equatoreal EOD
 *   ab_eq      --> ra/dec      geocentric apparent equatoreal EOD
 *                                      if (PREF_GEO)  --> output
 *   ta_par     --> ra/dec      topocentric apparent equatoreal EOD
 *                                      if (!PREF_GEO)  --> output
 *   hadec_aa   --> alt/az      topocentric horizontal
 *   refract    --> alt/az      observed --> output
 *
 * algorithm at fixed equinox:
 *   ecl_eq     --> ra/dec      geocentric mean equatoreal EOD (via mean obliq)
 *   deflect    --> ra/dec        relativistic deflection [for alt/az only]
 *   nut_eq     --> ra/dec      geocentric true equatoreal EOD [for aa only]
 *   ab_eq      --> ra/dec      geocentric apparent equatoreal EOD [for aa only]
 *   ta_par     --> ra/dec      topocentric apparent equatoreal EOD
 *     precess  --> ra/dec      topocentric equatoreal fixed equinox [eq only]
 *                                      --> output
 *   hadec_aa   --> alt/az      topocentric horizontal
 *   refract    --> alt/az      observed --> output
 */
static void
cir_pos (
Now *np,
double bet,     /* geo lat (mean ecliptic of date) */
double lam,     /* geo long (mean ecliptic of date) */
double *rho,    /* in: geocentric dist in AU; out: geo- or topocentic dist */
Obj *op)        /* object to set s_ra/dec as per equinox */
{
        double ra, dec;         /* apparent ra/dec, corrected for nut/ab */
        double tra, tdec;       /* astrometric ra/dec, no nut/ab */
        double lsn, rsn;        /* solar geocentric (mean ecliptic of date) */
        double ha_in, ha_out;   /* local hour angle before/after parallax */
        double dec_out;         /* declination after parallax */
        double dra, ddec;       /* parallax correction */
        double alt, az;         /* current alt, az */
        double lst;             /* local sidereal time */
        double rho_topo;        /* topocentric distance in earth radii */

        /* convert to equatoreal [mean equator, with mean obliquity] */
        ecl_eq (mjed, bet, lam, &ra, &dec);
        tra = ra;       /* keep mean coordinates */
        tdec = dec;

        /* get sun position */
        sunpos(mjed, &lsn, &rsn, NULL);

        /* allow for relativistic light bending near the sun.
         * (avoid calling deflect() for the sun itself).
         */
        if (!is_planet(op,SUN) && !is_planet(op,MOON))
            deflect (mjed, op->s_hlong, op->s_hlat, lsn, rsn, *rho, &ra, &dec);

        /* correct ra/dec to form geocentric apparent */
        nut_eq (mjed, &ra, &dec);
        if (!is_planet(op,MOON))
            ab_eq (mjed, lsn, &ra, &dec);
        op->s_gaera = (float)ra;
        op->s_gaedec = (float)dec;

        /* find parallax correction for equatoreal coords */
        now_lst (np, &lst);
        ha_in = hrrad(lst) - ra;
        rho_topo = *rho * MAU/ERAD;             /* convert to earth radii */
        ta_par (ha_in, dec, lat, elev, &rho_topo, &ha_out, &dec_out);

        /* transform into alt/az and apply refraction */
        hadec_aa (lat, ha_out, dec_out, &alt, &az);
        refract (pressure, temp, alt, &alt);
        op->s_alt = alt;
        op->s_az = az;

        /* Get parallax differences and apply to apparent or astrometric place
         * as needed.  For the astrometric place, rotating the CORRECTIONS
         * back from the nutated equator to the mean equator will be
         * neglected.  This is an effect of about 0.1" at moon distance.
         * We currently don't have an inverse nutation rotation.
         */
        if (pref_get(PREF_EQUATORIAL) == PREF_GEO) {
            /* no topo corrections to eq. coords */
            dra = ddec = 0.0;
        } else {
            dra = ha_in - ha_out;       /* ra sign is opposite of ha */
            ddec = dec_out - dec;
            *rho = rho_topo * ERAD/MAU; /* return topocentric distance in AU */
        }

        /* fill in ra/dec fields */
        if (epoch == EOD) {             /* apparent geo/topocentric */
            ra = ra + dra;
            dec = dec + ddec;
        } else {                        /* astrometric geo/topocent */
            ra = tra + dra;
            dec = tdec + ddec;
            precess (mjed, epoch, &ra, &dec);
        }
        range(&ra, 2*PI);
        op->s_ra = (float)ra;
        op->s_dec = (float)dec;
}

/* given geocentric ecliptic longitude and latitude, lam and bet, of some object
 * and the longitude of the sun, lsn, find the elongation, el. this is the
 * actual angular separation of the object from the sun, not just the difference
 * in the longitude. the sign, however, IS set simply as a test on longitude
 * such that el will be >0 for an evening object <0 for a morning object.
 * to understand the test for el sign, draw a graph with lam going from 0-2*PI
 *   down the vertical axis, lsn going from 0-2*PI across the hor axis. then
 *   define the diagonal regions bounded by the lines lam=lsn+PI, lam=lsn and
 *   lam=lsn-PI. the "morning" regions are any values to the lower left of the
 *   first line and bounded within the second pair of lines.
 * all angles in radians.
 */
static void
elongation (double lam, double bet, double lsn, double *el)
{
        *el = acos(cos(bet)*cos(lam-lsn));
        if (lam>lsn+PI || (lam>lsn-PI && lam<lsn)) *el = - *el;
}

/* apply relativistic light bending correction to ra/dec; stern
 *
 * The algorithm is from:
 * Mean and apparent place computations in the new IAU 
 * system. III - Apparent, topocentric, and astrometric 
 * places of planets and stars
 * KAPLAN, G. H.;  HUGHES, J. A.;  SEIDELMANN, P. K.;
 * SMITH, C. A.;  YALLOP, B. D.
 * Astronomical Journal (ISSN 0004-6256), vol. 97, April 1989, p. 1197-1210.
 *
 * This article is a very good collection of formulea for geocentric and
 * topocentric place calculation in general.  The apparent and
 * astrometric place calculation in this file currently does not follow
 * the strict algorithm from this paper and hence is not fully correct.
 * The entire calculation is currently based on the rotating EOD frame and
 * not the "inertial" J2000 frame.
 */
static void
deflect (
double mjd1,            /* equinox */
double lpd, double psi, /* heliocentric ecliptical long / lat */
double rsn, double lsn, /* distance and longitude of sun */
double rho,             /* geocentric distance */
double *ra, double *dec)/* geocentric equatoreal */
{
        double hra, hdec;       /* object heliocentric equatoreal */
        double el;              /* HELIOCENTRIC elongation object--earth */
        double g1, g2;          /* relativistic weights */
        double u[3];            /* object geocentric cartesian */
        double q[3];            /* object heliocentric cartesian unit vect */
        double e[3];            /* earth heliocentric cartesian unit vect */
        double qe, uq, eu;      /* scalar products */
        int i;                  /* counter */

#define G       1.32712438e20   /* heliocentric grav const; in m^3*s^-2 */
#define c       299792458.0     /* speed of light in m/s */

        elongation(lpd, psi, lsn-PI, &el);
        el = fabs(el);
        /* only continue if object is within about 10 deg around the sun,
         * not obscured by the sun's disc (radius 0.25 deg) and farther away
         * than the sun.
         *
         * precise geocentric deflection is:  g1 * tan(el/2)
         *      radially outwards from sun;  the vector munching below
         *      just applys this component-wise
         *      Note:   el = HELIOCENTRIC elongation.
         *              g1 is always about 0.004 arc seconds
         *              g2 varies from 0 (highest contribution) to 2
         */
        if (el<degrad(170) || el>degrad(179.75) || rho<rsn) return;

        /* get cartesian vectors */
        sphcart(*ra, *dec, rho, u, u+1, u+2);

        ecl_eq(mjd1, psi, lpd, &hra, &hdec);
        sphcart(hra, hdec, 1.0, q, q+1, q+2);

        ecl_eq(mjd1, 0.0, lsn-PI, &hra, &hdec);
        sphcart(hra, hdec, 1.0, e, e+1, e+2);

        /* evaluate scalar products */
        qe = uq = eu = 0.0;
        for(i=0; i<=2; ++i) {
            qe += q[i]*e[i];
            uq += u[i]*q[i];
            eu += e[i]*u[i];
        }

        g1 = 2*G/(c*c*MAU)/rsn;
        g2 = 1 + qe;

        /* now deflect geocentric vector */
        g1 /= g2;
        for(i=0; i<=2; ++i)
            u[i] += g1*(uq*e[i] - eu*q[i]);
        
        /* back to spherical */
        cartsph(u[0], u[1], u[2], ra, dec, &rho);       /* rho thrown away */
}

/* estimate size in arc seconds @ 1AU from absolute magnitude, H, and assuming
 * an albedo of 0.1. With this assumption an object with diameter of 1500m
 * has an absolute mag of 18.
 */
static double
h_albsize (double H)
{
        return (3600*raddeg(.707*1500*pow(2.51,(18-H)/2)/MAU));
}

/* For RCS Only -- Do Not Edit */
static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: circum.c,v $ $Date: 2009-07-16 10:34:36 $ $Revision: 1.9 $ $Name: not supported by cvs2svn $"};