#include <math.h>
#include <stdio.h>
#include <errno.h>
#include <stdarg.h>
#include <stdlib.h>
#include <strings.h>

#include "ephlib.h"

extern const char *__progname;

#define SIGN_WIDTH 6 /* degrees */

#define RADTODEG (180/M_PI)
#define DEGTORAD (M_PI/180)

#define JOFFSET_DAY 1721060
#define JOFFSET_SEC (12*60*60)
#define SECPERDAY (24*60*60)

#define EARTH_RADIUS (6371.01/149597870.691) /* rad in km / au in km */

#define ZERO_ARIES 0.040078256357732808257 /* radians */

typedef unsigned int PLANETSET;

/* Grand aspect parameters. */
#define GA_MAX_ASPS    6 /* maximum aspects involved */
#define GA_MAX_GAS     2 /* max # of grand aspects per array element */
#define GA_MAX_PLANETS 4 /* max # of planets involved */

typedef struct planet PLANET;
typedef struct gvs GVS;
typedef struct xyz XYZ;
typedef struct date DATE;
typedef struct time TIME;
typedef struct sign SIGN;
typedef struct aspect ASPECT;
typedef struct asplist ASPLIST;
typedef struct grandaspect GRANDASPECT;
typedef struct gaasp GAASP;
typedef struct gaga GAGA;
typedef struct gafound GAFOUND;

struct date {
  int y;
  int m;
  int d;
  int j;
  } ;

struct time {
  int h;
  int m;
  int s;
  } ;

struct xyz {
  double x;
  double y;
  double z;
  } ;

struct gvs {
  XYZ g;
  XYZ e;
  double r;
  double theta;
  double phi;
  } ;

struct planet {
  const char *name;
  int ephtag;
  GVS loc;
  GVS vel;
  } ;

struct sign {
  const char *name;
  int start;
  } ;

struct asplist {
  int i;
  int j;
  double d;
  } ;

struct aspect {
  const char *name;
  int ang;
  double strongorb;
  double weakorb;
  double noorb;
  int listn;
  ASPLIST *list;
  } ;

struct gaasp {
  int p1;
  int p2;
  double ang;
  double orb;
  } ;

struct gaga {
  int nasps;
  const char *name;
  } ;

struct gafound {
  GAGA *ga;
  PLANETSET planets;
  } ;

struct grandaspect {
  int nplanets;
  GAASP asps[GA_MAX_ASPS];
  int ngas;
  GAGA gas[GA_MAX_GAS];
  int nfound;
  GAFOUND *found;
  } ;

static void *h;
static double jtime;
static double stime;
static DATE ce_date;
static TIME ce_time;
static int ce_tz;
static int textmode = 0;
static int allaspects = 0;
static int colourful = 0;

static double ecliptic_axis[2];
static double ecliptic_ortho[2];
static double ecliptic_sin;
static double ecliptic_cos;

static double latitude;
static double longitude;
static int locgot;
#define LG_LAT  1
#define LG_LONG 2
static XYZ surface_loc;
static XYZ surface_vel;

static PLANET planets[]
 = { { "Sun", EPH_SUN },
#define PLANET_SUN 0
     { "Moon", EPH_MOON },
     { "Mercury", EPH_MERCURY },
     { "Venus", EPH_VENUS },
     { "Mars", EPH_MARS },
     { "Jupiter", EPH_JUPITER },
     { "Saturn", EPH_SATURN },
     { "Uranus", EPH_URANUS },
     { "Neptune", EPH_NEPTUNE },
     { "Pluto", EPH_PLUTO } };
#define NPLANETS (sizeof(planets)/sizeof(planets[0]))
#define PLANETPAIRS (NPLANETS*(NPLANETS-1)/2)

static SIGN signs[]
 = { { "Aries", 0 },
     { "Taurus", 30 },
     { "Gemini", 60 },
     { "Cancer", 90 },
     { "Leo", 120 },
     { "Virgo", 150 },
     { "Libra", 180 },
     { "Scorpio", 210 },
     { "Sagittarius", 240 },
     { "Capricorn", 270 },
     { "Aquarius", 300 },
     { "Pisces", 330 } };
#define NSIGNS (sizeof(signs)/sizeof(signs[0]))

static ASPECT aspects[]
  = { { "conjunction",   0, 4, 5, 10 },
      { "opposition",  180, 5, 7, 15 },
      { "trine",       120, 5, 7, 15 },
      { "square",       90, 5, 7, 15 },
      { "sextile",      60, 4, 5, 10 },
      { "quincunx",    150, 2, 3,  6 } };
#define NASPECTS (sizeof(aspects)/sizeof(aspects[0]))

static GRANDASPECT grandaspects[]
 = { { 4,
       { { 0, 1, 120, 7 },
	 { 1, 2, 120, 7 },
	 { 2, 0, 120, 7 },
	 { 3, 0, 180, 7 },
	 { 3, 1, 60, 5 },
	 { 3, 2, 60, 5 } },
       2,
       { { 3, "Grand Trine" },
	 { 6, "Kite" } } },
     { 4,
       { { 0, 1, 90, 7 },
	 { 1, 2, 90, 7 },
	 { 0, 2, 180, 7 },
	 { 0, 3, 90, 7 },
	 { 1, 3, 180, 7 },
	 { 3, 2, 90, 7 } },
       2,
       { { 3, "T-Square" },
	 { 6, "Grand Cross" } } },
     { 3,
       { { 0, 1, 150, 3 },
	 { 0, 2, 150, 3 },
	 { 1, 2, 60, 5 } },
       1,
       { { 3, "Yod" } } } };
