#include <cg6.h>
#include <poll.h>
#include <math.h>
#include <stdio.h>
#include <ctype.h>
#include <fcntl.h>
#include <errno.h>
#include <stdlib.h>
#include <signal.h>
#include <unistd.h>
#include <stdarg.h>
#include <strings.h>
#include <termios.h>
#include <sys/mman.h>
#include <sys/ioctl.h>
#include <machine/kbd.h>
#include <machine/fbio.h>
#include <machine/kbio.h>
#include <machine/vuid_event.h>

extern const char *__progname;

static const char *fbpath = "/dev/cgsix0";
static const char *kbpath = "/dev/kbd1";
static const char *model = 0;

static int fbfd;
static volatile struct cg6thc *thc;
static volatile struct cg6tec *tec;
static volatile struct cg6fbc *fbc;
static volatile struct brooktree *bt;
static volatile unsigned char *vram;
static int drawbit;
#define DT_NONE   1
#define DT_NORMAL 2
static int disptype = DT_NORMAL;

typedef struct xy XY;
typedef struct xyz XYZ;
typedef struct xyzw XYZW;

struct xy {
  double x;
  double y;
  } ;

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

struct xyzw {
  double x;
  double y;
  double z;
  double w;
  } ;

typedef struct color COLOR;
typedef struct line LINE;
typedef struct face FACE;
typedef struct pt PT;
typedef struct obj OBJ;
typedef struct arc ARC;
typedef struct objlist OBJLIST;
typedef struct xyzn XYZN;

struct xyzn {
  XYZ v;
  XYZ u;
  } ;

struct obj {
  unsigned char type;
#define OT_LINE 1
#define OT_FACE 2
  unsigned int flags;
#define OF_MARK 0x00000001
  union {
    LINE *l;
    FACE *f;
    } u;
  int x;
  ARC *out;
  ARC *in;
  } ;

struct arc {
  ARC *flink;
  ARC *tlink;
  OBJ *f;
  OBJ *t;
  } ;

struct objlist {
  OBJLIST *link;
  OBJ *o;
  } ;

struct color {
  COLOR *link;
  char *name;
  int namelen;
  unsigned char r;
  unsigned char g;
  unsigned char b;
  unsigned char inx;
  unsigned int flags;
#define CF_USED 0x00000001
  } ;

struct pt {
  PT *link;
  unsigned int flags;
#define PF_PROJECTED 0x00000001
  char *name;
  int namelen;
  XYZ xyz;
  XY scr[2];
  } ;

struct line {
  LINE *link;
  char *name;
  int namelen;
  PT *a;
  PT *b;
  COLOR *color;
  } ;

struct face {
  FACE *link;
  char *name;
  int namelen;
  PT *a;
  PT *b;
  PT *c;
  COLOR *color;
  XYZW plane;
  } ;

static PT *pts;
static LINE *lines;
static FACE *faces;
static COLOR *colors;
static void (*linefail)(const char *, va_list *)
				__attribute__((__noreturn__));
#define CMFIXED_BLACK 0
#define CMFIXED_WHITE 1
static int nextcx = 2;
static unsigned char cmap[16][3];
static int nobjs;
static OBJ **objs_l;
static OBJ **objs_r;
static signed char *occl1d;
static signed char **occl;
static unsigned int wantsort = 0;
#define SORT_ONCE   0x00000001
#define SORT_ALWAYS 0x00000002
#define SORT_DUMP   0x00000004
static OBJ **objv;
static int ex;
static XYZ ea;

static int kbdfd;
#define EVQ_SIZE 64 /* must be a power of two */
static int evq[EVQ_SIZE];
static int evqh;
static int evqt;
#define EV_TYPE  0x000f0000
#define EVT_NIL  0x00000000
#define EVT_KBD  0x00010000
#define EVK_CODE 0x0000007f /* must be low bits */
#define EVK_UPDN 0x00000080
#define EVK_UP   0x00000080
#define EVK_DN   0x00000000
#define EVT_TICK 0x00020000

static int curtime;
static volatile int rawtime;
static int timerunning;

static const XYZ xyz0 = { 0, 0, 0 };

static XYZ up;
static XYZ out;
static XYZ rt;
static XYZN lea;
static XYZ ler;
static XYZN rea;
static XYZ rer;
static double ed;
static double es;
static double ipd;
static int mvbits;
#define MV_UP  0x00000001
#define MV_LF  0x00000002
#define MV_RT  0x00000004
#define MV_DN  0x00000008
#define MV_CW  0x00000010
#define MV_CCW 0x00000020
#define MV_IN  0x00000040
#define MV_OUT 0x00000080
#define MV_AP  0x00000100
#define MV_TG  0x00000200
#define MV_ZI  0x00000400
#define MV_ZO  0x00000800
static int wantrender;

static void setcmap(int bit)
{
 int i;
 int v;

 switch (disptype)
  { default: abort(); break;
    case DT_NONE:
       break;
    case DT_NORMAL:
       bt->addr = 0;
       for (i=0;i<256;i++)
	{ v = bit ? (i >> 4) : (i & 15);
	  bt->cmap = cmap[v][0] * 0x01000000;
	  bt->cmap = cmap[v][1] * 0x01000000;
	  bt->cmap = cmap[v][2] * 0x01000000;
	}
       break;
  }
}

static void resethw(void)
{
 switch (disptype)
  { default: abort(); break;
    case DT_NONE:
       break;
    case DT_NORMAL:
       thc->cursxy = 0xffe0ffe0;
       tec->mv = 0;
       tec->clip = 0;
       tec->vdc = 0;
       fbc->s = 0;
       fbc->bg = 0;
       fbc->pixelm = ~0;
       fbc->s = 0;
       fbc->alu = ALU_NORMAL | ALU_FG;
       fbc->mode = MODE_NORMAL | MODE_COLOR8;
       fbc->clip = 0;
       fbc->offx = 0;
       fbc->offy = 0;
       fbc->clipminx = 0;
       fbc->clipminy = 0;
       fbc->clipmaxx = 1151;
       fbc->clipmaxy = 899;
       fbc->fg = 0;
       fbc->pm = ~0;
       break;
  }
}

static void initfb(void)
{
 void *mrv;
 struct fbgattr a;

 switch (disptype)
  { default: abort(); break;
    case DT_NONE:
       break;
    case DT_NORMAL:
       fbfd = open(fbpath,O_RDWR,0);
       if (fbfd < 0)
	{ fprintf(stderr,"%s: can't open %s: %s\n",__progname,fbpath,strerror(errno));
	  exit(1);
	}
       if ( (ioctl(fbfd,FBIOGATTR,&a) < 0) ||
	    (a.fbtype.fb_type != FBTYPE_SUNFAST_COLOR) ||
	    (a.fbtype.fb_width != 1152) ||
	    (a.fbtype.fb_height != 900) ||
	    (a.fbtype.fb_depth != 8) )
	{ fprintf(stderr,"%s: %s isn't an 1152x900 cgsix\n",__progname,fbpath);
	  exit(1);
	}
       mrv = mmap(0,0x16000+a.fbtype.fb_size,PROT_READ|PROT_WRITE,MAP_SHARED,fbfd,0x70000000);
       if (mrv == MAP_FAILED)
	{ fprintf(stderr,"%s: can't mmap %s: %s\n",__progname,fbpath,strerror(errno));
	  exit(1);
	}
       fbc = mrv;
       thc = (void *)(0x5000+(unsigned char *)mrv);
       tec = (void *)(0x1000+(unsigned char *)mrv);
       bt = (void *)(0x2000+(unsigned char *)mrv);
       vram = 0x16000 + (unsigned char *)mrv;
       setcmap(1);
       drawbit = 0;
       resethw();
       fbc->fg = CMFIXED_BLACK * 0x11;
       fbc->pm = ~0;
       fbc->arecty = 0;
       fbc->arectx = 0;
       fbc->arecty = 900;
       fbc->arectx = 1152;
       while (fbc->s & 0x10000000) ;
       fbc->clipminx = 0;
       fbc->clipminy = 0;
       fbc->clipmaxx = 1151;
       fbc->clipmaxy = 899;
       fbc->pm = 0x0f;
       break;
  }
}

