Index: trunk/SophyaExt/XephemAstroLib/circum.c
===================================================================
--- trunk/SophyaExt/XephemAstroLib/circum.c	(revision 1719)
+++ trunk/SophyaExt/XephemAstroLib/circum.c	(revision 2551)
@@ -9,28 +9,29 @@
 #include <stdio.h>
 #include <math.h>
-#if defined(__STDC__)
 #include <stdlib.h>
-#endif
-
-#include "P_.h"
+
 #include "astro.h"
-#include "circum.h"
 #include "preferences.h"
 
 
-static int obj_planet P_((Now *np, Obj *op));
-static int obj_fixed P_((Now *np, Obj *op));
-static int obj_elliptical P_((Now *np, Obj *op));
-static int obj_hyperbolic P_((Now *np, Obj *op));
-static int obj_parabolic P_((Now *np, Obj *op));
-static int sun_cir P_((Now *np, Obj *op));
-static int moon_cir P_((Now *np, Obj *op));
-static void cir_sky P_((Now *np, double lpd, double psi, double rp, double *rho,
-    double lam, double bet, double lsn, double rsn, Obj *op));
-static void cir_pos P_((Now *np, double bet, double lam, double *rho, Obj *op));
-static void elongation P_((double lam, double bet, double lsn, double *el));
-static void deflect P_((double mjd1, double lpd, double psi, double rsn,
-    double lsn, double rho, double *ra, double *dec));
-static double h_albsize P_((double 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.
@@ -38,9 +39,9 @@
  */
 int
-obj_cir (np, op)
-Now *np;
-Obj *op;
-{
+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));
@@ -50,6 +51,6 @@
 	case PLANET:     return (obj_planet (np, op));
 	default:
-	    printf ("obj_cir() called with type %d\n", op->o_type);
-	    exit(1);
+	    printf ("obj_cir() called with type %d %s\n", op->o_type, op->o_name);
+	    abort();
 	    return (-1);	/* just for lint */
 	}
@@ -57,7 +58,5 @@
 
 static int
-obj_planet (np, op)
-Now *np;
-Obj *op;
+obj_planet (Now *np, Obj *op)
 {
 	double lsn, rsn;	/* true geoc lng of sun; dist from sn to earth*/
@@ -67,16 +66,20 @@
 	double lam, bet;	/* geocentric ecliptic long and lat */
 	double dia, mag;	/* angular diameter at 1 AU and magnitude */
-	int p;
+	PLCode p;
 
 	/* validate code and check for a few special cases */
-	p = op->pl.pl_code;
+	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);
-	    exit(1);
-	}
-	else if (p == SUN)
-	    return (sun_cir (np, op));
-	else if (p == MOON)
-	    return (moon_cir (np, op));
+	    abort();
+	}
+
+	/* planet itself */
 
 	/* find solar ecliptical longitude and distance to sun from earth */
@@ -99,7 +102,109 @@
 
 static int
-obj_fixed (np, op)
-Now *np;
-Obj *op;
+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*/
@@ -108,23 +213,36 @@
 	double el;		/* elongation */
 	double alt, az;		/* current alt, az */
-	double ra, dec;		/* ra and dec at epoch of date */
+	double ra, dec;		/* ra and dec at equinox of date */
+	double rpm, dpm; 	/* astrometric ra and dec with PM to now */
 	double lst;
 
-	if (epoch != EOD && (float)epoch != op->f_epoch) {
-	    /* want a certain epoch -- if it's not what the database is at
-	     * we change the original to save time next time assuming the
-	     * user is likely to stick with this for a while.
-	     */
-	    double tra = op->f_RA, tdec = op->f_dec;
-	    float tepoch = (float)epoch;	/* compare w/float precision */
-	    precess (op->f_epoch, tepoch, &tra, &tdec);
-	    op->f_epoch = tepoch;
-	    op->f_RA = (float)tra;
-	    op->f_dec = (float)tdec;
-	}
-
-	/* set ra/dec to astrometric @ epoch of date */
-	ra = op->f_RA;
-	dec = op->f_dec;
+	/* 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);
 
@@ -154,8 +272,8 @@
 	} else {
 	    /* annual parallax at time mjd is to be added here, too, but
-	     * technically in the frame of epoch (usually different from mjd)
+	     * technically in the frame of equinox (usually different from mjd)
 	     */
-	    op->s_ra = op->f_RA;	/* already precessed */
-	    op->s_dec = op->f_dec;
+	    op->s_ra = rpm;
+	    op->s_dec = dpm;
 	}
 
@@ -184,7 +302,5 @@
  */
 static int