#define NGRANDASPECTS (sizeof(grandaspects)/sizeof(grandaspects[0]))

static double abs_angdiff(double a, double b)
{
 a = fabs(a-b);
 while (a > 180) a -= 360;
 return(fabs(a));
}

static void invert_3x3(double (*in)[3], double (*out)[3])
{
 int i;
 int j;
 int k;
 double s;
 double m[6][3];

 for (i=0;i<3;i++) for (j=0;j<3;j++)
  { m[i][j] = in[i][j];
    m[i+3][j] = (i == j);
  }
 for (i=0;i<3;i++)
  { k = i;
    if ((i < 1) && (fabs(m[i][1]) > fabs(m[i][0]))) k = 1;
    if ((i < 2) && (fabs(m[i][2]) > fabs(m[i][k]))) k = 2;
    if (k != i)
     { for (j=0;j<6;j++)
	{ s = m[j][i];
	  m[j][i] = m[j][k];
	  m[j][k] = s;
	}
     }
    s = m[i][i];
    for (j=0;j<6;j++) m[j][i] /= s;
    for (j=0;j<3;j++)
     { if (j == i) continue;
       s = m[i][j];
       for (k=i;k<6;k++) m[k][j] -= s * m[k][i];
     }
  }
 for (i=0;i<3;i++) for (j=0;j<3;j++) out[i][j] = m[i+3][j];
}

static void do_ecliptic(GVS *gvs)
{
 double f;
 double o;
 double no;

 f = (gvs->g.x * ecliptic_axis[0]) + (gvs->g.y * ecliptic_axis[1]);
 o = (gvs->g.x * ecliptic_ortho[0]) + (gvs->g.y * ecliptic_ortho[1]);
 no = (o * ecliptic_cos) - (gvs->g.z * ecliptic_sin);
 gvs->e.x = (f * ecliptic_axis[0]) + (no * ecliptic_ortho[0]);
 gvs->e.y = (f * ecliptic_axis[1]) + (no * ecliptic_ortho[1]);
 gvs->e.z = (o * ecliptic_sin) + (gvs->g.z * ecliptic_cos);
}

static void do_spherical(PLANET *p)
{
 double x;
 double y;
 double z;
 double vx;
 double vy;
 double vz;
 double xyr;
 double v[3][3];
 double vi[3][3];
 double r;
 double cos_theta;
 double sin_theta;
 double cos_phi;
 double sin_phi;

 x = p->loc.e.x - surface_loc.x;
 y = p->loc.e.y - surface_loc.y;
 z = p->loc.e.z - surface_loc.z;
 vx = p->vel.e.x - surface_vel.x;
 vy = p->vel.e.y - surface_vel.y;
 vz = p->vel.e.z - surface_vel.z;
 xyr = hypot(x,y);
 r = hypot(xyr,z);
 p->loc.r = r;
 p->loc.theta = atan2(y,x) * RADTODEG;
 p->loc.phi = atan2(z,xyr) * RADTODEG;
 cos_theta = x / xyr;
 sin_theta = y / xyr;
 cos_phi = xyr / r;
 sin_phi = z / r;
 v[0][0] =       cos_theta * cos_phi;
 v[1][0] = - r * sin_theta * cos_phi;
 v[2][0] = - r * cos_theta * sin_phi;
 v[0][1] =       sin_theta * cos_phi;
 v[1][1] =   r * cos_theta * cos_phi;
 v[2][1] = - r * sin_theta * sin_phi;
 v[0][2] =                   sin_phi;
 v[1][2] = 0;
 v[2][2] =   r *             cos_phi;
 invert_3x3(&v[0],&vi[0]);
 p->vel.r     =   (vi[0][0] * vx) +
		  (vi[1][0] * vy) +
		  (vi[2][0] * vz);
 p->vel.theta = ( (vi[0][1] * vx) +
		  (vi[1][1] * vy) +
		  (vi[2][1] * vz) ) * RADTODEG;
 p->vel.phi   = ( (vi[0][2] * vx) +
		  (vi[1][2] * vy) +
		  (vi[2][2] * vz) ) * RADTODEG;
 if (! textmode)
  { printf("%% %-7s  spherical  %15.9f %15.9f %15.9f\n",p->name,p->loc.r,p->loc.phi,p->loc.theta);
    printf("%%          sph. vel   %15.9f %15.9f %15.9f\n",p->vel.r,p->vel.phi,p->vel.theta);
  }
}