static void initkbd(void)
{
 int type;
 int mode;

 kbdfd = open(kbpath,O_RDWR,0);
 if (kbdfd < 0)
  { fprintf(stderr,"%s: can't open %s: %s\n",__progname,kbpath,strerror(errno));
    exit(1);
  }
 if ( (ioctl(kbdfd,KIOCTYPE,&type) < 0) ||
      (type != KB_SUN3) )
  { close(kbdfd);
    fprintf(stderr,"%s: %s isn't a type-3\n",__progname,kbpath);
    exit(1);
  }
 mode = 1;
 ioctl(kbdfd,KIOCSDIRECT,&mode);
 mode = TR_UNTRANS_EVENT;
 ioctl(kbdfd,KIOCTRANS,&mode);
 mode = 1;
 ioctl(kbdfd,FIONBIO,&mode);
}

static void initevq(void)
{
 evqh = 0;
 evqt = 0;
}

static void raw_tick(int sig __attribute__((__unused__)))
{
 rawtime ++;
}

static void starttiming(void)
{
 struct itimerval itv;

 if (timerunning) return;
 timerunning = 1;
 itv.it_value.tv_sec = 0;
 itv.it_value.tv_usec = 10000;
 itv.it_interval = itv.it_value;
 setitimer(ITIMER_REAL,&itv,0);
}

static void stoptiming(void)
{
 struct itimerval itv;

 if (! timerunning) return;
 timerunning = 0;
 itv.it_value.tv_sec = 0;
 itv.it_value.tv_usec = 0;
 itv.it_interval = itv.it_value;
 setitimer(ITIMER_REAL,&itv,0);
}

static void inittiming(void)
{
 struct sigaction sa;

 curtime = 0;
 rawtime = 0;
 sa.sa_handler = raw_tick;
 sigemptyset(&sa.sa_mask);
 sa.sa_flags = SA_RESTART;
 sigaction(SIGALRM,&sa,0);
 timerunning = 1;
 stoptiming();
}

static void queue_event(int ev)
{
 if (((evqh-evqt) & (EVQ_SIZE-1)) == (EVQ_SIZE-1))
  { write(1,"evq full\n",9);
  }
 else
  { evq[evqh] = ev;
    evqh = (evqh+1) & (EVQ_SIZE-1);
  }
}

static int dequeue_event(void)
{
 int ev;

 if (evqh == evqt) return(EVT_NIL);
 ev = evq[evqt];
 evqt = (evqt+1) & (EVQ_SIZE-1);
 return(ev);
}

static int checktick(void)
{
 int now;

 now = rawtime;
 if (curtime != now)
  { queue_event(EVT_TICK);
    curtime ++;
    return(1);
  }
 return(0);
}

static void fbline(int x1, int y1, int x2, int y2, int z)
{
 switch (disptype)
  { default: abort(); break;
    case DT_NONE:
       break;
    case DT_NORMAL:
       fbc->fg = z * 0x11;
       fbc->aliney = y1;
       fbc->alinex = x1;
       fbc->aliney = y2;
       fbc->alinex = x2;
       while ((fbc->draw & 0xa0000000) == 0xa0000000) ;
       break;
  }
}

static void fbtri(int x1, int y1, int x2, int y2, int x3, int y3, int z)
{
 switch (disptype)
  { default: abort(); break;
    case DT_NONE:
       break;
    case DT_NORMAL:
       fbc->fg = z * 0x11;
       fbc->atriy = y1;
       fbc->atrix = x1;
       fbc->atriy = y2;
       fbc->atrix = x2;
       fbc->atriy = y3;
       fbc->atrix = x3;
       while ((fbc->draw & 0xa0000000) == 0xa0000000) ;
       break;
  }
}

static void fbdrain(void)
{
 switch (disptype)
  { default: abort(); break;
    case DT_NONE:
       break;
    case DT_NORMAL:
       while (fbc->s & 0x10000000) ;
       break;
  }
}

static void fbclear(void)
{
 switch (disptype)
  { default: abort(); break;
    case DT_NONE:
       break;
    case DT_NORMAL:
       fbc->fg = CMFIXED_BLACK * 0x11;
       fbc->arecty = 0;
       fbc->arectx = 0;
       fbc->arecty = 900;
       fbc->arectx = 1152;
       while ((fbc->draw & 0xa0000000) == 0xa0000000) ;
       break;
  }
}

static void setpm(unsigned int pm)
{
 switch (disptype)
  { default: abort(); break;
    case DT_NONE:
       break;
    case DT_NORMAL:
       fbc->pm = pm;
       break;
  }
}

static double sqr(double) __attribute__((__const__));
static double sqr(double v)
{
 return(v*v);
}

static XY xysub(XY, XY) __attribute__((__const__));
static XY xysub(XY a, XY b)
{
 return((XY){.x=a.x-b.x,.y=a.y-b.y});
}

static double xyzlen(XYZ) __attribute__((__const__));
static double xyzlen(XYZ v)
{
 return(hypot(hypot(v.x,v.y),v.z));
}

static double xyzlen2(XYZ) __attribute__((__const__));
static double xyzlen2(XYZ v)
{
 return((v.x*v.x)+(v.y*v.y)+(v.z*v.z));
}

static XYZ xyzscale(XYZ, double) __attribute__((__const__));
static XYZ xyzscale(XYZ v, double s)
{
 return((XYZ){.x=v.x*s,.y=v.y*s,.z=v.z*s});
}

static XYZ norm(XYZ) __attribute__((__const__));
static XYZ norm(XYZ v)
{
 double l;

 l = xyzlen(v);
 if (l < 1e-20) return(xyz0);
 return(xyzscale(v,1/l));
}

static XYZ cross(XYZ, XYZ) __attribute__((__const__));
static XYZ cross(XYZ a, XYZ b)
{
 return((XYZ) { .x = (a.y * b.z) - (a.z * b.y),
		.y = (a.z * b.x) - (a.x * b.z),
		.z = (a.x * b.y) - (a.y * b.x) });
}

static double dot(XYZ, XYZ) __attribute__((__const__));
static double dot(XYZ a, XYZ b)
{
 return((a.x*b.x)+(a.y*b.y)+(a.z*b.z));
}

static XYZ xyzsub(XYZ, XYZ) __attribute__((__const__));
static XYZ xyzsub(XYZ a, XYZ b)
{
 return((XYZ){.x=a.x-b.x,.y=a.y-b.y,.z=a.z-b.z});
}

static XYZ xyzadd(XYZ, XYZ) __attribute__((__const__));
static XYZ xyzadd(XYZ a, XYZ b)
{
 return((XYZ){.x=a.x+b.x,.y=a.y+b.y,.z=a.z+b.z});
}

static XYZ xyzalpha(XYZ, double, XYZ) __attribute__((__const__));
static XYZ xyzalpha(XYZ a, double p, XYZ b)
{
 return((XYZ){ .x = (a.x * p) + (b.x * (1-p)),
	       .y = (a.y * p) + (b.y * (1-p)),
	       .z = (a.z * p) + (b.z * (1-p)) });
}

static XYZ rotunit(XYZ, XYZ, double) __attribute__((__const__));
static XYZ rotunit(XYZ v, XYZ axis, double angle)
{
 return(xyzadd(xyzscale(v,cos(angle)),xyzscale(norm(cross(axis,v)),sin(angle))));
}

static XY project(XYZ, XYZN, XYZ) __attribute__((__const__));
static XY project(XYZ p, XYZN ea, XYZ er)
{
 double d;

 p = xyzsub(p,ea.v);
 d = dot(p,ea.u);
 return((XY){.x=dot(p,er)/d,.y=dot(p,up)/d});
}