-obj_elliptical (np, op)
-Now *np;
-Obj *op;
+obj_elliptical (Now *np, Obj *op)
 {
 	double lsn, rsn;	/* true geoc lng of sun; dist from sn to earth*/
@@ -230,5 +346,6 @@
 
 	    tp = mjed - dt - (op->e_cepoch - op->e_M/e_n);
-	    vrc (&nu, &rp, tp, op->e_e, op->e_a*(1-op->e_e));
+	    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;
@@ -288,7 +405,5 @@
  */
 static int
-obj_hyperbolic (np, op)
-Now *np;
-Obj *op;
+obj_hyperbolic (Now *np, Obj *op)
 {
 	double lsn, rsn;	/* true geoc lng of sun; dist from sn to earth*/
@@ -309,5 +424,4 @@
 	double e;		/* fast eccentricity */
 	double ll=0, sll, cll;	/* helio angle between object and earth */
-	double n;		/* mean daily motion */
 	double mag;		/* magnitude */
 	double a;		/* mean distance */
@@ -323,5 +437,4 @@
 	e = op->h_e;
 	a = op->h_qp/(e - 1.0);
-	n = .98563/sqrt(a*a*a);
 
 	/* correct for light time by computing position at time mjd, then
@@ -337,5 +450,6 @@
 
 	    tp = mjed - dt - op->h_ep;
-	    vrc (&nu, &rp, tp, op->h_e, op->h_qp);
+	    if (vrc (&nu, &rp, tp, op->h_e, op->h_qp) < 0)
+		op->o_flags |= NOCIRCUM;
 	    nu = degrad(nu);
 	    lo = nu + om;
@@ -382,7 +496,5 @@
  */
 static int
-obj_parabolic (np, op)
-Now *np;
-Obj *op;
+obj_parabolic (Now *np, Obj *op)
 {
 	double lsn, rsn;	/* true geoc lng of sun; dist from sn to earth*/
@@ -422,7 +534,5 @@
  */
 static int
-sun_cir (np, op)
-Now *np;
-Obj *op;
+sun_cir (Now *np, Obj *op)
 {
 	double lsn, rsn;	/* true geoc lng of sun; dist from sn to earth*/
@@ -452,7 +562,5 @@
  */
 static int
-moon_cir (np, op)
-Now *np;
-Obj *op;
+moon_cir (Now *np, Obj *op)
 {
 	double lsn, rsn;	/* true geoc lng of sun; dist from sn to earth*/
@@ -480,5 +588,7 @@
 
 	/* TODO: improve mag; this is based on a flat moon model. */
-	set_smag (op, -12.7 + 2.5*(log10(PI) - log10(PI/2*(1+1.e-6-cos(el)))));
+	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 */
@@ -500,12 +610,15 @@
  */
 static void
-cir_sky (np, lpd, psi, rp, rho, lam, bet, lsn, rsn, op)
-Now *np;
-double lpd, psi;	/* heliocentric ecliptic long and lat */
-double rp;		/* dist from sun */
-double *rho;		/* dist from earth: in as geo, back as geo or topo */
-double lam, bet;	/* true geocentric ecliptic long and lat */
-double lsn, rsn;	/* true geoc lng of sun; dist from sn to earth*/
-Obj *op;
+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 */
@@ -546,5 +659,5 @@
  *   refract	--> alt/az	observed --> output
  *
- * algorithm at fixed epoch:
+ * algorithm at fixed equinox:
  *   ecl_eq	--> ra/dec	geocentric mean equatoreal EOD (via mean obliq)
  *   deflect	--> ra/dec	  relativistic deflection [for alt/az only]
@@ -558,9 +671,10 @@
  */
 static void
-cir_pos (np, bet, lam, rho, op)
-Now *np;
-double bet, lam;/* geo lat/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 epoch */
+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 */
@@ -649,7 +763,5 @@
  */
 static void
-elongation (lam, bet, lsn, el)
-double lam, bet, lsn;
-double *el;
+elongation (double lam, double bet, double lsn, double *el)
 {
 	*el = acos(cos(bet)*cos(lam-lsn));
@@ -675,10 +787,10 @@
  */
 static void
-deflect (mjd1, lpd, psi, lsn, rsn, rho, ra, dec)
-double mjd1;		/* epoch */
-double lpd, psi;	/* heliocentric ecliptical long / lat */
-double rsn, lsn;	/* distance and longitude of sun */
-double rho;		/* geocentric distance */
-double *ra, *dec;	/* geocentric equatoreal */
+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 */
@@ -742,6 +854,5 @@
  */
 static double
-h_albsize (H)
-double H;
+h_albsize (double H)
 {
 	return (3600*raddeg(.707*1500*pow(2.51,(18-H)/2)/MAU));
@@ -749,3 +860,3 @@
 
 /* For RCS Only -- Do Not Edit */
-static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: circum.c,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
+static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: circum.c,v $ $Date: 2004-06-15 16:52:37 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};