static void getloc(int x)
{
 PLANET *p;
 double ephr[6];

 p = &planets[x];
 switch (eph_lookup(h,jtime,p->ephtag,EPH_EARTH,&ephr[0]))
  { case 0:
       break;
    case EPHL_ERR_RANGE:
       fprintf(stderr,"%s: time %.5f out of range for ephemeris\n",__progname,jtime);
       exit(1);
       break;
    default:
       abort();
       break;
  }
 p->loc.g.x = ephr[0];
 p->loc.g.y = ephr[1];
 p->loc.g.z = ephr[2];
 p->vel.g.x = ephr[3];
 p->vel.g.y = ephr[4];
 p->vel.g.z = ephr[5];
 if (! textmode)
  { printf("%% %-7s  geocentric %15.9f %15.9f %15.9f\n",p->name,p->loc.g.x,p->loc.g.y,p->loc.g.z);
    printf("%%          geo. vel   %15.9f %15.9f %15.9f\n",p->vel.g.x,p->vel.g.y,p->vel.g.z);
  }
 do_ecliptic(&p->loc);
 do_ecliptic(&p->vel);
 if (! textmode)
  { printf("%%          ecliptic   %15.9f %15.9f %15.9f\n",p->loc.e.x,p->loc.e.y,p->loc.e.z);
    printf("%%          ecl. vel   %15.9f %15.9f %15.9f\n",p->vel.e.x,p->vel.e.y,p->vel.e.z);
  }
}

static void compute_ecliptic(void)
{
 int i;
 int j;
 int k;
 double s;
 double r;
 double a[3];
 double b[3];
 double v[6];
 double XpX[3][3];
 double XpXinv[3][3];
 double XpY[3];
 double B[3];
 double err2;
 double N[3];
 double X[365][3];
 double Y[365];

 if (! textmode) printf("%% Computing ecliptic fixup\n");
 for (i=0;i<365;i++)
  { eph_lookup(h,jtime-180+i,EPH_SUN,EPH_EARTH,&v[0]);
    X[i][0] = v[0];
    X[i][1] = v[1];
    X[i][2] = 1;
    Y[i] = v[2];
  }
 for (i=0;i<3;i++)
  { for (j=0;j<3;j++)
     { s = 0;
       for (k=0;k<365;k++) s += X[k][i] * X[k][j];
       XpX[i][j] = s;
     }
  }
 if (! textmode)
  { printf("%%           / %20.14f %20.14f %20.14f \\\n",XpX[0][0],XpX[1][0],XpX[2][0]);
    printf("%% XpX    = <  %20.14f %20.14f %20.14f  >\n",XpX[0][1],XpX[1][1],XpX[2][1]);
    printf("%%           \\ %20.14f %20.14f %20.14f /\n",XpX[0][2],XpX[1][2],XpX[2][2]);
  }
 invert_3x3(&XpX[0],&XpXinv[0]);
 if (! textmode)
  { printf("%%           / %20.16f %20.16f %20.16f \\\n",XpXinv[0][0],XpXinv[1][0],XpXinv[2][0]);
    printf("%% XpXinv = <  %20.16f %20.16f %20.16f  >\n",XpXinv[0][1],XpXinv[1][1],XpXinv[2][1]);
    printf("%%           \\ %20.16f %20.16f %20.16f /\n",XpXinv[0][2],XpXinv[1][2],XpXinv[2][2]);
  }
 for (i=0;i<3;i++)
  { s = 0;
    for (j=0;j<365;j++)
     { s += X[j][i] * Y[j];
     }
    XpY[i] = s;
  }
 if (! textmode) printf("%% XpY    = [ %20.16f %20.16f %20.16f ]\n",XpY[0],XpY[1],XpY[2]);
 for (i=0;i<3;i++)
  { s = 0;
    for (j=0;j<3;j++)
     { s += XpXinv[j][i] * XpY[j];
     }
    B[i] = s;
  }
 if (! textmode) printf("%% B      = [ %20.17f %20.17f %20.17f ]\n",B[0],B[1],B[2]);
 s = 0;
 for (i=0;i<365;i++)
  { r = (X[i][0] * B[0]) + (X[i][1] * B[1]) + B[2] - Y[i];
    s += r * r;
  }
 err2 = s / 365;
 if (! textmode) printf("%% err2   = %20.18f\n",err2);
 a[0] = X[140][0];
 a[1] = X[140][1];
 a[2] = (a[0] * B[0]) + (a[1] * B[1]) + B[2];
 b[0] = X[230][0];
 b[1] = X[230][1];
 b[2] = (b[0] * B[0]) + (b[1] * B[1]) + B[2];
 N[0] = (a[1] * b[2]) - (a[2] * b[1]);
 N[1] = (a[2] * b[0]) - (a[0] * b[2]);
 N[2] = (a[0] * b[1]) - (a[1] * b[0]);
 s = hypot(hypot(N[0],N[1]),N[2]);
 N[0] /= s;
 N[1] /= s;
 N[2] /= s;
 if (! textmode) printf("%% N      = [ %20.17f %20.17f %20.17f ]\n",N[0],N[1],N[2]);
 s = hypot(N[0],N[1]);
 ecliptic_ortho[0] = N[0] / s;
 ecliptic_ortho[1] = N[1] / s;
 ecliptic_axis[0] = ecliptic_ortho[1];
 ecliptic_axis[1] = -ecliptic_ortho[0];
 ecliptic_sin = sqrt(1-(N[2]*N[2]));
 ecliptic_cos = N[2];
 if (! textmode)
  { printf("%% fixup axis  [ %20.17f %20.17f ]\n",ecliptic_axis[0],ecliptic_axis[1]);
    printf("%%       ortho [ %20.17f %20.17f ]\n",ecliptic_ortho[0],ecliptic_ortho[1]);
    printf("%%       sincos  %20.17f %20.17f\n",ecliptic_sin,ecliptic_cos);
    printf("%% Ecliptic fixup done\n");
  }
}