static const char *objpstr(OBJ *o)
{
 static char *bufs[8];
 static int hand;
 char *s;

 switch (o->type)
  { default: abort(); break;
    case OT_LINE:
	{ LINE *l;
	  l = o->u.l;
	  if (l->name)
	   { asprintf(&s,"line %s",l->name);
	   }
	  else
	   { asprintf(&s,"line(%p)",(void *)o);
	   }
	}
       break;
    case OT_FACE:
	{ FACE *f;
	  f = o->u.f;
	  if (f->name)
	   { asprintf(&s,"face %s",f->name);
	   }
	  else
	   { asprintf(&s,"face(%p)",(void *)o);
	   }
	}
       break;
  }
 hand --;
 if (hand < 0) hand = 7;
 free(bufs[hand]);
 bufs[hand] = s;
 return(s);
}

static void drawobj(OBJ *o, int ex, int xoff)
{
 switch (o->type)
  { default: abort(); break;
    case OT_LINE:
       fbline( 288 + (o->u.l->a->scr[ex].x*es) + xoff,
	       450 - (o->u.l->a->scr[ex].y*es),
	       288 + (o->u.l->b->scr[ex].x*es) + xoff,
	       450 - (o->u.l->b->scr[ex].y*es),
	       o->u.l->color->inx );
       break;
    case OT_FACE:
       fbtri( 288 + (o->u.f->a->scr[ex].x*es) + xoff,
	      450 - (o->u.f->a->scr[ex].y*es),
	      288 + (o->u.f->b->scr[ex].x*es) + xoff,
	      450 - (o->u.f->b->scr[ex].y*es),
	      288 + (o->u.f->c->scr[ex].x*es) + xoff,
	      450 - (o->u.f->c->scr[ex].y*es),
	      o->u.f->color->inx );
       break;
  }
}

static void project_pts(void)
{
 int i;
 OBJ *o;
 PT *p;

 for (p=pts;p;p=p->link) p->flags &= ~PF_PROJECTED;
 for (i=nobjs-1;i>=0;i--)
  { static void proj(PT *p)
     { if (p->flags & PF_PROJECTED) return;
       p->scr[0] = project(p->xyz,lea,ler);
       p->scr[1] = project(p->xyz,rea,rer);
       p->flags |= PF_PROJECTED;
     }
    o = objs_l[i];
    switch (o->type)
     { default: abort(); break;
       case OT_LINE:
	  proj(o->u.l->a);
	  proj(o->u.l->b);
	  break;
       case OT_FACE:
	  proj(o->u.f->a);
	  proj(o->u.f->b);
	  proj(o->u.f->c);
	  break;
     }
  }
}

/* line-intersect-line */
typedef struct lil LIL;
struct lil {
  int inrange;
  double A;
  double B;
  double D;
  } ;

/*
 * Solving for the point of intersection is easy, mathematically.
 *
 * In METAFONT terms, A[a1,a2] = B[b1,b2], which turns into
 *
 *	A a1.x + (1-A) a2.x = B b1.x + (1-B) b2.x
 *	A a1.y + (1-A) a2.y = B b1.y + (1-B) b2.y
 *
 * Collecting terms and rearranging,
 *
 *	A (a1.x - a2.x) + a2.x = B (b1.x - b2.x) + b2.x
 *	A (a1.y - a2.y) + a2.y = B (b1.y - b2.y) + b2.y
 *
 *	A (a1.x-a2.x) + B (b2.x-b1.x) = b2.x-a2.x
 *	A (a1.y-a2.y) + B (b2.y-b1.y) = b2.y-a2.y
 *
 * Solving A c + B d = e and A f + B g = h gives
 *
 *	A = (e g - d h) / (c g - d f)
 *	B = (c h - e f) / (c g - d f)
 *
 * This is not quite suitable for numeric use, though.  If the matrix
 *  is singular or nearly so, the division will blow up.  To avoid
 *  these problems, we actually return eg-dh, ch-ef, and eg-df in three
 *  different variables (A, B, and D, respectively).  We also negate
 *  all three of them if necessary so that D>=0.
 *
 * Note that A in the return value corresponds to the first two
 *  arguments and B to the second.
 */
static LIL line_intersect_line(XY, XY, XY, XY) __attribute__((__const__));
static LIL line_intersect_line(XY l1p1, XY l1p2, XY l2p1, XY l2p2)
{
 double A;
 double B;
 double c;
 double d;
 double e;
 double f;
 double g;
 double h;
 double D;

 c = l1p1.x - l1p2.x;
 d = l2p2.x - l2p1.x;
 e = l2p2.x - l1p2.x;
 f = l1p1.y - l1p2.y;
 g = l2p2.y - l2p1.y;
 h = l2p2.y - l1p2.y;
 D = (c * g) - (d * f);
 A = (e * g) - (d * h);
 B = (c * h) - (e * f);
 /* If D<0, then A/D>1 turns into A<D; if D>0, A>D.
    To avoid trouble, negate all three variables in this case. */
 if (D < 0)
  { D = - D;
    A = - A;
    B = - B;
  }
 return((LIL){ .A = A, .B = B, .D = D,
	       .inrange = (A > 0) && (B > 0) && (A < D) && (B < D) });
}

/*
 * Conceptually, to see if a point P is inside a triangle with points
 *  A, B, and C, we compute (A-P)×(B-P), (B-P)×(C-P), and (C-P)×(A-P),
 *  where P, A, B, and C are considered as three-vectors with zero Z
 *  components.  The point is inside the triangle iff all cross-product
 *  Z components have the same sign.
 *
 * We don't actually call cross(), because that is rather wasteful for
 *  vectors with known-zero Z components.  Instead, we just compute the
 *  Z component by hand.
 *
 * The triangles are considered open: points exactly on an edge are
 *  considered outside.
 */
static int pt_inside_triangle(XY, XY, XY, XY) __attribute__((__const__));
static int pt_inside_triangle(XY p, XY ta, XY tb, XY tc)
{
 double zab;
 double zbc;
 double zca;

 ta = xysub(ta,p);
 tb = xysub(tb,p);
 tc = xysub(tc,p);
 zab = (ta.x * tb.y) - (ta.y * tb.x);
 zbc = (tb.x * tc.y) - (tb.y * tc.x);
 zca = (tc.x * ta.y) - (tc.y * ta.x);
 return((zab==0)?0:(zab<0)?((zbc<0)&&(zca<0)):((zbc>0)&&(zca>0)));
}

/*
 * We have a line, given by two points on it; we have a plane, given by
 *  an XYZW holding the X, Y, Z, and W values for the Xx+Yy+Zz=W form
 *  of the line.  Return the distance from the first of the two points
 *  to the point where the line intersects the plane.  Assume the line
 *  is not close enough to parallel to the plane to cause trouble.
 *
 * We compute u, a unit vector from p towards p2; then we write
 *
 * (p + A u) · XYZ = W
 *
 * (where XYZ is the <X,Y,Z> three-vector) and solve for A:
 *
 * p·XYZ + A u·XYZ = W
 * A = (W - p·XYZ) / (u·XYZ)
 *
 * This A is our return value, since |u| = 1.
 */
static double line2p_planexyzw_distance2(XYZ, XYZ, XYZW)__attribute__((__const__));
static double line2p_planexyzw_distance2(XYZ p, XYZ p2, XYZW plane)
{
 XYZ pxyz;

 pxyz = (XYZ){.x=plane.x,.y=plane.y,.z=plane.z};
 return(sqr((plane.w-dot(pxyz,p))/dot(pxyz,norm(xyzsub(p2,p)))));
}

static void sdp_(const char *, ...)
	__attribute__((__format__(__printf__,1,2)));