inline static unsigned int ydays(unsigned int) __attribute__ ((const));
inline static unsigned int ydays(unsigned int y)
{
 return( (y % 4)
	 ? 365
	 : (y % 100)
	 ? 366
	 : (y % 400)
	 ? 365
	 : 366 );
}

inline static unsigned int mdays(unsigned int, unsigned int)
						__attribute__ ((const));
inline static unsigned int mdays(unsigned int y, unsigned int m)
{
 return( (m == 2)
	 ? ( (y % 4)
	     ? 28
	     : (y % 100)
	     ? 29
	     : (y % 400)
	     ? 28
	     : 29 )
	 : "\0\37\00\37\36\37\36\37\37\36\37\36\37"[m] );
}

static void convert_ymd_j(DATE *d)
{
 int j;
 int y;
 int m;

 y = d->y;
 j = 146097 * (y / 400);
 y %= 400;
 if (y > 0) j ++;
 j += 36524 * (y / 100);
 y %= 100;
 if (y > 0) j --;
 j += 1461 * (y / 4);
 y %= 4;
 if (y > 0) j ++;
 j += 365 * y;
 y = d->y;
 for (m=d->m-1;m>0;m--) j += mdays(y,m);
 d->j = j + d->d - 1;
}

static void convert_j_ymd(DATE *d)
{
 int n;
 int j;

 j = d->j;
 d->y = 400 * (j / 146097); /* 146097d = 400y */
 j %= 146097;
 if (j > 1100)
  { j -= 500;
    n = j / 365;
    j %= 365;
    d->y += n;
    j += 499 + ((n + 99) / 100) - ((n + 3) / 4);
  }
 while (1)
  { n = ydays(d->y);
    if (j < n) break;
    d->y ++;
    j -= n;
  }
 d->m = 1;
 while (1)
  { n = mdays(d->y,d->m);
    if (j < n) break;
    d->m ++;
    j -= n;
    if (d->m > 12) abort();
  }
 d->d = j + 1;
}

static const char *tzstr(int tz)
{
 static char buf[6];
 char *bp;

 bp = &buf[0];
 if (tz < 0)
  { *bp++ = '-';
    tz = - tz;
  }
 else
  { *bp++ = '+';
  }
 sprintf(bp,"%02d",tz/60);
 sprintf(bp+2,"%02d",tz%60);
 return(&buf[0]);
}

static __inline__ int planetset_tst(PLANETSET s, int p)
{
 return((s&(1U<<p))?1:0);
}

#define planetset_empty() ((PLANETSET)0)

#define planetset_full() ((PLANETSET)((1U<<NPLANETS)-1))

static __inline__ PLANETSET planetset_add(PLANETSET s, int p)
{
 return(s|(1U<<p));
}

static __inline__ PLANETSET planetset_del(PLANETSET s, int p)
{
 return(s&~(1U<<p));
}

static __inline__ int planetset_subset(PLANETSET big, PLANETSET small)
{
 return((small&~big)==0);
}

static void find_grands(GRANDASPECT *g)
{
 int pv[GA_MAX_PLANETS];
 PLANETSET ptaken;

 static void check(void)
  { int gan;
    GAGA *ga;
    int a;
    GAASP *gasp;
    double d;
    PLANETSET ps;
    int i;
    for (gan=0;gan<g->ngas;gan++)
     { ga = &g->gas[gan];
       ps = planetset_empty();
       for (a=0;a<ga->nasps;a++)
	{ gasp = &g->asps[a];
	  d = abs_angdiff(planets[pv[gasp->p1]].loc.theta,planets[pv[gasp->p2]].loc.theta);
	  if (fabs(d-gasp->ang) > gasp->orb) goto nextga;
	  ps = planetset_add(planetset_add(ps,pv[gasp->p1]),pv[gasp->p2]);
	}
       for (i=g->nfound-1;i>=0;i--)
	{ if (planetset_subset(g->found[i].planets,ps))
	   { goto nextga;
	   }
	}
       for (i=g->nfound-1;i>=0;i--)
	{ if (planetset_subset(ps,g->found[i].planets))
	   { g->nfound --;
	     if (i < g->nfound) g->found[i] = g->found[g->nfound];
	   }
	}
       g->found = realloc(g->found,(g->nfound+1)*sizeof(GAFOUND));
       g->found[g->nfound].planets = ps;
       g->found[g->nfound].ga = ga;
       g->nfound ++;
nextga:;
     }
  }

 static void enumerate(int x)
  { int i;
    for (i=0;i<NPLANETS;i++)
     { if (planetset_tst(ptaken,i)) continue;
       pv[x] = i;
       if (x > 0)
	{ ptaken = planetset_add(ptaken,i);
	  enumerate(x-1);
	  ptaken = planetset_del(ptaken,i);
	}
       else
	{ check();
	}
     }
  }

 if (g->nplanets > GA_MAX_PLANETS) abort();
 g->nfound = 0;
 g->found = 0;
 ptaken = planetset_empty();
 enumerate(g->nplanets-1);
}

static void find_aspects(void)
{
 int c;
 int i;
 int j;
 int x;
 ASPLIST list[PLANETPAIRS];
 ASPLIST lt;
 double err;

 x = 0;
 for (i=0;i<NPLANETS;i++) for (j=i+1;j<NPLANETS;j++)
  { list[x].i = i;
    list[x].j = j;
    list[x].d = abs_angdiff(planets[i].loc.theta,planets[j].loc.theta);
    x ++;
  }
 for (c=0;c<NASPECTS;c++)
  { do
     { j = 0;
       for (i=1;i<PLANETPAIRS;i++)
	{ if ( fabs(list[i].d-aspects[c].ang) <
	       fabs(list[i-1].d-aspects[c].ang) )
	   { lt = list[i];
	     list[i] = list[i-1];
	     list[i-1] = lt;
	     j ++;
	   }
	}
     } while (j > 0);
    for (i=0;i<PLANETPAIRS;i++)
     { err = fabs(list[i].d-aspects[c].ang);
       if (!allaspects && (err > aspects[c].noorb)) break;
     }
    if (i < PLANETPAIRS) i ++;
    aspects[c].listn = i;
    aspects[c].list = malloc(i*sizeof(ASPLIST));
    bcopy(&list[0],aspects[c].list,i*sizeof(ASPLIST));
  }
 for (i=0;i<NGRANDASPECTS;i++) find_grands(&grandaspects[i]);
}

static void print_grands(GRANDASPECT *g, void (*fn)(GAFOUND *))
{
 int i;
 GAGA *ga;
 int j;
 GAFOUND *f;

 for (i=0;i<g->ngas;i++)
  { ga = &g->gas[i];
    for (j=0;j<g->nfound;j++)
     { f = &g->found[j];
       if (f->ga == ga) (*fn)(f);
     }
  }
}

static void print_planet_list(PLANETSET showvec, void (*fn)(PLANET *, double))
{
 int i;
 int i0;
 int i1;
 int n;
 int v[NPLANETS];
 int px[NPLANETS];
 double d;
 double maxd;
 int maxx;
 double stheta[NPLANETS];

 static void sort_px(void)
  { int m;
    int v;
    do
     { m = 0;
       for (i=1;i<n;i++)
	{ if (stheta[px[i]] < stheta[px[i-1]])
	   { v = px[i];
	     px[i] = px[i-1];
	     px[i-1] = v;
	     m ++;
	   }
	}
     } while (m > 0);
  }

 n = 0;
 for (i=0;i<NPLANETS;i++) if (planetset_tst(showvec,i)) v[n++] = i;
 for (i=0;i<n;i++)
  { px[i] = i;
    stheta[i] = planets[v[i]].loc.theta;
  }
 sort_px();
 maxd = stheta[px[0]] + 360 - stheta[px[n-1]];
 maxx = 0;
 for (i=1;i<n;i++)
  { d = stheta[px[i]] - stheta[px[i-1]];
    if (d > maxd)
     { maxd = d;
       maxx = i;
     }
  }
 for (i=maxx-1;i>=0;i--) stheta[px[i]] += 360;
 sort_px();
 i0 = -1;
 for (i=1;i<n;i++)
  { if (stheta[px[i]]-stheta[px[i-1]] < SIGN_WIDTH)
     { double a0;
       double a1;
       double s0;
       double s1;
       double w;
       i0 = i - 1;
       i1 = i;
       while (1)
	{ s0 = stheta[px[i0]];
	  s1 = stheta[px[i1]];
	  w = (i1 - i0) * SIGN_WIDTH;
	  a0 = (s0 + s1 - w) / 2;
	  a1 = (s0 + s1 + w) / 2;
	  if ((i0 > 0) && (a0-stheta[px[i0-1]] < SIGN_WIDTH))
	   { i0 --;
	     continue;
	   }
	  if ((i1 < n-1) && (stheta[px[i1+1]]-a1 < SIGN_WIDTH))
	   { i1 ++;
	     continue;
	   }
	  break;
	}
       for (i=i0;i<=i1;i++)
	{ stheta[px[i]] = a0;
	  a0 += SIGN_WIDTH;
	}
     }
  }
 for (i=0;i<n;i++) (*fn)(&planets[v[i]],stheta[i]);
}

static void ps_print_planet(PLANET *p, double stheta, const char *psfn)
{
 printf("%.9f %.9f /%s %s %s\n",
	p->loc.theta,
	stheta,
	p->name,
	(p->vel.theta<0)?"true":"false",
	psfn );
}
#define PS_PRINT_PLANET(s) ({ static void foo(PLANET *p, double sa)	\
			       { ps_print_planet(p,sa,s);		\
			       }					\
			      &foo;					\
			    })