static void sdp_(const char *fmt, ...)
{
 static FILE *f = 0;
 va_list ap;

 if (! f)
  { f = fopen("sortdump","w");
    setlinebuf(f);
  }
 va_start(ap,fmt);
 vfprintf(f,fmt,ap);
 va_end(ap);
}
#define sdp(fmt, args...) \
	do { if (wantsort & SORT_DUMP) sdp_(fmt ,##args); } while (0)

static int occlude_line_line(LINE *la, LINE *lb)
{
 LIL lil;
 double da;
 double db;

 sdp("occlude_line_line: ");
 lil = line_intersect_line( la->a->scr[ex], la->b->scr[ex],
			    lb->a->scr[ex], lb->b->scr[ex] );
 if (! lil.inrange)
  { sdp("no\n");
    return(0);
  }
 /* Intersection.  Work out the threespace points corresponding to
    the intersection point and see which is closer to the eye. */
 da = xyzlen2(xyzsub(ea,xyzalpha(la->a->xyz,lil.A/lil.D,la->b->xyz)));
 db = xyzlen2(xyzsub(ea,xyzalpha(lb->a->xyz,lil.B/lil.D,lb->b->xyz)));
 sdp("da %g, db %g\n",da,db);
 return((da<db)?1:-1);
}

/*
 * To tell whether a line and face occlude one way or the other, we
 *  first test to see whether they overlap in screen coordinates.  This
 *  tells whether they occlude.  Then, to find out which way the
 *  occlusion goes, we find a point where both triangle and line are
 *  present and check which is in front.  (If the threespace line
 *  segment intersects the threespace triangle, we don't care which we
 *  claim is in front of the other.)
 */
static int occlude_line_face(LINE *a, FACE *b)
{
 XY la;
 XY lb;
 XY ta;
 XY tb;
 XY tc;
 LIL i;
 double alpha;
 XYZ lp;
 double dl;
 double dt;

 sdp("occlude_line_face: ");
 la = a->a->scr[ex];
 lb = a->b->scr[ex];
 ta = b->a->scr[ex];
 tb = b->b->scr[ex];
 tc = b->c->scr[ex];
 /* If either end of the line is inside the triangle, we have occlusion,
    and can use that line end as the depth-check point.  Otherwise, we
    have to check to see if the line cuts across the triangle. */
 do
  { if (pt_inside_triangle(la,ta,tb,tc))
     { alpha = 0;
       sdp("end 1 inside triangle: ");
       break;
     }
    if (pt_inside_triangle(lb,ta,tb,tc))
     { alpha = 1;
       sdp("end 2 inside triangle: ");
       break;
     }
    /* The line cuts across the triangle if it intersects any of the
       three line segments forming the triangle's sides. */
    do
     { i = line_intersect_line(la,lb,ta,tb);
       if (i.inrange)
	{ sdp("line cuts a-b: ");
	  break;
	}
       i = line_intersect_line(la,lb,tb,tc);
       if (i.inrange)
	{ sdp("line cuts b-c: ");
	  break;
	}
       i = line_intersect_line(la,lb,tc,ta);
       if (i.inrange)
	{ sdp("line cuts c-a: ");
	  break;
	}
       sdp("no overlap\n");
       return(0);
     } while (0);
    alpha = i.A / i.D;
    sdp("alpha=%g: ",alpha);
  } while (0);
 lp = xyzalpha(a->a->xyz,alpha,a->b->xyz);
 dl = xyzlen2(xyzsub(ea,lp));
 dt = line2p_planexyzw_distance2(ea,lp,b->plane);
 sdp("dl %g, dt %g\n",dl,dt);
 return((dl<dt)?1:-1);
}

/*
 * To see if two faces occlude one another, first we do a quick check
 *  to see if their on-screen bounding rectangles overlap; if not, we
 *  don't need to compute further.  Then we check if any edge segment
 *  of one intersects any edge segment of the other, if so, we have
 *  overlap.  If not, the only overlap possible is if one is wholly
 *  inside the other; to check for this, we check to see if points of
 *  one triangle are inside the other.  Unfortunately we have to check
 *  all three points, because otherwise two faces that share an edge
 *  and project such that one lies within the other can lose, 'cause
 *  the only point that works is the interior point on the smaller
 *  triangle.
 *
 * Even this has problems; if one triangle is wholly contained within
 *  the other on-screen, but all verticies of the smaller triangle are
 *  exactly on edges of the larger, it will hiccup.  This should be
 *  fixed, but for now ignore it; that's a set of measure zero.
 */
static int occlude_face_face(FACE *a, FACE *b)
{
 XY aa;
 XY ab;
 XY ac;
 XY ba;
 XY bb;
 XY bc;
 double a1;
 double a2;
 double b1;
 double b2;
 LIL i;
 double da;
 double db;

 sdp("occlude_face_face: ");
 aa = a->a->scr[ex];
 ab = a->b->scr[ex];
 ac = a->c->scr[ex];
 ba = b->a->scr[ex];
 bb = b->b->scr[ex];
 bc = b->c->scr[ex];
 a1 = aa.x;
 a2 = a1;
 if (ab.x < a1) a1 = ab.x;
 if (ab.x > a2) a2 = ab.x;
 if (ac.x < a1) a1 = ac.x;
 if (ac.x > a2) a2 = ac.x;
 b1 = ba.x;
 b2 = b1;
 if (bb.x < b1) b1 = bb.x;
 if (bb.x > b2) b2 = bb.x;
 if (bc.x < b1) b1 = bc.x;
 if (bc.x > b2) b2 = bc.x;
 if ((a2 < b1) || (a1 > b2))
  { sdp("no X overlap\n");
    return(0);
  }
 a1 = aa.y;
 a2 = a1;
 if (ab.y < a1) a1 = ab.y;
 if (ab.y > a2) a2 = ab.y;
 if (ac.y < a1) a1 = ac.y;
 if (ac.y > a2) a2 = ac.y;
 b1 = ba.y;
 b2 = b1;
 if (bb.y < b1) b1 = bb.y;
 if (bb.y > b2) b2 = bb.y;
 if (bc.y < b1) b1 = bc.y;
 if (bc.y > b2) b2 = bc.y;
 if ((a2 < b1) || (a1 > b2))
  { sdp("no Y overlap\n");
    return(0);
  }
 do <"occlude">
  {
#define FOO(E,F,G,H) do {\
	i = line_intersect_line(a##E,a##F,b##G,b##H);			\
	if (i.inrange)							\
	 { da = xyzlen2(xyzsub(xyzalpha(a->E->xyz,i.A/i.D,a->F->xyz),ea));\
	   db = xyzlen2(xyzsub(xyzalpha(b->G->xyz,i.B/i.D,b->H->xyz),ea));\
	   sdp("a "#E"-"#F" x b "#G"-"#H": ");				\
	   break <"occlude">;						\
	 } } while (0)
    FOO(a,b,a,b);
    FOO(a,b,b,c);
    FOO(a,b,c,a);
    FOO(b,c,a,b);
    FOO(b,c,b,c);
    FOO(b,c,c,a);
    FOO(c,a,a,b);
    FOO(c,a,b,c);
    FOO(c,a,c,a);
#undef FOO
#define FOO(t1,t2,p1) do {\
	if (pt_inside_triangle(t1##p1,t2##a,t2##b,t2##c))		\
	 { d##t1 = xyzlen2(xyzsub(t1->p1->xyz,ea));			\
	   d##t2 = line2p_planexyzw_distance2(ea,t1->p1->xyz,t2->plane);\
	   sdp(#t1#p1" inside "#t2": ");				\
	   break <"occlude">;						\
	 } } while(0)
    FOO(a,b,a);
    FOO(a,b,b);
    FOO(a,b,c);
    FOO(b,a,a);
    FOO(b,a,b);
    FOO(b,a,c);
#undef FOO
    sdp("no overlap\n");
    return(0);
  } while (0);
 sdp("da %g, db %g\n",da,db);
 return((da<db)?1:-1);
}

static int occlusion(OBJ *a, OBJ *b)
{
 switch (a->type)
  { default: abort(); break;
    case OT_LINE:
       switch (b->type)
	{ default: abort(); break;
	  case OT_LINE:
	     return(occlude_line_line(a->u.l,b->u.l));
	     break;
	  case OT_FACE:
	     return(occlude_line_face(a->u.l,b->u.f));
	     break;
	}
    case OT_FACE:
       switch (b->type)
	{ default: abort(); break;
	  case OT_LINE:
	     return(-occlude_line_face(b->u.l,a->u.f));
	     break;
	  case OT_FACE:
	     return(occlude_face_face(a->u.f,b->u.f));
	     break;
	}
  }
}

static void walk_down(OBJ *o, void (*fn)(OBJ *))
{
 static void step1(OBJ *o)
  { ARC *a;
    if (o->flags & OF_MARK) return;
    o->flags |= OF_MARK;
    (*fn)(o);
    for (a=o->out;a;a=a->flink) step1(a->t);
  }

 static void step2(OBJ *o)
  { ARC *a;
    if (o->flags & OF_MARK)
     { o->flags &= ~OF_MARK;
       for (a=o->out;a;a=a->flink) step2(a->t);
     }
  }

 step1(o);
 step2(o);
}

static void walk_up(OBJ *o, void (*fn)(OBJ *))
{
 static void step1(OBJ *o)
  { ARC *a;
    if (o->flags & OF_MARK) return;
    o->flags |= OF_MARK;
    (*fn)(o);
    for (a=o->in;a;a=a->tlink) step1(a->f);
  }

 static void step2(OBJ *o)
  { ARC *a;
    if (o->flags & OF_MARK)
     { o->flags &= ~OF_MARK;
       for (a=o->in;a;a=a->tlink) step2(a->f);
     }
  }

 step1(o);
 step2(o);
}

static void add_arc(OBJ *f, OBJ *t)
{
 ARC *a;

 static void foo1(OBJ *o)
  { occl[f->x][o->x] = 1;
    occl[o->x][f->x] = -1;
  }

 static void foo2(OBJ *o)
  { occl[o->x][t->x] = 1;
    occl[t->x][o->x] = -1;
  }

 a = malloc(sizeof(ARC));
 a->f = f;
 a->t = t;
 a->flink = f->out;
 f->out = a;
 a->tlink = t->in;
 t->in = a;
 walk_down(t,foo1);
 walk_up(f,foo2);
}

static void sort_objs(void)
{
 int i;
 int j;
 OBJLIST *root;
 OBJLIST **tail;

 memset(occl1d,0,nobjs*nobjs);
 for (i=0;i<nobjs;i++)
  { OBJ *o;
    o = objv[i];
    o->x = i;
    o->in = 0;
    o->out = 0;
  }
 for (i=0;i<nobjs;i++)
  { for (j=0;j<i;j++)
     { int o;
       if (occl[i][j]) continue;
       o = occlusion(objv[i],objv[j]);
       sdp("occlusion %s <-> %s = %d\n",objpstr(objv[i]),objpstr(objv[j]),o);
       if (o == 0) continue;
       add_arc(objv[(o<0)?j:i],objv[(o<0)?i:j]);
     }
  }
 /* We now have a forest of DAGs.  Find all the roots in this forest
    and collect them into a list "root".  See below for OF_MARK stuff. */
 tail = &root;
 for (i=0;i<nobjs;i++)
  { OBJ *o;
    OBJLIST *l;
    o = objv[i];
    o->x = -1;
    if (o->in) continue;
    l = malloc(sizeof(OBJLIST));
    l->o = o;
    o->flags |= OF_MARK;
    *tail = l;
    tail = &l->link;
  }
 *tail = 0;
 /* Now we want to assign indices from 0 to nobjs-1, with closer
    objects getting lower indices.  The way we do this amounts to a
    sort of cross between breadth-first and depth-first walk of the
    DAG we've just constructed.  We loop, pulling nodes off this "root"
    list.  For each node, if all its ins have had numbers assigned
    (including the vacuous case), we assign it the next number and push
    all its outs that aren't yet on the root list on the front.  (We
    use OF_MARK to indicate OBJs that are on the list.)  Otherwise, we
    move it to the end of the list.  Repeat until the list is empty. */
 j = 0;
 while <"assign"> (root)
  { OBJ *o;
    OBJLIST *l;
    ARC *a;
    l = root;
    root = l->link;
    if (! root) tail = &root;
    o = l->o;
    free(l);
    for (a=o->in;a;a=a->tlink)
     { if (a->f->x < 0)
	{ l = malloc(sizeof(OBJLIST));
	  l->o = o;
	  l->link = 0;
	  *tail = l;
	  tail = &l->link;
	  continue <"assign">;
	}
     }
    if (o->x >= 0) abort();
    sdp("assign %d to %s\n",j,objpstr(o));
    o->x = j;
    objv[j] = o;
    o->flags &= ~OF_MARK;
    j ++;
    for (a=o->out;a;a=a->flink)
     { if (a->t->flags & OF_MARK) continue;
       l = malloc(sizeof(OBJLIST));
       l->o = a->t;
       if (! root) tail = &l->link;
       l->link = root;
       root = l;
       a->t->flags |= OF_MARK;
     }
  }
 /* Sanity-check. */
 if (j != nobjs) abort();
 for (i=0;i<nobjs;i++)
  { OBJ *o;
    o = objv[i];
    if (o->flags & OF_MARK) abort();
    if (o->x != i) abort();
    while (o->out)
     { ARC *a;
       a = o->out;
       o->out = a->flink;
       free(a);
     }
  }
}

static void sdp_pt(PT *p)
{
 sdp("%g,%g,%g/%.0f,%.0f/%.0f,%.0f",
	p->xyz.x, p->xyz.y, p->xyz.z,
	p->scr[0].x*es, p->scr[0].y*es,
	p->scr[1].x*es, p->scr[1].y*es );
}

static void drawmodel(void)
{
 static int x;
 int i;

 fbline(0,0,x++&31,31,CMFIXED_WHITE);
 project_pts();
 if (wantsort)
  { if (wantsort & SORT_DUMP)
     { sdp("================ object list\n");
       for (i=nobjs-1;i>=0;i--)
	{ OBJ *o;
	  o = objs_l[i];
	  sdp("%s: ",objpstr(o));
	  switch (o->type)
	   { default: abort(); break;
	     case OT_LINE:
		sdp_pt(o->u.l->a);
		sdp("  ");
		sdp_pt(o->u.l->b);
		break;
	     case OT_FACE:
		sdp_pt(o->u.f->a);
		sdp("  ");
		sdp_pt(o->u.f->b);
		sdp("  ");
		sdp_pt(o->u.f->c);
		break;
	   }
	  sdp("\n");
	}
     }
    objv = objs_l;
    ex = 0;
    ea = lea.v;
    sdp("================ left eye view\n");
    sort_objs();
    objv = objs_r;
    ex = 1;
    ea = rea.v;
    sdp("================ right eye view\n");
    sort_objs();
    wantsort &= ~(SORT_ONCE|SORT_DUMP);
  }
 for (i=nobjs-1;i>=0;i--)
  { drawobj(objs_l[i],0,0);
    drawobj(objs_r[i],1,450);
  }
}

static void mayberender(void)
{
 if (! wantrender) return;
 wantrender = 0;
 fbdrain();
 fbclear();
 rt = norm(cross(up,out));
 rea.v = xyzadd(xyzscale(out,ed),xyzscale(rt,ipd/2));
 rea.u = norm(rea.v);
 rer = norm(cross(up,rea.v));
 lea.v = xyzadd(xyzscale(out,ed),xyzscale(rt,-ipd/2));
 lea.u = norm(lea.v);
 ler = norm(cross(up,lea.v));
 drawmodel();
 fbdrain();
 setcmap(drawbit);
 drawbit = ! drawbit;
 setpm(drawbit?0xf0:0x0f);
}

static int checkkbd(void)
{
 struct firm_event ev;
 int r;

 r = read(kbdfd,&ev,sizeof(ev));
 if (r == sizeof(ev))
  { queue_event(EVT_KBD|(ev.id&0x7f)|((ev.value==VKEY_UP)?EVK_UP:EVK_DN));
    return(1);
  }
 return(0);
}

static void set_wantsort(int bit)
{
 wantsort ^= bit;
 wantrender = 1;
}

static void process_evq(void)
{
 int ev;

 while (1)
  { ev = dequeue_event();
    switch (ev & EV_TYPE)
     { case EVT_NIL:
	  return;
	  break;
       case EVT_KBD:
	  switch (ev & (EVK_CODE|EVK_UPDN))
	   { default:
		printf("kbd: %d ",ev&EVK_CODE);
		switch (ev & EVK_UPDN)
		 { case EVK_UP: printf("up\n"); break;
		   case EVK_DN: printf("dn\n"); break;
		 }
		break;
/* L1 */     case   1|EVK_DN: mvbits |= MV_IN;  break;
	     case   1|EVK_UP: mvbits &=~MV_IN;  break;
/* L2 */     case   3|EVK_DN: mvbits |= MV_OUT; break;
	     case   3|EVK_UP: mvbits &=~MV_OUT; break;
/* L3 */     case  25|EVK_DN: mvbits |= MV_AP;  break;
	     case  25|EVK_UP: mvbits &=~MV_AP;  break;
/* L4 */     case  26|EVK_DN: mvbits |= MV_TG;  break;
	     case  26|EVK_UP: mvbits &=~MV_TG;  break;
/* L5 */     case  49|EVK_DN: mvbits |= MV_ZI;  break;
	     case  49|EVK_UP: mvbits &=~MV_ZI;  break;
/* L6 */     case  51|EVK_DN: mvbits |= MV_ZO;  break;
	     case  51|EVK_UP: mvbits &=~MV_ZO;  break;
/* L7 */     case  72|EVK_DN:                   break;
	     case  72|EVK_UP:                   break;
/* L8 */     case  73|EVK_DN:                   break;
	     case  73|EVK_UP:                   break;
/* L9 */     case  95|EVK_DN:                   break;
	     case  95|EVK_UP:                   break;
/* L10 */    case  97|EVK_DN:                   break;
	     case  97|EVK_UP:                   break;
/* F1 */     case   5|EVK_DN: set_wantsort(1);  break;
	     case   5|EVK_UP:                   break;
/* F2 */     case   6|EVK_DN: set_wantsort(2);  break;
	     case   6|EVK_UP:                   break;
/* F3 */     case   8|EVK_DN: set_wantsort(4);  break;
	     case   8|EVK_UP:                   break;
/* F4 */     case  10|EVK_DN:                   break;
	     case  10|EVK_UP:                   break;
/* F5 */     case  12|EVK_DN:                   break;
	     case  12|EVK_UP:                   break;
/* F6 */     case  14|EVK_DN:                   break;
	     case  14|EVK_UP:                   break;
/* F7 */     case  16|EVK_DN: resethw();        break;
	     case  16|EVK_UP:                   break;
/* F8 */     case  17|EVK_DN:                   break;
	     case  17|EVK_UP:                   break;
/* F9 */     case  18|EVK_DN:                   break;
	     case  18|EVK_UP:                   break;
/* ESC */    case  29|EVK_DN:                   break;
	     case  29|EVK_UP: exit(0);          break;
/* BS */     case  43|EVK_DN:                   break;
	     case  43|EVK_UP:                   break;
/* R7 */     case  68|EVK_DN: mvbits |= MV_CCW; break;
	     case  68|EVK_UP: mvbits &=~MV_CCW; break;
/* R8 */     case  69|EVK_DN: mvbits |= MV_UP;  break;
	     case  69|EVK_UP: mvbits &=~MV_UP;  break;
/* R9 */     case  70|EVK_DN: mvbits |= MV_CW;  break;
	     case  70|EVK_UP: mvbits &=~MV_CW;  break;
/* R10 */    case  91|EVK_DN: mvbits |= MV_LF;  break;
	     case  91|EVK_UP: mvbits &=~MV_LF;  break;
/* R11 */    case  92|EVK_DN:                   break;
	     case  92|EVK_UP:                   break;
/* R12 */    case  93|EVK_DN: mvbits |= MV_RT;  break;
	     case  93|EVK_UP: mvbits &=~MV_RT;  break;
/* R13 */    case 112|EVK_DN:                   break;
	     case 112|EVK_UP:                   break;
/* R14 */    case 113|EVK_DN: mvbits |= MV_DN;  break;
	     case 113|EVK_UP: mvbits &=~MV_DN;  break;
/* R15 */    case 114|EVK_DN:                   break;
	     case 114|EVK_UP:                   break;
	   }
	  break;
       case EVT_TICK:
	  switch (mvbits & (MV_LF|MV_RT))
	   { case MV_RT:
		out = rotunit(out,up,.01);
		wantrender = 1;
		break;
	     case MV_LF:
		out = rotunit(out,up,-.01);
		wantrender = 1;
		break;
	   }
	  switch (mvbits & (MV_UP|MV_DN))
	   { case MV_DN:
		rt = norm(cross(up,out));
		up = rotunit(up,rt,.01);
		out = norm(cross(rt,up));
		wantrender = 1;
		break;
	     case MV_UP:
		rt = norm(cross(up,out));
		up = rotunit(up,rt,-.01);
		out = norm(cross(rt,up));
		wantrender = 1;
		break;
	   }
	  switch (mvbits & (MV_CW|MV_CCW))
	   { case MV_CW:
		rt = norm(cross(up,out));
		up = rotunit(up,out,.01);
		wantrender = 1;
		break;
	     case MV_CCW:
		rt = norm(cross(up,out));
		up = rotunit(up,out,-.01);
		wantrender = 1;
		break;
	   }
	  switch (mvbits & (MV_IN|MV_OUT))
	   { case MV_IN:
		ed -= .02;
		wantrender = 1;
		break;
	     case MV_OUT:
		ed += .02;
		wantrender = 1;
		break;
	   }
	  switch (mvbits & (MV_AP|MV_TG))
	   { case MV_AP:
		ipd += .002;
		wantrender = 1;
		break;
	     case MV_TG:
		ipd -= .002;
		wantrender = 1;
		break;
	   }
	  switch (mvbits & (MV_ZI|MV_ZO))
	   { case MV_ZI:
		es += 1;
		wantrender = 1;
		break;
	     case MV_ZO:
		es -= 1;
		wantrender = 1;
		break;
	   }
	  break;
     }
  }
}

static void initsim(void)
{
 up.x = 0;
 up.y = 0;
 up.z = -1;
 out.x = 1;
 out.y = .2;
 out.z = .3;
 out = norm(out);
 rt = norm(cross(up,out));
 up = norm(cross(out,rt));
 ed = 10;
 es = 1000;
 ipd = 1;
 mvbits = 0;
}

static void run(void)
{
 struct pollfd pfd;

 process_evq();
 mayberender();
 if (checkkbd()) return;
 if (mvbits) starttiming(); else stoptiming();
 if (checktick()) return;
 pfd.fd = kbdfd;
 pfd.events = POLLIN | POLLRDNORM;
 if (poll(&pfd,1,INFTIM) < 0)
  { if (errno == EINTR) return;
    fprintf(stderr,"%s: poll: %s\n",__progname,strerror(errno));
    exit(1);
  }
}

static int setdt(const char *s)
{
 if (! strcmp(s,"none"))
  { disptype = DT_NONE;
  }
 else if (! strcmp(s,"normal"))
  { disptype = DT_NORMAL;
  }
 else
  { fprintf(stderr,"%s: bad display type `%s'\n",__progname,s);
    return(1);
  }
 return(0);
}

static void usage(void)
{
 fprintf(stderr,"Usage: %s [-fb fbdev] [-kb|-kbd kbdev] modelfile\n",__progname);
}

static void handleargs(int ac, char **av)
{
 int skip;
 int errs;

 skip = 0;
 errs = 0;
 for (ac--,av++;ac;ac--,av++)
  { if (skip > 0)
     { skip --;
       continue;
     }
    if (**av != '-')
     { if (! model)
	{ model = *av;
	}
       else
	{ fprintf(stderr,"%s: extra argument `%s'\n",__progname,*av);
	  errs ++;
	}
       continue;
     }
    if (0)
     {
needarg:;
       fprintf(stderr,"%s: %s needs a following argument\n",__progname,*av);
       errs ++;
       continue;
     }
#define WANTARG() do { if (++skip >= ac) goto needarg; } while (0)
    if (!strcmp(*av,"-fb"))
     { WANTARG();
       fbpath = av[skip];
       continue;
     }
    if (!strcmp(*av,"-kb") || !strcmp(*av,"-kbd"))
     { WANTARG();
       kbpath = av[skip];
       continue;
     }
    if (!strcmp(*av,"-dt") || !strcmp(*av,"-displaytype"))
     { WANTARG();
       errs += setdt(av[skip]);
       continue;
     }
#undef WANTARG
    fprintf(stderr,"%s: unrecognized option `%s'\n",__progname,*av);
    errs ++;
  }
 if (! model)
  { usage();
    errs ++;
  }
 if (errs) exit(1);
}

static char *readaline(FILE *f)
{
 char *b;
 int a;
 int l;
 int c;

 static void add(char ch)
  { if (l >= a) b = realloc(b,a=l+16);
    b[l++] = ch;
  }

 b = 0;
 a = 0;
 l = 0;
 while (1)
  { c = getc(f);
    if (c == EOF)
     { if (l > 0) add('\0');
       return(b);
     }
    if (c == '\n')
     { add('\0');
       return(b);
     }
    add(c);
  }
}

static const char *skipws(const char *s)
{
 while (*s && isspace(*s)) s ++;
 return(s);
}

static const char *skipnws(const char *s)
{
 while (*s && !isspace(*s)) s ++;
 return(s);
}

static void badline(const char *, ...)
	__attribute__((__format__(__printf__,1,2),__noreturn__));
static void badline(const char *fmt, ...)
{
 va_list ap;

 va_start(ap,fmt);
 (*linefail)(fmt,&ap);
 va_end(ap);
 abort();
}

static void model_line_pt(const char *lp)
{
 const char *ep;
 char *newname;
 int l;
 PT *pt;
 double x;
 double y;
 double z;

 ep = skipnws(lp);
 l = ep - lp;
 for (pt=pts;pt;pt=pt->link)
  { if ((l == pt->namelen) && !strncmp(lp,pt->name,l))
     { badline("redefinition of point `%s'",pt->name);
     }
  }
 newname = malloc(l+1);
 bcopy(lp,newname,l);
 newname[l] = '\0';
 if (sscanf(ep,"%lg%lg%lg",&x,&y,&z) != 3)
  { free(newname);
    badline("can't scan coordinates");
  }
 pt = malloc(sizeof(PT));
 pt->name = newname;
 pt->namelen = l;
 pt->xyz.x = x;
 pt->xyz.y = y;
 pt->xyz.z = z;
 pt->link = pts;
 pts = pt;
}

static const char *getpoint(const char *lp, PT **pp, int musthave)
{
 const char *ep;

 lp = skipws(lp);
 if (*lp == '(')
  { double x;
    double y;
    double z;
    PT *pt;
    lp ++;
    ep = lp;
    while (*ep && (*ep != ')')) ep ++;
    if (! *ep) badline("missing ) for anonymous point");
    if (sscanf(lp,"%lg%lg%lg",&x,&y,&z) != 3) badline("can't scan anonymous point coords");
    pt = malloc(sizeof(PT));
    pt->name = 0;
    pt->namelen = -1;
    pt->xyz.x = x;
    pt->xyz.y = y;
    pt->xyz.z = z;
    pt->link = pts;
    pts = pt;
    *pp = pt;
    return(ep+1);
  }
 else if (*lp == '\0')
  { if (musthave) badline("missing point");
    *pp = 0;
    return(lp);
  }
 else
  { int l;
    PT *pt;
    ep = skipnws(lp);
    l = ep - lp;
    for (pt=pts;pt;pt=pt->link)
     { if ((l == pt->namelen) && !strncmp(lp,pt->name,l)) break;
     }
    if (! pt) badline("reference to undefined point `%.*s'",l,lp);
    *pp = pt;
    return(ep);
  }
}

static const char *getcolor(const char *lp, COLOR **cp)
{
 const char *ep;
 int l;
 COLOR *c;

 lp = skipws(lp);
 if (*lp == '\0') badline("missing color name");
 ep = skipnws(lp);
 l = ep - lp;
 for (c=colors;c;c=c->link)
  { if ((l == c->namelen) && !strncmp(lp,c->name,l)) break;
  }
 if (! c) badline("reference to undefined color `%.*s'",l,lp);
 if (! (c->flags & CF_USED))
  { c->flags |= CF_USED;
    c->inx = nextcx++;
  }
 *cp = c;
 return(ep);
}

static void line_precompute(LINE *l __attribute__((__unused__)))
{
}

static void model_line_line(const char *lp)
{
 const char *ep;
 static char *newname = 0;
 int l;
 LINE *ln;
 PT *pt1;
 PT *pt2;
 COLOR *c;

 if (newname) free(newname);
 ep = skipnws(lp);
 l = ep - lp;
 if ((l == 1) && (*lp == '-'))
  { newname = 0;
    l = -1;
  }
 else
  { for (ln=lines;ln;ln=ln->link)
     { if ((l == ln->namelen) && !strncmp(lp,ln->name,l))
	{ badline("redefinition of line `%s'",ln->name);
	}
     }
    newname = malloc(l+1);
    bcopy(lp,newname,l);
    newname[l] = '\0';
  }
 ep = getcolor(ep,&c);
 ep = getpoint(ep,&pt1,1);
 ep = getpoint(ep,&pt2,1);
 ln = malloc(sizeof(LINE));
 ln->name = newname;
 ln->namelen = l;
 ln->a = pt1;
 ln->b = pt2;
 ln->color = c;
 ln->link = lines;
 lines = ln;
 line_precompute(ln);
 newname = 0;
}

/*
 * To solve for the plane, we find its normal by computing (p2-p1)×(p3-p1).
 *  This gives us A, B, and C in the Ax+By+Cz+D=0 form of the plane.
 *  Then we get D by dotting the normal vector with p1.
 *
 * We check for collinearity by checking to see if the cross product
 *  vector is very short.  If so, we fail.
 */
static int face_precompute(FACE *f)
{
 XYZ cprod;
 double cplen;

 cprod = cross( norm(xyzsub(f->b->xyz,f->a->xyz)),
		norm(xyzsub(f->c->xyz,f->a->xyz)) );
 cplen = xyzlen(cprod);
 if (cplen < 1e-5) return(0);
 cprod = xyzscale(cprod,1/cplen);
 f->plane.x = cprod.x;
 f->plane.y = cprod.y;
 f->plane.z = cprod.z;
 f->plane.w = dot(cprod,f->a->xyz);
 return(1);
}

static void model_line_face(const char *lp, int autotri)
{
 const char *ep;
 static char *newname = 0;
 static int verta = 0;
 static PT **verts = 0;
 int vertn;
 int l;
 FACE *f;
 PT *pt;
 COLOR *c;

 static void add_face(FACE *f, char *name)
  { if (face_precompute(f))
     { f->color = c;
       f->link = faces;
       f->name = name;
       f->namelen = strlen(name);
       faces = f;
     }
    else
     { badline("points too close to collinear");
     }
  }

 if (newname) free(newname);
 ep = skipnws(lp);
 l = ep - lp;
 if ((l == 1) && (*lp == '-'))
  { newname = 0;
    l = -1;
  }
 else
  {
#if 0
    for (fc=faces;fc;fc=fc->link)
     { if ((l == fc->namelen) && !strncmp(lp,fc->name,l))
	{ badline("redefinition of face `%s'",fc->name);
	}
     }
#endif
    newname = malloc(l+1);
    bcopy(lp,newname,l);
    newname[l] = '\0';
  }
 ep = getcolor(ep,&c);
 vertn = 0;
 while (1)
  { lp = ep;
    ep = getpoint(lp,&pt,0);
    if (! pt)
     { if (vertn < 3)
	{ badline("too few vertices");
	}
       else if (vertn > 3)
	{ int inx;
	  char *n;
	  if (! autotri) badline("too many vertices");
	  inx = 1;
	  while (vertn > 3)
	   { vertn --;
	     f = malloc(sizeof(FACE));
	     f->a = verts[0];
	     f->b = verts[1];
	     f->c = verts[vertn];
	     asprintf(&n,"%s (%d)",newname,inx++);
	     add_face(f,n);
	     verts[0] = verts[vertn];
	   }
	  f = malloc(sizeof(FACE));
	  f->a = verts[0];
	  f->b = verts[1];
	  f->c = verts[2];
	  asprintf(&n,"%s (%d)",newname,inx);
	  add_face(f,n);
	}
       else
	{ f = malloc(sizeof(FACE));
	  f->a = verts[0];
	  f->b = verts[1];
	  f->c = verts[2];
	  add_face(f,newname);
	  newname = 0;
	}
       break;
     }
    if (vertn >= verta) verts = realloc(verts,(verta=vertn+4)*sizeof(PT *));
    verts[vertn++] = pt;
  }
 free(newname);
 newname = 0;
}

static void model_line_color(const char *lp)
{
 const char *ep;
 char *newname;
 int l;
 COLOR *c;
 int r;
 int g;
 int b;

 ep = skipnws(lp);
 l = ep - lp;
 if ((l == 1) && (*lp == '-'))
  { newname = 0;
    l = -1;
  }
 else
  { for (c=colors;c;c=c->link)
     { if ((l == c->namelen) && !strncmp(lp,c->name,l))
	{ badline("redefinition of color `%s'",c->name);
	}
     }
    newname = malloc(l+1);
    bcopy(lp,newname,l);
    newname[l] = '\0';
  }
 if (sscanf(ep,"%d%d%d",&r,&g,&b) != 3)
  { free(newname);
    badline("can't scan rgb triple");
  }
 if ((r < 0) || (r > 255) || (g < 0) || (g > 255) || (b < 0) || (b > 255))
  { free(newname);
    badline("rgb value out of range 0..255");
  }
 c = malloc(sizeof(COLOR));
 c->name = newname;
 c->namelen = l;
 c->r = r;
 c->g = g;
 c->b = b;
 c->flags = 0;
 c->link = colors;
 colors = c;
}

static int model_line(const char *line)
{
 __label__ lret;

 const char *lp;
 const char *ep;
 int l;

 static void line_fail(const char *, va_list *)
		__attribute__((__format__(__printf__,1,0),__noreturn__));
 static void line_fail(const char *fmt, va_list *app)
  { fprintf(stderr,"%s: bad model line [",__progname);
    vfprintf(stderr,fmt,*app);
    fprintf(stderr,"]: %s\n",line);
    goto lret;
  }

 if (0)
  {
lret:;
    return(1);
  }
 linefail = &line_fail;
 lp = skipws(line);
 if (*lp == '#') return(0);
 if (*lp == '\0') return(0);
 ep = skipnws(lp);
 l = ep - lp;
 ep = skipws(ep);
 if ((l == 2) && !strncmp(lp,"pt",2))
  { model_line_pt(ep);
  }
 else if ((l == 4) && !strncmp(lp,"line",4))
  { model_line_line(ep);
  }
 else if ((l == 4) && !strncmp(lp,"face",4))
  { model_line_face(ep,0);
  }
 else if ((l == 5) && !strncmp(lp,"face*",5))
  { model_line_face(ep,1);
  }
 else if ((l == 5) && !strncmp(lp,"color",5))
  { model_line_color(ep);
  }
 else
  { badline("unrecognized keyword `%.*s'",l,lp);
  }
 return(0);
}

static void initmodel(void)
{
 FILE *f;
 char *line;
 int errs;
 LINE *ln;
 FACE *fc;
 int i;

 f = fopen(model,"r");
 if (f == 0)
  { fprintf(stderr,"%s: can't read %s: %s\n",__progname,model,strerror(errno));
    exit(1);
  }
 pts = 0;
 lines = 0;
 faces = 0;
 errs = 0;
 while ((line=readaline(f)))
  { errs += model_line(line);
    free(line);
  }
 fclose(f);
 if (errs > 0) exit(1);
 nobjs = 0;
 for (ln=lines;ln;ln=ln->link) nobjs ++;
 for (fc=faces;fc;fc=fc->link) nobjs ++;
 objs_l = malloc(2*nobjs*sizeof(OBJ *));
 objs_r = objs_l + nobjs;
 objs_l[0] = malloc(nobjs*sizeof(OBJ));
 for (i=1;i<nobjs;i++) objs_l[i] = objs_l[0]++;
 i = 0;
 for (ln=lines;ln;ln=ln->link)
  { objs_l[i]->type = OT_LINE;
    objs_l[i]->flags = 0;
    objs_l[i]->u.l = ln;
    i ++;
  }
 for (fc=faces;fc;fc=fc->link)
  { objs_l[i]->type = OT_FACE;
    objs_l[i]->flags = 0;
    objs_l[i]->u.f = fc;
    i ++;
  }
 bcopy(&objs_l[0],&objs_r[0],nobjs*sizeof(OBJ *));
 occl = malloc(nobjs*sizeof(char *));
 occl1d = malloc(nobjs*nobjs);
 occl[0] = occl1d;
 for (i=1;i<nobjs;i++) occl[i] = occl[i-1] + nobjs;
}

static void initcmap(void)
{
 int i;
 COLOR *c;

 if (nextcx > 16)
  { fprintf(stderr,"%s: too many distinct colors\n",__progname);
    exit(1);
  }
 for (i=0;i<16;i++)
  { cmap[i][0] = 0;
    cmap[i][1] = 0;
    cmap[i][2] = 0;
  }
 cmap[CMFIXED_BLACK][0] = 0;
 cmap[CMFIXED_BLACK][1] = 0;
 cmap[CMFIXED_BLACK][2] = 0;
 cmap[CMFIXED_WHITE][0] = 255;
 cmap[CMFIXED_WHITE][1] = 255;
 cmap[CMFIXED_WHITE][2] = 255;
 for (c=colors;c;c=c->link)
  { if (c->flags & CF_USED)
     { cmap[c->inx][0] = c->r;
       cmap[c->inx][1] = c->g;
       cmap[c->inx][2] = c->b;
     }
  }
 printf("colormap:\n");
 for (i=0;i<16;i++)
  { printf("%2d: %3d %3d %3d\n",i,cmap[i][0],cmap[i][1],cmap[i][2]);
  }
}

int main(int, char **);
int main(int ac, char **av)
{
 handleargs(ac,av);
 initmodel();
 initsim();
 initevq();
 initfb();
 initkbd();
 initcmap();
 inittiming();
 wantrender = 1;
 while (1) run();
}