static void write_ps_output(void)
{
 int c;
 int i;
 int didbar;
 double err;
 ASPECT *a;
 ASPLIST *al;
 PLANETSET pshow;

 printf("/colourful %s def\n",colourful?"true":"false");
 printf("gsave 400 575 zodiac\n");
 printf("(%.6f) showtime1\n(%04d-%02d-%02d %02d:%02d:%02d %s) showtime2\n",
	jtime,
	ce_date.y,ce_date.m,ce_date.d,
	ce_time.h,ce_time.m,ce_time.s,
	tzstr(ce_tz));
 printf("newsigns\n");
 print_planet_list(planetset_full(),PS_PRINT_PLANET("mainplanet"));
 for (c=0;c<NASPECTS;c++)
  { a = &aspects[c];
    if (fabs(a->list[0].d-a->ang) > a->weakorb) continue;
    for (i=0,al=a->list;i<a->listn;i++,al++)
     { err = fabs(al->d-a->ang);
       if (err <= a->strongorb)
	{ printf("%.9f %.9f /%s strongaspect\n",planets[al->i].loc.theta,planets[al->j].loc.theta,a->name);
	}
       else if (err <= a->weakorb)
	{ printf("%.9f %.9f /%s weakaspect\n",planets[al->i].loc.theta,planets[al->j].loc.theta,a->name);
	}
     }
  }
 printf("grestore\n");
 for (c=0;c<NASPECTS;c++)
  { a = &aspects[c];
    if (fabs(a->list[0].d-a->ang) > a->weakorb) continue;
    printf("gsave /%s aspzodiac\n",a->name);
    pshow = planetset_empty();
    for (i=0,al=a->list;i<a->listn;i++,al++)
     { err = fabs(al->d-a->ang);
       if (err <= a->strongorb)
	{ printf("%.9f %.9f /%s strongaspect\n",planets[al->i].loc.theta,planets[al->j].loc.theta,a->name);
	  pshow = planetset_add(planetset_add(pshow,al->i),al->j);
	}
       else if (err <= a->weakorb)
	{ printf("%.9f %.9f /%s weakaspect\n",planets[al->i].loc.theta,planets[al->j].loc.theta,a->name);
	  pshow = planetset_add(planetset_add(pshow,al->i),al->j);
	}
     }
    print_planet_list(pshow,PS_PRINT_PLANET("aspplanet"));
    printf("grestore\n");
  }
 printf("showpage\n");
 for (c=0;c<NASPECTS;c++)
  { a = &aspects[c];
    printf("/%s newaspects\n",a->name);
    didbar = 0;
    for (i=0,al=a->list;i<a->listn;i++,al++)
     { err = fabs(al->d-a->ang);
       if ((err > a->strongorb) && (didbar < 1))
	{ printf("aspectbar\n");
	  didbar = 1;
	}
       if ((err > a->weakorb) && (didbar < 2))
	{ printf("aspectbar\n");
	  didbar = 2;
	}
       if (allaspects && (err > a->noorb) && (didbar < 3))
	{ printf("aspectbar\n");
	  didbar = 3;
	}
       printf("/%s /%s %.9f aspect\n",planets[al->i].name,planets[al->j].name,err);
     }
    printf("endaspects\n");
  }
 printf("/grand newaspects\n");
 for (i=0;i<NGRANDASPECTS;i++)
  { print_grands(&grandaspects[i],({ static void foo(GAFOUND *f)
				      { int np;
					int k;
					printf("(%s) ",f->ga->name);
					np = 0;
					for (k=NPLANETS-1;k>=0;k--)
					 { if (planetset_tst(f->planets,k))
					    { printf("/%s ",planets[k].name);
					      np ++;
					    }
					 }
					printf("%d grandaspect\n",np);
				      }
				     &foo;
				   }));
  }
 printf("endaspects\n");
 printf("showpage\n");
}

static void write_text_output(void)
{
 int c;
 int i;
 int didbar;
 double err;
 ASPECT *a;
 ASPLIST *al;

 printf("%04d-%02d-%02d %02d:%02d:%02d %s (ephemeris time %.6f)\n\n",
	ce_date.y,ce_date.m,ce_date.d,
	ce_time.h,ce_time.m,ce_time.s,
	tzstr(ce_tz),
	jtime);
 print_planet_list(planetset_full(),({ static void foo(PLANET *p, double sa __attribute__((__unused__)))
					{ double a;
					  int s;
					  a = p->loc.theta;
					  while (a < 0) a += 360;
					  while (a >= 360) a -= 360;
					  for (s=NSIGNS-1;(s>0)&&(signs[s].start>a);s--) ;
					  printf("%s\t%8.5f° %s%s\n",
						p->name,
						a-signs[s].start,
						signs[s].name,
						(p->vel.theta<0)?", retrograde":"");
					}
				       &foo;
				     }));
 for (c=0;c<NASPECTS;c++)
  { a = &aspects[c];
    printf("\n");
    didbar = 0;
    for (i=0,al=a->list;i<a->listn;i++,al++)
     { err = fabs(al->d-a->ang);
       if ((err > a->strongorb) && (didbar < 1))
	{ printf("%s --------\n",a->name);
	  didbar = 1;
	}
       if ((err > a->weakorb) && (didbar < 2))
	{ printf("%s --------\n",a->name);
	  didbar = 2;
	}
       printf("%s %.5f %s %s\n",a->name,err,planets[al->i].name,planets[al->j].name);
     }
  }
 didbar = 0;
 for (i=0;i<NGRANDASPECTS;i++)
  { print_grands(&grandaspects[i],({ static void foo(GAFOUND *f)
				      { int k;
					if (! didbar)
					 { printf("\n");
					   didbar = 1;
					 }
					printf("%s:",f->ga->name);
					for (k=NPLANETS-1;k>=0;k--)
					 { if (planetset_tst(f->planets,k))
					    { printf(" %s",planets[k].name);
					    }
					 }
					printf("\n");
				      }
				     &foo;
				   }));
  }
}

static void parse_tzstr(const char *tzstr)
{
 char sign;
 int zh;
 int zm;

 if (tzstr == 0)
  { ce_tz = 0;
    return;
  }
 sscanf(tzstr,"%c%2d%2d",&sign,&zh,&zm);
 ce_tz = (zh * 60) + zm;
 if (sign == '-') ce_tz = - ce_tz;
}

static void set_ce_time(int s)
{
 ce_time.h = s / 3600;
 ce_time.m = (s / 60) % 60;
 ce_time.s = s % 60;
}

static void setup_jd(char *jtstr, char *tz)
{
 int subday;
 char *dot;

 dot = index(jtstr,'.');
 if (dot == 0)
  { ce_date.j = atoi(jtstr);
    subday = 0;
  }
 else
  { *dot = '\0';
    ce_date.j = atoi(jtstr);
    *dot = '.';
    subday = floor((SECPERDAY*atof(dot))+.5);
  }
 while (subday < 0)
  { subday += SECPERDAY;
    ce_date.j --;
  }
 while (subday >= SECPERDAY)
  { subday -= SECPERDAY;
    ce_date.j ++;
  }
 jtime = ce_date.j + (subday / (double)SECPERDAY);
 parse_tzstr(tz);
 subday += (ce_tz * 60) + JOFFSET_SEC;
 while (subday >= SECPERDAY)
  { subday -= SECPERDAY;
    ce_date.j ++;
  }
 ce_date.j -= JOFFSET_DAY;
 if (ce_date.j < 0)
  { fprintf(stderr,"%s: jdate out of range\n",__progname);
    exit(1);
  }
 convert_j_ymd(&ce_date);
 set_ce_time(subday);
}

static void setup_ymdhms(char *yr, char *mo, char *dy, char *hr, char *mn, char *sc, char *tz)
{
 int subday;

 ce_date.y = atoi(yr);
 ce_date.m = atoi(mo);
 ce_date.d = atoi(dy);
 convert_ymd_j(&ce_date);
 subday = (atoi(hr) * 3600) + (atoi(mn) * 60) + atoi(sc);
 while (subday < 0)
  { subday += SECPERDAY;
    ce_date.j --;
  }
 while (subday >= SECPERDAY)
  { subday -= SECPERDAY;
    ce_date.j ++;
  }
 set_ce_time(subday);
 subday -= JOFFSET_SEC;
 parse_tzstr(tz);
 subday -= ce_tz * 60;
 while (subday < 0)
  { subday += SECPERDAY;
    ce_date.j --;
  }
 jtime = ce_date.j + JOFFSET_DAY + (subday / (double)SECPERDAY);
}

static void setup_locvar(char *s)
{
 char *e;
 char c;
 int neg;
 double v;
 int gotdot;
 int bit;
 double *loc;
 const char *what;

 gotdot = 0;
 for (e=s;*e;e++)
  { switch (*e)
     { case '.':
	  gotdot = 1;
	  /* fall through */
       case '0': case '1': case '2': case '3': case '4':
       case '5': case '6': case '7': case '8': case '9':
	  continue;
	  break;
     }
    break;
  }
 if ((e == s) || (gotdot && (e == s+1)))
  { fprintf(stderr,"%s: %s: missing number\n",__progname,s);
    exit(1);
  }
 else if (! *e)
  { fprintf(stderr,"%s: %s: missing direction\n",__progname,s);
    exit(1);
  }
 switch (*e)
  { case 'n': case 'N': neg = 0; loc = &latitude;  what = "latitude";  bit = LG_LAT;  break;
    case 's': case 'S': neg = 1; loc = &latitude;  what = "latitude";  bit = LG_LAT;  break;
    case 'e': case 'E': neg = 0; loc = &longitude; what = "longitude"; bit = LG_LONG; break;
    case 'w': case 'W': neg = 1; loc = &longitude; what = "longitude"; bit = LG_LONG; break;
    default:
       fprintf(stderr,"%s: %s: invalid direction\n",__progname,s);
       exit(1);
       break;
  }
 if (gotdot)
  { if (e[1])
     { fprintf(stderr,"%s: %s: junk after %s\n",__progname,s,what);
       exit(1);
     }
    c = e[-1];
    e[-1] = '\0';
    v = atof(s);
    e[-1] = c;
  }
 else
  { c = e[-1];
    e[-1] = '\0';
    v = atoi(s);
    e[-1] = c;
    v += atof(e) / 60;
  }
 *loc = neg ? -v : v;
 locgot |= bit;
}

static void setup_location(char *str1, char *str2)
{
 locgot = 0;
 setup_locvar(str1);
 setup_locvar(str2);
 if (locgot != (LG_LAT|LG_LONG))
  { fprintf(stderr,"%s: %s %s: need one each latitude and longitude\n",__progname,str1,str2);
    exit(1);
  }
}

static void surface_offset(void)
{
/*
 double x;
 double y;
 double z;
 double theta;

 a day is approximately 
 theta = longitude + (2*M_PI*fmod(jtime,1));
 x = EARTH_RADIUS * 
*/
 surface_loc.x = 0;
 surface_loc.y = 0;
 surface_loc.z = 0;
 surface_vel.x = 0;
 surface_vel.y = 0;
 surface_vel.z = 0;
}

static void compute_sidereal(void)
{
 double solar_angle;
 int ss;

 solar_angle = atan2(planets[PLANET_SUN].loc.e.y,planets[PLANET_SUN].loc.e.x) - ZERO_ARIES;
 if (! textmode) printf("%% solar_angle = %26.21f\n",solar_angle);
 stime = (fmod(jtime-.5,1) * 360)
	+ longitude
	+ (solar_angle * RADTODEG)
	- 180;
 if (! textmode) printf("%% stime = %g -> ",stime);
 while (stime < 0) stime += 360;
 while (stime > 360) stime -= 360;
 ss = stime * 240;
 if (ss > 86400) ss -= 86400;
 if (! textmode) printf("%g (%02d:%02d:%02d)\n",stime,ss/3600,(ss/60)%60,ss%60);
}

static void ephemeris_compute(const char *eph_file)
{
 int i;

 h = eph_open(eph_file);
 if (h == 0)
  { fprintf(stderr,"%s: can't open %s: %s\n",__progname,eph_file,strerror(errno));
    exit(1);
  }
 compute_ecliptic();
 for (i=0;i<NPLANETS;i++) getloc(i);
 compute_sidereal();
 surface_offset();
 for (i=0;i<NPLANETS;i++) do_spherical(&planets[i]);
}

static void setup_data_one(PLANET *p, const char *sname, const char *angle)
{
 int i;
 int sign;
 int sl;
 int retrograde;
 double a;
 char *aend;

 sign = -1;
 sl = strlen(sname);
 for (i=0;i<NSIGNS;i++)
  { if (!strncasecmp(signs[i].name,sname,sl))
     { if (sign < 0)
	{ sign = i;
	}
       else
	{ fprintf(stderr,"%s: sign `%s' is ambiguous\n",__progname,sname);
	  exit(1);
	}
     }
  }
 if (sign < 0)
  { fprintf(stderr,"%s: sign `%s' is unknown\n",__progname,sname);
    exit(1);
  }
 a = strtod(angle,&aend);
 if (aend == angle)
  { fprintf(stderr,"%s: invalid angle `%s'\n",__progname,angle);
    exit(1);
  }
 retrograde = 0;
 if ((*aend == 'r') || (*aend == 'R'))
  { retrograde = 1;
    aend ++;
  }
 if (*aend)
  { fprintf(stderr,"%s: junk after angle in `%s'\n",__progname,angle);
    exit(1);
  }
 if ((a < -360) || (a > 360))
  { fprintf(stderr,"%s: angle `%s' out of range\n",__progname,angle);
    exit(1);
  }
 p->loc.theta = signs[i].start + a;
 p->vel.theta = retrograde ? -1 : 1;
}

static void setup_data(char **args)
{
 int i;

 for (i=0;i<NPLANETS;i++) setup_data_one(&planets[i],args[i+i],args[i+i+1]);
}

int main(int, char **);
int main(int ac, char **av)
{
 int i;

 for (;ac>2;ac--,av++)
  { if (!strcmp(av[1],"-text"))
     { textmode = 1;
       continue;
     }
    if (!strcmp(av[1],"-all"))
     { allaspects = 1;
       continue;
     }
    if (!strcmp(av[1],"-colour") || !strcmp(av[1],"-color"))
     { colourful = 1;
       continue;
     }
    break;
  }
 switch (ac)
  { case 5:
    case 6:
       setup_location(av[2],av[3]);
       setup_jd(av[4],av[5]);
       ephemeris_compute(av[1]);
       break;
    case 10:
    case 11:
       setup_location(av[2],av[3]);
       setup_ymdhms(av[4],av[5],av[6],av[7],av[8],av[9],av[10]);
       ephemeris_compute(av[1]);
       break;
    case 21:
       setup_data(av+1);
       break;
    default:
       fprintf(stderr,"Usage: %s ephemeris-file latitude longitude jdate [[+-]zzzz]\n",__progname);
       fprintf(stderr,"       %s ephemeris-file latitude longitude yyyy mm dd hh mm ss [[+-]zzzz]\n",__progname);
       fprintf(stderr,"       %s <data>\n",__progname);
       fprintf(stderr,"\
Latitude must have a suffix N or S; longitude E or W.\n\
<data> consists of a sign name and an angle, for each planet.  Sign names\n\
may be abbreviated to any unambiguous abbrevation; an R (or r) suffix on\n\
the angle indicates retrograde motion.  Planets must be given in order:\n\
   ");
       for (i=0;i<NPLANETS;i++) fprintf(stderr," %s",planets[i].name);
       fprintf(stderr,"\n");
       fprintf(stderr,"Possible options (immediately after %s):\n",__progname);
       fprintf(stderr,"-text\n\tProduce text output (default: PostScript)\n");
       fprintf(stderr,"-all\n\tList all aspects regardless of orb (default: reasonable cutoff)\n");
       fprintf(stderr,"-colour\n\tProduce colour output (default: B/W; ignored if -text)\n");
       exit(1);
       break;
  }
 find_aspects();
 if (textmode) write_text_output(); else write_ps_output();
 exit(0);
}