Imported from /srv/lorry/lorry-area/swig-tarball/swig-1.3.40.tar.gz.HEADswig-1.3.40master

9 files changed, 5594 insertions, 0 deletions

diff --git a/Examples/GIFPlot/Lib/Makefile.in b/Examples/GIFPlot/Lib/Makefile.in
new file mode 100644
index 0000000..9db828e
--- /dev/null
+++ b/Examples/GIFPlot/Lib/Makefile.in
@@ -0,0 +1,22 @@
+CC      = @CC@
+CCSHARED= @CCSHARED@ 
+INCLUDES= -I../Include
+CFLAGS  = -O
+SRCS    = frame.c color.c plot2d.c plot3d.c font.c pixmap.c matrix.c gif.c
+OBJS    = $(SRCS:.c=.@OBJEXT@)
+AR      = @AR@
+RANLIB  = @RANLIB@
+TARGET  = ../libgifplot.a
+
+.c.@OBJEXT@: 
+	$(CC) $(CCSHARED) $(INCLUDES) $(CFLAGS) -c -o $*.@OBJEXT@ $<
+
+all:  $(OBJS)
+	@rm -f ../libgifplot.a
+	$(AR) cr $(TARGET) $(OBJS)
+	$(RANLIB) $(TARGET)
+
+clean:
+	rm -f *.@OBJEXT@ *~ $(TARGET)
+
+check: all
diff --git a/Examples/GIFPlot/Lib/color.c b/Examples/GIFPlot/Lib/color.c
new file mode 100644
index 0000000..7d79bac
--- /dev/null
+++ b/Examples/GIFPlot/Lib/color.c
@@ -0,0 +1,143 @@
+/* ----------------------------------------------------------------------------- 
+ * color.c
+ *
+ *     Colormaps
+ * 
+ * Author(s) : David Beazley (beazley@cs.uchicago.edu)
+ * Copyright (C) 1995-1996
+ *
+ * See the file LICENSE for information on usage and redistribution.	
+ * ----------------------------------------------------------------------------- */
+
+#define COLORMAP
+#include "gifplot.h"
+#include <string.h>
+
+/* ------------------------------------------------------------------------
+   ColorMap *new_ColorMap(char *filename)
+
+   Read a colormap from a file.
+   ------------------------------------------------------------------------ */
+
+ColorMap *new_ColorMap(char *filename) {
+  ColorMap *c;
+  FILE *cm;
+  void ColorMap_default(ColorMap *);
+  
+  if (!filename) {
+    c = (ColorMap *) malloc(sizeof(ColorMap));
+    c->cmap = (unsigned char *) malloc(768*sizeof(char));
+    c->name = 0;
+    ColorMap_default(c);
+    return c;
+  }
+  if (strlen(filename) == 0) {
+    c = (ColorMap *) malloc(sizeof(ColorMap));
+    c->cmap = (unsigned char *) malloc(768*sizeof(char));
+    ColorMap_default(c);
+    return c;
+  }
+  if ((cm = fopen(filename,"rb")) == NULL) {
+    return (ColorMap *) 0;
+  }
+  c = (ColorMap *) malloc(sizeof(ColorMap));
+  c->cmap = (unsigned char *) malloc(768*sizeof(char));
+  if (fread(c->cmap, 768, 1, cm) != 1) {
+    free((char *)c->cmap);
+    free((char *)c);
+    fclose(cm);
+    return (ColorMap *) 0;
+  }
+  fclose(cm);
+  c->name = (char *) malloc(strlen(filename)+1);
+  strcpy(c->name, filename);
+  ColorMap_default(c);
+  return c;
+}
+
+/* --------------------------------------------------------------------------
+   delete_ColorMap(ColorMap *cm)
+
+   Destroy a ColorMap
+   -------------------------------------------------------------------------- */
+
+void delete_ColorMap(ColorMap *cm) {
+  if (cm) {
+    free((char *)cm->cmap);
+    if (cm->name)
+      free((char *)cm->name);
+    free((char *)cm);
+  }
+}
+
+/* --------------------------------------------------------------------------
+   ColorMap_default(ColorMap *cm)
+
+   This function fills in default values for the first 8 entries of the
+   colormap.  These are *fixed* values---used internally.
+   -------------------------------------------------------------------------- */
+
+void ColorMap_default(ColorMap *cm) {
+  unsigned char *r,*g,*b;
+  if (cm) {
+    r = cm->cmap;
+    g = cm->cmap+256;
+    b = cm->cmap+512;
+
+    r[0] = 0;   g[0] = 0;   b[0] = 0;      /* BLACK */
+    r[1] = 255; g[1] = 255; b[1] = 255;    /* WHITE */
+    r[2] = 255; g[2] = 0;   b[2] = 0;      /* RED   */
+    r[3] = 0;   g[3] = 255; b[3] = 0;      /* GREEN */
+    r[4] = 0;   g[4] = 0;   b[4] = 255;    /* BLUE  */
+    r[5] = 255; g[5] = 255; b[5] = 0;      /* YELLOW */
+    r[6] = 0;   g[6] = 255; b[6] = 255;    /* CYAN   */
+    r[7] = 255; g[7] = 0;   b[7] = 255;    /* MAGENTA */
+  }
+}
+
+void ColorMap_assign(ColorMap *cm, int index, int r, int g, int b) {
+  unsigned char *rb,*gb,*bb;
+
+  rb = cm->cmap;
+  gb = cm->cmap+256;
+  bb = cm->cmap+512;
+
+  rb[index] = r;
+  gb[index] = g;
+  bb[index] = b;
+}
+
+int ColorMap_getitem(ColorMap *cm, int index) {
+  if ((index < 0) && (index >= 768)) return -1;
+  return cm->cmap[index];
+}
+
+void ColorMap_setitem(ColorMap *cm, int index, int value) {
+  if ((index < 0) && (index >= 768)) return;
+  cm->cmap[index]=value;
+}  
+
+/* --------------------------------------------------------------------
+   ColorMap_write(ColorMap *cm, char *filename)
+
+   Write out a colormap to disk.
+   -------------------------------------------------------------------- */
+
+int ColorMap_write(ColorMap *cm, char *filename) {
+
+  FILE *f;
+  if (!cm) return -1;
+  if (!filename) return -1;
+  if (strlen(filename) == 0) return -1;
+
+  f = fopen(filename,"w");
+  
+  if (fwrite(cm->cmap,768,1,f) != 1) {
+    fclose(f);
+    return -1;
+  }
+  fclose(f);
+  return 0;
+}
+
+#undef COLORMAP
diff --git a/Examples/GIFPlot/Lib/font.c b/Examples/GIFPlot/Lib/font.c
new file mode 100644
index 0000000..b30ee9b
--- /dev/null
+++ b/Examples/GIFPlot/Lib/font.c
@@ -0,0 +1,705 @@
+/* ----------------------------------------------------------------------------- 
+ * font.c
+ *
+ *     Some basic fonts.
+ * 
+ * Author(s) : David Beazley (beazley@cs.uchicago.edu)
+ * Copyright (C) 1995-1996
+ *
+ * See the file LICENSE for information on usage and redistribution.	
+ * ----------------------------------------------------------------------------- */
+
+#define FONT
+#include "gifplot.h"
+
+static char Char_A[80] = "\
+...x....\
+...x....\
+..x.x...\
+..x.x...\
+.x...x..\
+.xxxxx..\
+x.....x.\
+x.....x.\
+x.....x.\
+........";
+
+static char Char_B[80] = "\
+xxxxxx..\
+x.....x.\
+x.....x.\
+x.....x.\
+xxxxxx..\
+x.....x.\
+x.....x.\
+x.....x.\
+xxxxxx..\
+........";
+
+static char Char_C[80] = "\
+..xxxx..\
+.x....x.\
+x.......\
+x.......\
+x.......\
+x.......\
+x.......\
+.x....x.\
+..xxxx..\
+........";
+
+static char Char_D[80] = "\
+xxxxx...\
+x....x..\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+x....x..\
+xxxxx...\
+........";
+static char Char_E[80] = "\
+xxxxxxx.\
+x.......\
+x.......\
+x.......\
+xxxxx...\
+x.......\
+x.......\
+x.......\
+xxxxxxx.\
+........";
+static char Char_F[80] = "\
+xxxxxxx.\
+x.......\
+x.......\
+x.......\
+xxxxx...\
+x.......\
+x.......\
+x.......\
+x.......\
+........";
+static char Char_G[80] = "\
+.xxxxx..\
+x.....x.\
+x.......\
+x.......\
+x...xxx.\
+x.....x.\
+x.....x.\
+x.....x.\
+.xxxxx..\
+........";
+static char Char_H[80] = "\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+xxxxxxx.\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+........";
+static char Char_I[80] = "\
+xxxxxxx.\
+...x....\
+...x....\
+...x....\
+...x....\
+...x....\
+...x....\
+...x....\
+xxxxxxx.\
+........";
+static char Char_J[80] = "\
+......x.\
+......x.\
+......x.\
+......x.\
+......x.\
+......x.\
+x.....x.\
+.x...x..\
+..xxx...\
+........";
+static char Char_K[80] = "\
+x.....x.\
+x....x..\
+x...x...\
+x..x....\
+xxx.....\
+x..x....\
+x...x...\
+x....x..\
+x.....x.\
+........";
+static char Char_L[80] = "\
+x.......\
+x.......\
+x.......\
+x.......\
+x.......\
+x.......\
+x.......\
+x.......\
+xxxxxxx.\
+........";
+static char Char_M[80] = "\
+x.....x.\
+xx...xx.\
+xx...xx.\
+x.x.x.x.\
+x.x.x.x.\
+x..x..x.\
+x..x..x.\
+x.....x.\
+x.....x.\
+........";
+static char Char_N[80] = "\
+x.....x.\
+xx....x.\
+x.x...x.\
+x.x...x.\
+x..x..x.\
+x...x.x.\
+x...x.x.\
+x....xx.\
+x.....x.\
+........";
+static char Char_O[80] = "\
+.xxxxx..\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+.xxxxx..\
+........";
+static char Char_P[80] = "\
+xxxxxx..\
+x.....x.\
+x.....x.\
+x.....x.\
+xxxxxx..\
+x.......\
+x.......\
+x.......\
+x.......\
+........";
+static char Char_Q[80] = "\
+.xxxxx..\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+x...x.x.\
+x....x..\
+.xxxx.x.\
+........";
+static char Char_R[80] = "\
+xxxxxx..\
+x.....x.\
+x.....x.\
+x.....x.\
+xxxxxx..\
+x..x....\
+x...x...\
+x....x..\
+x.....x.\
+........";
+static char Char_S[80] = "\
+.xxxxx..\
+x.....x.\
+x.......\
+x.......\
+.xxxxx..\
+......x.\
+......x.\
+x.....x.\
+.xxxxx..\
+........";
+static char Char_T[80] = "\
+xxxxxxx.\
+...x....\
+...x....\
+...x....\
+...x....\
+...x....\
+...x....\
+...x....\
+...x....\
+........";
+static char Char_U[80] = "\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+.xxxxx..\
+........";
+static char Char_V[80] = "\
+x.....x.\
+x.....x.\
+.x...x..\
+.x...x..\
+..x.x...\
+..x.x...\
+...x....\
+...x....\
+...x....\
+........";
+static char Char_W[80] = "\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+x.....x.\
+x..x..x.\
+x..x..x.\
+x.x.x.x.\
+.x...x..\
+........";
+static char Char_X[80] = "\
+x.....x.\
+x.....x.\
+.x...x..\
+..x.x...\
+...x....\
+..x.x...\
+.x...x..\
+x.....x.\
+x.....x.\
+........";
+static char Char_Y[80] = "\
+x.....x.\
+x.....x.\
+.x...x..\
+..x.x...\
+...x....\
+...x....\
+...x....\
+...x....\
+...x....\
+........";
+static char Char_Z[80] = "\
+xxxxxxx.\
+......x.\
+.....x..\
+....x...\
+...x....\
+..x.....\
+.x......\
+x.......\
+xxxxxxx.\
+........";
+static char Char_0[80] = "\
+.xxxxx..\
+x....xx.\
+x...x.x.\
+x..x..x.\
+x..x..x.\
+x.x...x.\
+x.x...x.\
+xx....x.\
+.xxxxx..\
+........";
+static char Char_1[80] = "\
+...x....\
+..xx....\
+...x....\
+...x....\
+...x....\
+...x....\
+...x....\
+...x....\
+..xxx...\
+........";
+static char Char_2[80] = "\
+..xxxx..\
+.x....x.\
+x.....x.\
+.....x..\
+....x...\
+...x....\
+..x.....\
+.x......\
+xxxxxxx.\
+........";
+static char Char_3[80] = "\
+.xxxxx..\
+x.....x.\
+......x.\
+......x.\
+...xxx..\
+......x.\
+......x.\
+x.....x.\
+.xxxxx..\
+........";
+static char Char_4[80] = "\
+....xx..\
+...x.x..\
+..x..x..\
+.x...x..\
+xxxxxxx.\
+.....x..\
+.....x..\
+.....x..\
+.....x..\
+........";
+static char Char_5[80] = "\
+xxxxxxx.\
+x.......\
+x.......\
+x.......\
+xxxxxx..\
+......x.\
+......x.\
+x.....x.\
+.xxxxx..\
+........";
+static char Char_6[80] = "\
+....xxx.\
+..xx....\
+.x......\
+x.......\
+x.xxx...\
+xx...x..\
+x.....x.\
+.x...x..\
+..xxx...\
+........";
+static char Char_7[80] = "\
+xxxxxxx.\
+x.....x.\
+.....x..\
+....x...\
+...x....\
+..x.....\
+.x......\
+x.......\
+x.......\
+........";
+static char Char_8[80] = "\
+.xxxxx..\
+x.....x.\
+x.....x.\
+x.....x.\
+.xxxxx..\
+x.....x.\
+x.....x.\
+x.....x.\
+.xxxxx..\
+........";
+static char Char_9[80] = "\
+..xxxx..\
+.x....x.\
+x.....x.\
+x....xx.\
+.xxxx.x.\
+......x.\
+......x.\
+....xx..\
+.xxx....\
+........";
+static char Char_MINUS[80] = "\
+........\
+........\
+........\
+........\
+.xxxxxx.\
+........\
+........\
+........\
+........\
+........";
+static char Char_PLUS[80] = "\
+........\
+........\
+...x....\
+...x....\
+.xxxxx..\
+...x....\
+...x....\
+........\
+........\
+........";
+static char Char_EQUAL[80] = "\
+........\
+........\
+........\
+.xxxxx..\
+........\
+.xxxxx..\
+........\
+........\
+........\
+........";
+static char Char_LPAREN[80] = "\
+....x...\
+...x....\
+..x.....\
+.x......\
+.x......\
+.x......\
+..x.....\
+...x....\
+....x...\
+........";
+static char Char_RPAREN[80] = "\
+...x....\
+....x...\
+.....x..\
+......x.\
+......x.\
+......x.\
+.....x..\
+....x...\
+...x....\
+........";
+static char Char_QUOTE[80] = "\
+..x.x...\
+..x.x...\
+........\
+........\
+........\
+........\
+........\
+........\
+........\
+........";
+static char Char_COLON[80] = "\
+........\
+........\
+...xx...\
+...xx...\
+........\
+...xx...\
+...xx...\
+........\
+........\
+........";
+static char Char_PERIOD[80] = "\
+........\
+........\
+........\
+........\
+........\
+........\
+........\
+...xx...\
+...xx...\
+........";
+static char Char_COMMA[80] = "\
+........\
+........\
+........\
+........\
+........\
+........\
+...xx...\
+...xx...\
+....x...\
+...x....";
+
+static char Char_SLASH[80] = "\
+........\
+......x.\
+.....x..\
+....x...\
+...x....\
+..x.....\
+.x......\
+x.......\
+........\
+........";
+
+static char Char_SPACE[80] = "\
+........\
+........\
+........\
+........\
+........\
+........\
+........\
+........\
+........\
+........";
+
+static char  *GP_Font[128];
+static int    InitGP_Font = 0;
+
+static void initGP_Fonts(void) {
+
+  int i;
+  for (i = 0; i < 128; i++)
+    GP_Font[i] = (char *) 0;
+  
+  GP_Font['A'] = GP_Font['a'] = Char_A;
+  GP_Font['B'] = GP_Font['b'] = Char_B;
+  GP_Font['C'] = GP_Font['c'] = Char_C;
+  GP_Font['D'] = GP_Font['d'] = Char_D;
+  GP_Font['E'] = GP_Font['e'] = Char_E;
+  GP_Font['F'] = GP_Font['f'] = Char_F;
+  GP_Font['G'] = GP_Font['g'] = Char_G;
+  GP_Font['H'] = GP_Font['h'] = Char_H;
+  GP_Font['I'] = GP_Font['i'] = Char_I;
+  GP_Font['J'] = GP_Font['j'] = Char_J;
+  GP_Font['K'] = GP_Font['k'] = Char_K;
+  GP_Font['L'] = GP_Font['l'] = Char_L;
+  GP_Font['M'] = GP_Font['m'] = Char_M;
+  GP_Font['N'] = GP_Font['n'] = Char_N;
+  GP_Font['O'] = GP_Font['o'] = Char_O;
+  GP_Font['P'] = GP_Font['p'] = Char_P;
+  GP_Font['Q'] = GP_Font['q'] = Char_Q;
+  GP_Font['R'] = GP_Font['r'] = Char_R;
+  GP_Font['S'] = GP_Font['s'] = Char_S;
+  GP_Font['T'] = GP_Font['t'] = Char_T;
+  GP_Font['U'] = GP_Font['u'] = Char_U;
+  GP_Font['V'] = GP_Font['v'] = Char_V;
+  GP_Font['W'] = GP_Font['w'] = Char_W;
+  GP_Font['X'] = GP_Font['x'] = Char_X;
+  GP_Font['Y'] = GP_Font['y'] = Char_Y;
+  GP_Font['Z'] = GP_Font['z'] = Char_Z;
+  GP_Font['0'] = Char_0;
+  GP_Font['1'] = Char_1;
+  GP_Font['2'] = Char_2;
+  GP_Font['3'] = Char_3;
+  GP_Font['4'] = Char_4;
+  GP_Font['5'] = Char_5;
+  GP_Font['6'] = Char_6;
+  GP_Font['7'] = Char_7;
+  GP_Font['8'] = Char_8;
+  GP_Font['9'] = Char_9;
+  GP_Font['.'] = Char_PERIOD;
+  GP_Font[','] = Char_COMMA;
+  GP_Font['='] = Char_EQUAL;
+  GP_Font['-'] = Char_MINUS;
+  GP_Font['+'] = Char_PLUS;
+  GP_Font['\"'] = Char_QUOTE;
+  GP_Font['('] = Char_LPAREN;
+  GP_Font[')'] = Char_RPAREN;
+  GP_Font[':'] = Char_COLON;
+  GP_Font['/'] = Char_SLASH;
+  GP_Font[' '] = Char_SPACE;
+  InitGP_Font = 1;
+}
+
+/* -----------------------------------------------------------------------
+   void FrameBuffer_drawchar(FrameBuffer *f, int x, int y, Pixel fgcolor, Pixel bgcolor, char chr, int orientation)
+
+   Draws a character at location x, y with given color and orientation parameters.
+   Orientation can either be HORIZONTAL or VERTICAL
+   ----------------------------------------------------------------------- */
+void FrameBuffer_drawchar(FrameBuffer *f, int x, int y, int fgcolor,
+			  int bgcolor, char chr, int orientation) {
+
+  Pixel c, bc,*p,*p1;
+  char *ch;
+  int i,j;
+  int xpixels,ypixels;
+
+  if (!InitGP_Font) initGP_Fonts();
+  
+  c =  (Pixel) fgcolor;
+  bc = (Pixel) bgcolor;
+  xpixels = f->width;
+  ypixels = f->height;
+
+  if (orientation == HORIZONTAL) {
+    if ((x < f->xmin) || (x > f->xmax-8) ||
+	(y < f->ymin) || (y > f->ymax-10)) return;
+  
+    ch = GP_Font[(int) chr];
+    if (!ch) return;
+    p = &f->pixels[y+9][x];
+    for (i = 0; i < 10; i++) {
+      p1 = p;
+      for (j = 0; j< 8; j++) {
+	if (*(ch++) == 'x') *p= c;
+	else if (bgcolor >= 0)
+	  *p = bc;
+	p++;
+      }
+      p = (p1 - xpixels);
+    }
+  } else {
+    if ((x < f->xmin+10) || (x >= f->xmax) ||
+	(y < f->ymin) || (y > f->ymax-8)) return;
+  
+    ch = GP_Font[(int) chr];
+    if (!ch) return;
+    p = &f->pixels[y][x-9];
+    for (i = 0; i < 10; i++) {
+      p1 = p;
+      for (j = 0; j< 8; j++) {
+	if (*(ch++) == 'x') *p= c;
+	else if (bgcolor >= 0)
+	  *p = bc;
+	p+=xpixels;
+      }
+      p = p1 + 1;
+    }
+  }    
+}
+
+/* ----------------------------------------------------------------------
+   void FrameBuffer_drawstring(FrameBuffer *f, int x, int y, int fgcolor,
+                               int bgcolor, char *text, int orientation)
+
+   Draws an ASCII string on the framebuffer.   Can be oriented either horizontally
+   or vertically.  
+   ---------------------------------------------------------------------- */
+
+void FrameBuffer_drawstring(FrameBuffer *f, int x, int y, int fgcolor, int bgcolor, char *text, int orientation) {
+
+  char *c;
+  int  x1, y1;
+  int xpixels, ypixels;
+  
+  x1 = x;
+  y1 = y;
+  xpixels = f->width;
+  ypixels = f->height;
+  c = text;
+  while (*c) {
+    if (*c == '\n') {
+      if (orientation == HORIZONTAL) {
+	x1 = x; y1= y1- 10*xpixels;
+      } else {
+	y1 = y; x1= x1 + 10*ypixels;
+      }
+    } else {
+      FrameBuffer_drawchar(f, x1,y1,fgcolor, bgcolor,*c, orientation);
+      if (orientation == HORIZONTAL) {
+	x1+=8;
+	if (x1 >= (xpixels-8)) {
+	  x1 = x; y1= y1- 10;}
+	if (y1 < 0) return;
+      } else {
+	y1 += 8;
+	if (y1 >= (ypixels-8)) {
+	  y1 = y; x1 = x1 + 10;}
+	if (x1 > (xpixels-10)) return;
+      }
+    }
+    c++;
+  }
+}
+
+
+  
+  
+  
+  
+  
+
diff --git a/Examples/GIFPlot/Lib/frame.c b/Examples/GIFPlot/Lib/frame.c
new file mode 100644
index 0000000..e006f0d
--- /dev/null
+++ b/Examples/GIFPlot/Lib/frame.c
@@ -0,0 +1,924 @@
+/* ----------------------------------------------------------------------------- 
+ * frame.c
+ *
+ *     Frame buffer management
+ * 
+ * Author(s) : David Beazley (beazley@cs.uchicago.edu)
+ * Copyright (C) 1995-1996
+ *
+ * See the file LICENSE for information on usage and redistribution.	
+ * ----------------------------------------------------------------------------- */
+
+#define FRAME
+#include "gifplot.h"
+#include <float.h>
+
+/* ------------------------------------------------------------------------
+   FrameBuffer *new_FrameBuffer(int width, int height)
+
+   Creates a new framebuffer for storing data.
+   ------------------------------------------------------------------------ */
+
+FrameBuffer *new_FrameBuffer(unsigned int width, unsigned int height) {
+
+  FrameBuffer *f;
+  int  FrameBuffer_resize(FrameBuffer *f, int width, int height);
+
+  /* Create a new frame buffer */
+  
+  f = (FrameBuffer *) malloc(sizeof(FrameBuffer));
+  f->pixels = (Pixel **) 0;
+  f->zbuffer = (Zvalue **) 0;
+  /* Set its size */
+  
+  if (FrameBuffer_resize(f, width, height) == -1) {
+    free((char *) f);
+    return (FrameBuffer *) 0;
+  }
+
+  f->xmin = 0;
+  f->ymin = 0;
+  f->xmax = width;
+  f->ymax = height;
+  return f;
+}
+
+/* ------------------------------------------------------------------------
+   void delete_FrameBuffer(FrameBuffer *f)
+
+   Destroys the given framebuffer
+   ------------------------------------------------------------------------ */
+
+void delete_FrameBuffer(FrameBuffer *f) {
+
+  if (f) {
+    if (f->pixels) {
+      free((char *) f->pixels[0]);
+      free((char *) f->pixels);
+    }
+    if (f->zbuffer) {
+      free((char *) f->zbuffer[0]);
+      free((char *) f->zbuffer);
+    }
+    free((char *)f);
+  }
+}
+
+/* ------------------------------------------------------------------------
+   int *FrameBuffer_resize(FrameBuffer *f, int width, int height)
+
+   Resize the given framebuffer. Returns 0 on success, -1 on failure.
+   ------------------------------------------------------------------------ */
+
+int FrameBuffer_resize(FrameBuffer *f, int width, int height) {
+  int i;
+  if ((f) && (width > 0) && (height > 0)) {
+    if (f->pixels) {
+      free((char *)f->pixels[0]);
+      free((char *)f->pixels);
+    }
+    f->pixels = (Pixel **) malloc (height*sizeof(Pixel *));
+    if (!f->pixels) return -1;
+    f->pixels[0] = (Pixel *) malloc(height*width*sizeof(Pixel));
+    if (!f->pixels[0]) {
+      free((char *)f->pixels);
+      return -1;
+    }
+    for (i = 0; i < height; i++) 
+      f->pixels[i] = f->pixels[0] + i*width;
+    f->width = width;
+    f->height = height;
+    if (f->zbuffer) {
+      FrameBuffer_zresize(f,width,height);
+    }
+    return 0;
+  } else {
+    return -1;
+  }
+}
+
+/* ------------------------------------------------------------------------
+   void FrameBuffer_clear(FrameBuffer *f, Pixel color)
+
+   Clears the current FrameBuffer
+   ------------------------------------------------------------------------ */
+
+void FrameBuffer_clear(FrameBuffer *f, Pixel color) {
+  Pixel *p;
+  unsigned int i;
+  p = &f->pixels[0][0];
+
+  for (i = 0; i < f->width*f->height; i++, p++)
+    *p = color;
+}
+
+/* ------------------------------------------------------------------------
+   void FrameBuffer_plot(FrameBuffer *f, int x1, int y1, Pixel color)
+
+   Plots a point and does a bounds check.
+   ------------------------------------------------------------------------ */
+
+void FrameBuffer_plot(FrameBuffer *f, int x1, int y1, Pixel color) {
+
+  if ((x1 < f->xmin) || (x1 >= f->xmax) || (y1 < f->ymin) || (y1 >= f->ymax))
+    return;
+  f->pixels[y1][x1] = color;
+
+}
+
+/* ------------------------------------------------------------------------
+   FrameBuffer_horizontal(Framebuffer *f, int xmin, int xmax, int y, Pixel color)
+
+   Draw a horizontal line (clipped)
+   ------------------------------------------------------------------------ */
+
+void FrameBuffer_horizontal(FrameBuffer *f, int xmin, int xmax, int y, Pixel color) {
+  
+  Pixel *p;
+  int    i;
+  
+  if ((y < f->ymin) || (y >= f->ymax)) return;
+  if (xmin < f->xmin) xmin = f->xmin;
+  if (xmax >= f->xmax) xmax = f->xmax - 1;
+
+  p = &f->pixels[y][xmin];
+  for (i = xmin; i <= xmax; i++, p++)
+    *p = color;
+
+}
+
+/* ------------------------------------------------------------------------
+   FrameBuffer_horizontalinterp(Framebuffer *f, int xmin, int xmax, int y,
+                                Pixel c1, Pixel c2)
+
+   Draw a horizontal line (clipped) with color interpolation.
+   ------------------------------------------------------------------------ */
+
+void FrameBuffer_horizontalinterp(FrameBuffer *f, int xmin, int xmax, int y,
+				  Pixel c1, Pixel c2) {
+  
+  Pixel *p;
+  int    i;
+  double mc;
+  int    x1;
+  if ((y < f->ymin) || (y >= f->ymax)) return;
+
+  x1 = xmin;
+  if (xmin < f->xmin) xmin = f->xmin;
+  if (xmax >= f->xmax) xmax = f->xmax - 1;
+  if (xmax < f->xmin) return;
+  if (xmin >= f->xmax) return;
+
+  if (xmin != xmax)
+    mc = (double)(c2 - c1)/(double) (xmax - xmin);
+  else
+    mc = 0.0;
+
+  p = &f->pixels[y][xmin];
+  for (i = xmin; i <= xmax; i++, p++) 
+    *p = (Pixel) (mc*(i-x1) + c1);
+  
+}
+
+
+/* ------------------------------------------------------------------------
+   FrameBuffer_vertical(Framebuffer *f, int xmin, int xmax, int y, Pixel color)
+
+   Draw a Vertical line (clipped)
+   ------------------------------------------------------------------------ */
+
+void FrameBuffer_vertical(FrameBuffer *f, int ymin, int ymax, int x, Pixel color) {
+  
+  Pixel *p;
+  int    i;
+  
+  if ((x < f->xmin) || (x >= f->xmax)) return;
+  if (ymax < f->ymin) return;
+  if (ymin > f->ymax) return;
+  if (ymin < f->ymin) ymin = f->ymin;
+  if (ymax >= f->ymax) ymax = f->ymax - 1;
+
+  p = &f->pixels[ymin][x];
+  for (i = 0; i <= (ymax - ymin); i++, p+=f->width)
+    *p = color;
+
+}
+  
+/* ------------------------------------------------------------------------
+   void FrameBuffer_box(FrameBuffer *f, int x1, int y1, int x2, int y2, Pixel color)
+
+   Makes an outline box.
+   ------------------------------------------------------------------------ */
+
+void FrameBuffer_box(FrameBuffer *f, int x1, int y1, int x2, int y2, Pixel color) {
+
+  int  xt, yt;
+
+  /* Make sure points are in correct order */
+
+  if (x2 < x1) {
+    xt = x2;
+    x2 = x1;
+    x1 = xt;
+  }
+  if (y2 < y1) {
+    yt = y2;
+    y2 = y1;
+    y1 = yt;
+  }
+
+  /* Draw lower edge */
+
+  FrameBuffer_horizontal(f,x1,x2,y1,color);
+
+  /* Draw upper edge */
+
+  FrameBuffer_horizontal(f,x1,x2,y2,color);
+
+  /* Draw left side */
+
+  FrameBuffer_vertical(f,y1,y2,x1,color);
+
+  /* Draw right side */
+
+  FrameBuffer_vertical(f,y1,y2,x2,color);
+  
+}
+
+/* ------------------------------------------------------------------------
+   void FrameBuffer_solidbox(FrameBuffer *f, int x1, int y1, int x2, int y2, Pixel color)
+
+   Makes an solid box.
+   ------------------------------------------------------------------------ */
+
+void FrameBuffer_solidbox(FrameBuffer *f, int x1, int y1, int x2, int y2, Pixel color) {
+
+  int  xt, yt;
+
+  /* Make sure points are in correct order */
+
+  if (x2 < x1) {
+    xt = x2;
+    x2 = x1;
+    x1 = xt;
+  }
+  if (y2 < y1) {
+    yt = y2;
+    y2 = y1;
+    y1 = yt;
+  }
+
+  /* Now perform some clipping */
+
+  if (y1 < f->ymin) y1 = f->ymin;
+  if (y2 >= f->ymax) y2 = f->ymax - 1;
+
+  /* Fill it in using horizontal lines */
+  
+  for (yt = y1; yt <= y2; yt++)
+    FrameBuffer_horizontal(f,x1,x2,yt,color);
+
+}
+
+/* ------------------------------------------------------------------------
+   void FrameBuffer_interpbox(FrameBuffer *f, int x1, int y1, int x2, int y2
+                              Pixel c1, Pixel c2, Pixel c3, Pixel c4)
+
+   Makes a box with interpolated color.   Colors are assigned as follows :
+   (x1,y1) = c1
+   (x1,y2) = c2
+   (x2,y1) = c3
+   (x2,y2) = c4
+   ------------------------------------------------------------------------ */
+
+void FrameBuffer_interpbox(FrameBuffer *f, int x1, int y1, int x2, int y2,
+			   Pixel c1, Pixel c2, Pixel c3, Pixel c4) {
+
+  int  xt, yt;
+  Pixel ct;
+  double  mc1,mc2;
+  int  ystart;
+  /* Make sure points are in correct order */
+
+  if (x2 < x1) {
+    xt = x2;
+    x2 = x1;
+    x1 = xt;
+    ct = c1;
+    c1 = c3;
+    c3 = ct;
+    ct = c2;
+    c2 = c4;
+    c4 = ct;
+  }
+  if (y2 < y1) {
+    yt = y2;
+    y2 = y1;
+    y1 = yt;
+    ct = c1;
+    c1 = c2;
+    c2 = ct;
+    ct = c3;
+    c3 = c4;
+    c4 = ct;
+  }
+
+  /* Now perform some clipping */
+
+  ystart = y1;
+  mc1 = (double) (c2 - c1)/(double) (y2 - y1);
+  mc2 = (double) (c4 - c3)/(double) (y2 - y1);
+  if (y1 < f->ymin) y1 = f->ymin;
+  if (y2 >= f->ymax) y2 = f->ymax - 1;
+
+  /* Fill it in using horizontal lines */
+  
+  for (yt = y1; yt <= y2; yt++)
+    FrameBuffer_horizontalinterp(f,x1,x2,yt,(Pixel) ((mc1*(yt - ystart)) + c1),
+				  (Pixel) ((mc2*(yt-ystart))+c3));
+
+}
+
+/* ---------------------------------------------------------------------------
+   FrameBuffer_line(FrameBuffer *f, int x1, int y1, int x2, int y2, color)
+
+   Draws a line on the framebuffer using the Bresenham line algorithm.  The
+   line is clipped to fit within the current view window.
+   ---------------------------------------------------------------------------- */
+
+void FrameBuffer_line(FrameBuffer *f, int x1, int y1, int x2, int y2, Pixel c) {
+
+  int  dx,dy,dxneg,dyneg, inc1,inc2,di;
+  int  x, y, xpixels, ypixels, xt, yt;
+  Pixel *p;
+  double m;
+  int  end1 = 0, end2 = 0;
+  
+  /* Need to figure out where in the heck this line is */
+
+  dx = x2 - x1;
+  dy = y2 - y1;
+
+  if (dx == 0) {
+    /* Draw a Vertical Line */
+    if (y1 < y2)
+      FrameBuffer_vertical(f,y1,y2,x1,c);
+    else
+      FrameBuffer_vertical(f,y2,y1,x1,c);
+    return;
+  }
+  if (dy == 0) {
+    /* Draw a Horizontal Line */
+    if (x1 < x2)
+      FrameBuffer_horizontal(f,x1,x2,y1,c);
+    else
+      FrameBuffer_horizontal(f,x2,x1,y1,c);
+    return;
+  }
+
+  /* Figure out where in the heck these lines are using the
+     Cohen-Sutherland Line Clipping Scheme. */
+
+  end1 = ((x1 - f->xmin) < 0) |
+    (((f->xmax- 1 - x1) < 0) << 1) |
+    (((y1 - f->ymin) < 0) << 2) |
+    (((f->ymax-1 - y1) < 0) << 3);
+
+  end2 = ((x2 - f->xmin) < 0) |
+    (((f->xmax-1 - x2) < 0) << 1) |
+    (((y2 - f->ymin) < 0) << 2) |
+    (((f->ymax-1 - y2) < 0) << 3);
+
+  if (end1 & end2) return;      /* Nope : Not visible */
+
+  /* Make sure points have a favorable orientation */
+  
+  if (x1 > x2) {
+    xt = x1;
+    x1 = x2;
+    x2 = xt;
+    yt = y1;
+    y1 = y2;
+    y2 = yt;
+  }
+  
+  /* Clip against the boundaries */
+  m = (y2 - y1)/(double) (x2-x1);
+  if (x1 < f->xmin) {
+    y1 = (int) ((f->xmin - x1)*m + y1);
+    x1 = (int) f->xmin;
+  }
+  if (x2 >= f->xmax) {
+    y2 = (int) ((f->xmax -1 -x1)*m + y1);
+    x2 = (int) (f->xmax - 1);
+  }
+
+  if (y1 > y2) {
+    xt = x1;
+    x1 = x2;
+    x2 = xt;
+    yt = y1;
+    y1 = y2;
+    y2 = yt;
+  }
+
+  m = 1/m;
+  if (y1 < f->ymin) {
+    x1 = (int) ((f->ymin - y1)*m + x1);
+    y1 = (int) f->ymin;
+  }
+  if (y2 >= f->ymax) {
+    x2 = (int) ((f->ymax-1-y1)*m + x1);
+    y2 = (int) (f->ymax-1);
+  }
+
+  if ((x1 < f->xmin) || (x1 >= f->xmax) || (y1 < f->ymin) || (y1 >= f->ymax) ||
+      (x2 < f->xmin) || (x2 >= f->xmax) || (y2 < f->ymin) || (y2 >= f->ymax)) return;
+
+  dx = x2 - x1;
+  dy = y2 - y1;
+  xpixels = f->width;
+  ypixels = f->height;
+  
+  dxneg = (dx < 0) ? 1 : 0;
+  dyneg = (dy < 0) ? 1 : 0;
+  
+  dx = abs(dx);
+  dy = abs(dy);
+  if (dx >= dy) {
+    /* Slope between -1 and 1. */
+    if (dxneg) {
+      x = x1;
+      y = y1;
+      x1 = x2;
+      y1 = y2;
+      x2 = x;
+      y2 = y;
+      dyneg = !dyneg;
+    }
+    inc1 = 2*dy;
+    inc2 = 2*(dy-dx);
+    di = 2*dy-dx;
+
+    /* Draw a line using x as independent variable */
+    
+    p = &f->pixels[y1][x1];
+    x = x1;
+    while (x <= x2) {
+      *(p++) = c;
+      if (di < 0) {
+	di = di + inc1;
+      } else {
+	if (dyneg) {
+	  p = p - xpixels;
+	  di = di + inc2;
+	} else {
+	  p = p + xpixels;
+	  di = di + inc2;
+	}
+      }
+      x++;
+    }
+  } else {
+    /* Slope < -1 or > 1 */
+    if (dyneg) {
+      x = x1;
+      y = y1;
+      x1 = x2;
+      y1 = y2;
+      x2 = x;
+      y2 = y;
+      dxneg = !dxneg;
+    }
+    inc1 = 2*dx;
+    inc2 = 2*(dx-dy);
+    di = 2*dx-dy;
+
+    /* Draw a line using y as independent variable */
+    
+    p = &f->pixels[y1][x1];
+    y = y1;
+    while (y <= y2) {
+      *p = c;
+      p = p + xpixels;
+      if (di < 0) {
+	di = di + inc1;
+      } else {
+	if (dxneg) {
+	  p = p - 1;
+	  di = di + inc2;
+	} else {
+	  p = p + 1;
+	  di = di + inc2;
+	}
+      }
+      y++;
+    }
+  }    
+}
+
+
+/* -------------------------------------------------------------------------
+   FrameBuffer_circle(FrameBuffer f, int xc, int yc, int radius, Pixel c)
+
+   Create an outline circle
+   ------------------------------------------------------------------------- */
+
+#define plot_circle(x,y,c) \
+   if ((x >= xmin) && (x < xmax) && \
+       (y >= ymin) && (y < ymax)) \
+        pixels[y][x] = c;
+	
+void FrameBuffer_circle(FrameBuffer *f, int xc, int yc, int radius, Pixel c) {
+
+  int xpixels, ypixels, x, y, p;
+  int xmin, ymin, xmax, ymax;
+  Pixel **pixels;
+
+  if (radius <= 0) return;
+  xpixels = f->width;
+  ypixels = f->height;
+  pixels = f->pixels;
+  xmin = f->xmin;
+  ymin = f->ymin;
+  xmax = f->xmax;
+  ymax = f->ymax;
+  x = 0;
+  y = radius;
+  p = 3-2*radius;
+  while (x <= y) {
+    plot_circle(xc+x,yc+y,c);
+    plot_circle(xc-x,yc+y,c);
+    plot_circle(xc+x,yc-y,c);
+    plot_circle(xc-x,yc-y,c);
+    plot_circle(xc+y,yc+x,c);
+    plot_circle(xc-y,yc+x,c);
+    plot_circle(xc+y,yc-x,c);
+    plot_circle(xc-y,yc-x,c);
+    if (p < 0) p = p + 4*x + 6;
+    else {
+      p = p + 4*(x-y) + 10;
+      y = y -1;
+    }
+    x++;
+  }
+}
+
+
+/* -------------------------------------------------------------------------
+   FrameBuffer_solidcircle(FrameBuffer f, int xc, int yc, int radius, Pixel c)
+
+   Create an filled circle
+   ------------------------------------------------------------------------- */
+
+       
+#define fill_circle(x,y,c) \
+       x1 = xc - x; \
+       x2 = xc + x; \
+       FrameBuffer_horizontal(f,x1,x2,y,c);
+	
+void FrameBuffer_solidcircle(FrameBuffer *f, int xc, int yc, int radius, Pixel c) {
+
+  int xpixels, ypixels, x, y, p;
+  int x1,x2;
+  int xmin, ymin, xmax, ymax;
+  Pixel **pixels;
+
+  if (radius <= 0) return;
+  xpixels = f->width;
+  ypixels = f->height;
+  pixels = f->pixels;
+  xmin = f->xmin;
+  ymin = f->ymin;
+  xmax = f->xmax;
+  ymax = f->ymax;
+  x = 0;
+  y = radius;
+  p = 3-2*radius;
+  while (x <= y) {
+    fill_circle(x,yc+y,c);
+    fill_circle(x,yc-y,c);
+    fill_circle(y,yc+x,c);
+    fill_circle(y,yc-x,c);
+    if (p < 0) p = p + 4*x + 6;
+    else {
+      p = p + 4*(x-y) + 10;
+      y = y -1;
+    }
+    x++;
+  }
+}
+
+/* ------------------------------------------------------------------------
+   void FrameBuffer_setclip(f,xmin,ymin,xmax,ymax)
+
+   Set clipping region for plotting
+   ------------------------------------------------------------------------ */
+
+void FrameBuffer_setclip(FrameBuffer *f, int xmin, int ymin, int xmax, int ymax) {
+
+  if (xmin >= xmax) return;
+  if (ymin >= ymax) return;
+
+  if (xmin < 0) xmin = 0;
+  if (ymin < 0) ymin = 0;
+  if (xmax > (int) f->width) xmax = f->width;
+  if (ymax > (int) f->height) ymax = f->height;
+
+  f->xmin = xmin;
+  f->ymin = ymin;
+  f->xmax = xmax;
+  f->ymax = ymax;
+}
+
+/* ------------------------------------------------------------------------
+   void FrameBuffer_noclip(f)
+
+   Disable clipping region
+   ------------------------------------------------------------------------ */
+
+void FrameBuffer_noclip(FrameBuffer *f) {
+  f->xmin = 0;
+  f->ymin = 0;
+  f->xmax = f->width;
+  f->ymax = f->height;
+}
+
+
+/* ------------------------------------------------------------------------
+   FrameBuffer_zresize(FrameBuffer *f, int width, int height)
+
+   This function resizes the framebuffer's zbuffer.  If none exist, it
+   creates a new one.
+   ------------------------------------------------------------------------ */
+
+void FrameBuffer_zresize(FrameBuffer *f, int width, int height) {
+  int i;
+
+  if (f->zbuffer) {
+    free((char *)f->zbuffer[0]);
+    free((char *)f->zbuffer);
+  }
+  f->zbuffer = (Zvalue **) malloc(height*sizeof(Zvalue *));
+  f->zbuffer[0] = (Zvalue *) malloc(height*width*sizeof(Zvalue));
+  for (i = 0; i < height; i++)
+    f->zbuffer[i] = f->zbuffer[0]+i*width;
+}
+
+/* ------------------------------------------------------------------------
+   FrameBuffer_zclear(FrameBuffer *f)
+
+   Clears the z-buffer for a particular frame.   Sets all of the z-values to
+   ZMIN.
+   ------------------------------------------------------------------------- */
+
+void FrameBuffer_zclear(FrameBuffer *f) {
+  unsigned int i,j;
+  if (f) {
+    if (f->zbuffer) {
+      for (i = 0; i < f->width; i++)
+	for (j = 0; j < f->height; j++)
+	  f->zbuffer[j][i] = ZMIN;
+    }
+  }
+}
+   
+
+
+/* -------------------------------------------------------------------------
+   FrameBuffer_solidtriangle(FrameBuffer *f, int tx1, int ty2, 
+                                 int tx2, int ty2, 
+				 int tx3, int ty3, Pixel color) 
+
+   This function draws a 2D filled triangle.   
+
+   General idea :
+         1.   Transform the three points into screen coordinates
+	 2.   Order three points vertically on screen.
+	 3.   Check for degenerate cases (where 3 points are colinear).
+	 4.   Fill in the resulting triangle using horizontal lines.
+   -------------------------------------------------------------------------- */
+
+void FrameBuffer_solidtriangle(FrameBuffer *f, int tx1, int ty1, 
+			       int tx2, int ty2,
+			       int tx3, int ty3, Pixel color) {
+  int        tempx, tempy;
+  double     m1,m2,m3;
+  int        y;
+  int        ix1, ix2;
+
+  /* Figure out which point has the greatest "y" value */
+
+  if (ty2 > ty1) {   /* Swap points 1 and 2 if 2 is higher */
+    tempx = tx1;
+    tempy = ty1;
+    tx1 = tx2;
+    ty1 = ty2;
+    tx2 = tempx;
+    ty2 = tempy;
+  }
+  if (ty3 > ty1) {  /* Swap points 1 and 3 if 3 is higher */
+    tempx = tx1;
+    tempy = ty1;
+    tx1 = tx3;
+    ty1 = ty3;
+    tx3 = tempx;
+    ty3 = tempy;
+  }
+  if (ty3 > ty2) { /* Swap points 2 and 3 if 3 is higher */ 
+    tempx = tx2;
+    tempy = ty2;
+    tx2 = tx3;
+    ty2 = ty3;
+    tx3 = tempx;
+    ty3 = tempy;
+  }
+
+  /* Points are now order so that t_1 is the highest point, t_2 is the
+     middle point, and t_3 is the lowest point */
+
+  /* Check for degenerate cases here */
+
+  if ((ty1 == ty2) && (ty2 == ty3)) {
+
+    /* Points are aligned horizontally.   Handle as a special case */
+    /* Just draw three lines using the outline color */
+
+    FrameBuffer_line(f,tx1,ty1,tx2,ty2,color);
+    FrameBuffer_line(f,tx1,ty1,tx3,ty3,color);
+    FrameBuffer_line(f,tx2,ty2,tx3,ty3,color);
+
+  } else {
+
+    if (ty2 < ty1) {
+      /* First process line segments between (x1,y1)-(x2,y2)
+	 And between (x1,y1),(x3,y3) */
+
+      m1 = (double) (tx2 - tx1)/(double) (ty2 - ty1);
+      m2 = (double) (tx3 - tx1)/(double) (ty3 - ty1);
+
+      y = ty1;
+      while (y >= ty2) {
+	/* Calculate x values from slope */
+	ix1 = (int) (m1*(y-ty1)+0.5) + tx1;
+	ix2 = (int) (m2*(y-ty1)+0.5) + tx1;
+	if (ix1 > ix2) 
+	  FrameBuffer_horizontal(f,ix2,ix1,y,color);
+	else
+	  FrameBuffer_horizontal(f,ix1,ix2,y,color);
+	y--;
+      }
+    }
+    if (ty3 < ty2) {
+      /* Draw lower half of the triangle */
+      m2 = (double) (tx3 - tx1)/(double) (ty3 - ty1);
+      m3 = (double) (tx3 - tx2)/(double)(ty3 - ty2);
+      y = ty2;
+      while (y >= ty3) {
+	ix1 = (int) (m3*(y-ty2)+0.5)+tx2;
+	ix2 = (int) (m2*(y-ty1)+0.5)+tx1;
+	if (ix1 > ix2)
+	  FrameBuffer_horizontal(f,ix2,ix1,y,color);
+	else
+	  FrameBuffer_horizontal(f,ix1,ix2,y,color);
+	y--;
+      }
+    }
+  }
+}
+
+/* -------------------------------------------------------------------------
+   FrameBuffer_interptriangle(FrameBuffer *f,
+   int tx1, int ty2, Pixel c1,
+   int tx2, int ty2, Pixel c2,
+   int tx3, int ty3, Pixel c3)
+
+   This function draws a filled triangle with color
+   interpolation.  
+
+   General idea :
+         1.   Transform the three points into screen coordinates
+	 2.   Order three points vertically on screen.
+	 3.   Check for degenerate cases (where 3 points are colinear).
+	 4.   Fill in the resulting triangle using horizontal lines.
+	 5.   Colors are interpolated between end points
+   -------------------------------------------------------------------------- */
+
+void FrameBuffer_interptriangle(FrameBuffer *f,
+				int tx1, int ty1, Pixel c1,
+				int tx2, int ty2, Pixel c2,
+				int tx3, int ty3, Pixel c3) {
+  int        tempx, tempy;
+  double     m1,m2,m3;
+  double     mc1,mc2,mc3;
+  Pixel      ic1,ic2,tempc;
+  int        y;
+  int        ix1, ix2;
+
+  /* Figure out which point has the greatest "y" value */
+
+  if (ty2 > ty1) {   /* Swap points 1 and 2 if 2 is higher */
+    tempx = tx1;
+    tempy = ty1;
+    tempc = c1;
+    tx1 = tx2;
+    ty1 = ty2;
+    c1 = c2;
+    tx2 = tempx;
+    ty2 = tempy;
+    c2 = tempc;
+  }
+  if (ty3 > ty1) {  /* Swap points 1 and 3 if 3 is higher */
+    tempx = tx1;
+    tempy = ty1;
+    tempc = c1;
+    tx1 = tx3;
+    ty1 = ty3;
+    c1 = c3;
+    tx3 = tempx;
+    ty3 = tempy;
+    c3 = tempc;
+  }
+  if (ty3 > ty2) { /* Swap points 2 and 3 if 3 is higher */ 
+    tempx = tx2;
+    tempy = ty2;
+    tempc = c2;
+    tx2 = tx3;
+    ty2 = ty3;
+    c2 = c3;
+    tx3 = tempx;
+    ty3 = tempy;
+    c3 = tempc;
+  }
+
+  /* Points are now order so that t_1 is the highest point, t_2 is the
+     middle point, and t_3 is the lowest point */
+
+  /* Check for degenerate cases here */
+
+  if ((ty1 == ty2) && (ty2 == ty3)) {
+
+    /* Points are aligned horizontally.   Handle as a special case */
+    /* Just draw three lines using the outline color */
+
+    if (tx2 > tx1)
+      FrameBuffer_horizontalinterp(f,tx1,tx2,ty1,c1,c2);
+    else
+      FrameBuffer_horizontalinterp(f,tx2,tx1,ty1,c2,c1);
+    if (tx3 > tx1)
+      FrameBuffer_horizontalinterp(f,tx1,tx3,ty1,c1,c3);
+    else
+      FrameBuffer_horizontalinterp(f,tx3,tx1,ty1,c3,c1);
+    if (tx3 > tx2)
+      FrameBuffer_horizontalinterp(f,tx2,tx3,ty2,c2,c3);
+    else
+      FrameBuffer_horizontalinterp(f,tx3,tx2,ty2,c3,c2);
+    
+  } else {
+
+    /* First process line segments between (x1,y1)-(x2,y2)
+       And between (x1,y1),(x3,y3) */
+
+    if (ty2 < ty1) {
+      m1 = (double) (tx2 - tx1)/(double) (ty2 - ty1);
+      m2 = (double) (tx3 - tx1)/(double) (ty3 - ty1);
+      mc1 = (c2 - c1)/(double) (ty2 - ty1);
+      mc2 = (c3 - c1)/(double) (ty3 - ty1);
+
+      y = ty1;
+      while (y >= ty2) {
+	/* Calculate x values from slope */
+	ix1 = (int) (m1*(y-ty1)+0.5) + tx1;
+	ix2 = (int) (m2*(y-ty1)+0.5) + tx1;
+	ic1 = (int) (mc1*(y-ty1) + c1);
+	ic2 = (int) (mc2*(y-ty1) + c1);
+	if (ix1 > ix2) 
+	  FrameBuffer_horizontalinterp(f,ix2,ix1,y,ic2,ic1);
+	else
+	  FrameBuffer_horizontalinterp(f,ix1,ix2,y,ic1,ic2);
+	y--;
+      }
+    }
+    if (ty3 < ty2) {
+      /* Draw lower half of the triangle */
+      m2 = (double) (tx3 - tx1)/(double) (ty3 - ty1);
+      mc2 = (c3 - c1)/(double) (ty3 - ty1);
+      m3 = (double) (tx3 - tx2)/(double)(ty3 - ty2);
+      mc3 = (c3 - c2)/(double) (ty3 - ty2);
+      y = ty2;
+      while (y >= ty3) {
+	ix1 = (int) (m3*(y-ty2)+0.5)+tx2;
+	ix2 = (int) (m2*(y-ty1)+0.5)+tx1;
+	ic1 = (int) (mc3*(y-ty2)+c2);
+	ic2 = (int) (mc2*(y-ty1)+c1);
+	if (ix1 > ix2)
+	  FrameBuffer_horizontalinterp(f,ix2,ix1,y,ic2,ic1);
+	else
+	  FrameBuffer_horizontalinterp(f,ix1,ix2,y,ic1,ic2);
+	y--;
+      }
+    }
+  }
+}
+
+
diff --git a/Examples/GIFPlot/Lib/gif.c b/Examples/GIFPlot/Lib/gif.c
new file mode 100644
index 0000000..7953e5c
--- /dev/null
+++ b/Examples/GIFPlot/Lib/gif.c
@@ -0,0 +1,672 @@
+
+/**********************************************************************
+ * GIFPlot 0.0
+ * 
+ * Dave Beazley
+ * 
+ * Department of Computer Science        Theoretical Division (T-11)        
+ * University of Utah                    Los Alamos National Laboratory     
+ * Salt Lake City, Utah 84112            Los Alamos, New Mexico  87545      
+ * beazley@cs.utah.edu                   beazley@lanl.gov                   
+ *
+ * Copyright (c) 1996
+ * The Regents of the University of California and the University of Utah
+ * All Rights Reserved
+ *
+ * Permission is hereby granted, without written agreement and without
+ * license or royalty fees, to use, copy, modify, and distribute this
+ * software and its documentation for any purpose, provided that 
+ * (1) The above copyright notice and the following two paragraphs
+ * appear in all copies of the source code and (2) redistributions
+ * including binaries reproduces these notices in the supporting
+ * documentation.   Substantial modifications to this software may be
+ * copyrighted by their authors and need not follow the licensing terms
+ * described here, provided that the new terms are clearly indicated in
+ * all files where they apply.
+ * 
+ * IN NO EVENT SHALL THE AUTHOR, THE UNIVERSITY OF CALIFORNIA, THE 
+ * UNIVERSITY OF UTAH OR DISTRIBUTORS OF THIS SOFTWARE BE LIABLE TO ANY
+ * PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL
+ * DAMAGES ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION,
+ * EVEN IF THE AUTHORS OR ANY OF THE ABOVE PARTIES HAVE BEEN ADVISED OF
+ * THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * THE AUTHOR, THE UNIVERSITY OF CALIFORNIA, AND THE UNIVERSITY OF UTAH
+ * SPECIFICALLY DISCLAIM ANY WARRANTIES,INCLUDING, BUT NOT LIMITED TO, 
+ * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND 
+ * THE AUTHORS AND DISTRIBUTORS HAVE NO OBLIGATION TO PROVIDE MAINTENANCE,
+ * SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+ *
+ **************************************************************************/
+
+/*******************************************************************
+ * Creates a GIF format file.  
+ *
+ * Dave Beazley (T-11)
+ * August 11, 1995
+ *
+ * Rather than writing directly to files, this module fills out
+ * output buffer.
+ * 
+ * Note : To save memory, this routine uses approximately 50K of the
+ * output buffer as temporary storage (for hash tables and compression codes).
+ * The remainder of the output buffer is used to store the final image.
+ * This feature allows GIF images to be created with no additional
+ * memory overhead.
+ *
+ * -- Revision History 
+ * $Log$
+ * Revision 1.2  2003/09/01 16:23:31  beazley
+ * Restored the 'mojo'.
+ *
+ * Revision 1.2  1996/09/25 22:39:30  dmb
+ * Fixed prototypes and use of void pointers for compatibility with the Cray T3D
+ *
+ * Revision 1.1  1996/09/10 17:44:00  dmb
+ * Initial revision
+ *
+ * Revision 1.2  1995/08/31 14:46:07  beazley
+ * Minor changes to support comments and a few bug fixes.
+ *
+ *
+ ******************************************************************/
+
+
+/*
+ * xvgifwr.c  -  handles writing of GIF files.  based on flgife.c and
+ *               flgifc.c from the FBM Library, by Michael Maudlin
+ *
+ * Contains: 
+ *   WriteGIF(fp, pic, ptype, w, h, rmap, gmap, bmap, numcols, colorstyle,
+ *            comment)
+ *
+ * Note: slightly brain-damaged, in that it'll only write non-interlaced 
+ *       GIF files (in the interests of speed, or something)
+ *
+ */
+
+
+
+/*****************************************************************
+ * Portions of this code Copyright (C) 1989 by Michael Mauldin.
+ * Permission is granted to use this file in whole or in
+ * part for any purpose, educational, recreational or commercial,
+ * provided that this copyright notice is retained unchanged.
+ * This software is available to all free of charge by anonymous
+ * FTP and in the UUNET archives.
+ *
+ *
+ * Authors:  Michael Mauldin (mlm@cs.cmu.edu)
+ *           David Rowley (mgardi@watdcsu.waterloo.edu)
+ *
+ * Based on: compress.c - File compression ala IEEE Computer, June 1984.
+ *
+ *	Spencer W. Thomas       (decvax!harpo!utah-cs!utah-gr!thomas)
+ *	Jim McKie               (decvax!mcvax!jim)
+ *	Steve Davies            (decvax!vax135!petsd!peora!srd)
+ *	Ken Turkowski           (decvax!decwrl!turtlevax!ken)
+ *	James A. Woods          (decvax!ihnp4!ames!jaw)
+ *	Joe Orost               (decvax!vax135!petsd!joe)
+ *****************************************************************/
+
+#include "gifplot.h"
+#include <string.h> 
+typedef long int        count_int;
+typedef unsigned char   byte;
+
+static int  gif_error;
+static unsigned char *op;
+static int  Width, Height;
+static int  curx, cury;
+static int  Interlace;
+
+static void putgifword(int);
+static void compress(int, byte **, int);
+static void output_GIF(int);
+static void cl_block(void);
+static void cl_hash(count_int);
+static void char_init(void);
+static void char_out(int);
+static void flush_char(void);
+static void *OutBuffer;
+static int   OutBufSize;
+static FrameBuffer *GIF_frame;
+
+static unsigned char pc2nc[256],r1[256],g1[256],b1[256];
+
+/*************************************************************/
+int FrameBuffer_makeGIF(FrameBuffer *f, ColorMap *c, void *outbuffer, unsigned int outbufsize)
+{
+  int   RWidth, RHeight;
+  int   LeftOfs, TopOfs;
+  int   ColorMapSize, InitCodeSize, Background, BitsPerPixel;
+  int   i,j,nc;
+  char *rmap, *gmap, *bmap;
+  char *cmap;
+  int   count;
+
+  Interlace = 0;
+  Background = 0;
+  OutBuffer = outbuffer;
+  OutBufSize = outbufsize;
+  GIF_frame = f;
+  cmap = (char *) c->cmap;
+  
+  op = (unsigned char *) outbuffer;
+  gif_error = 0;
+  for (i=0; i<256; i++) { pc2nc[i] = r1[i] = g1[i] = b1[i] = 0; }
+
+  /* compute number of unique colors */
+  nc = 0;
+  rmap = &cmap[0];
+  gmap = &cmap[256];
+  bmap = &cmap[512];
+
+  for (i=0; i<256; i++) {
+    /* see if color #i is already used */
+    for (j=0; j<i; j++) {
+      if (rmap[i] == rmap[j] && gmap[i] == gmap[j] && 
+	  bmap[i] == bmap[j]) break;
+    }
+
+    if (j==i) {  /* wasn't found */
+      pc2nc[i] = nc;
+      r1[nc] = rmap[i];
+      g1[nc] = gmap[i];
+      b1[nc] = bmap[i];
+      nc++;
+    }
+    else pc2nc[i] = pc2nc[j];
+  }
+
+  /* figure out 'BitsPerPixel' */
+  for (i=1; i<8; i++)
+    if ( (1<<i) >= nc) break;
+  
+  BitsPerPixel = i;
+
+  ColorMapSize = 1 << BitsPerPixel;
+	
+  RWidth  = Width  = f->width;
+  RHeight = Height = f->height;
+  LeftOfs = TopOfs = 0;
+
+  if (BitsPerPixel <= 1) InitCodeSize = 2;
+                    else InitCodeSize = BitsPerPixel;
+
+  curx = 0;
+  cury = f->height - 1;
+
+  strcpy((char *) op,"GIF89a");        /* Put in GIF magic number */
+  op+=6;
+  putgifword(RWidth);           /* screen descriptor */
+  putgifword(RHeight);
+
+  i = 0x80;	                 /* Yes, there is a color map */
+  i |= (8-1)<<4;                 /* OR in the color resolution (hardwired 8) */
+  i |= (BitsPerPixel - 1);       /* OR in the # of bits per pixel */
+  *(op++) = i;
+  *(op++) = Background;          /* background color */
+  *(op++) = 0;
+  for (i=0; i<ColorMapSize; i++) {       /* write out Global colormap */
+    *(op++) = r1[i];
+    *(op++) = g1[i];
+    *(op++) = b1[i];
+  }
+
+  *(op++) = ',';                /* image separator */
+
+  /* Write the Image header */
+  putgifword(LeftOfs);
+  putgifword(TopOfs);
+  putgifword(Width);
+  putgifword(Height);
+  *(op++) = 0;
+  *(op++) = InitCodeSize;
+
+  compress(InitCodeSize+1, f->pixels, f->width*f->height);
+
+  *(op++) = 0;
+  *(op++) = ';';
+  
+  count = (op - (unsigned char *) OutBuffer);
+  if (gif_error) return -1;
+  else return count;
+}
+
+/******************************/
+static void putgifword(w)
+int w;
+{
+  /* writes a 16-bit integer in GIF order (LSB first) */
+  *(op++) = w & 0xff;
+  *(op++) = (w>>8)&0xff;
+}
+
+/***********************************************************************/
+
+
+static unsigned long cur_accum = 0;
+static int           cur_bits = 0;
+
+
+
+
+#define GP_BITS	12    /* BITS was already defined on some systems */
+
+#define HSIZE  5003            /* 80% occupancy */
+
+typedef unsigned char   char_type;
+
+static int n_bits;                    /* number of bits/code */
+static int maxbits = GP_BITS;         /* user settable max # bits/code */
+static int maxcode;                   /* maximum code, given n_bits */
+static int maxmaxcode = 1 << GP_BITS; /* NEVER generate this */
+
+#define MAXCODE(n_bits)     ( (1 << (n_bits)) - 1)
+
+static  count_int      *htab;
+static  unsigned short *codetab;
+static  int GIFOutBufSize;
+
+/* static  count_int      htab [HSIZE];
+static  unsigned short codetab [HSIZE]; */
+
+#define HashTabOf(i)   htab[i]
+#define CodeTabOf(i)   codetab[i]
+
+static int hsize = HSIZE;            /* for dynamic table sizing */
+
+/*
+ * To save much memory, we overlay the table used by compress() with those
+ * used by decompress().  The tab_prefix table is the same size and type
+ * as the codetab.  The tab_suffix table needs 2**BITS characters.  We
+ * get this from the beginning of htab.  The output stack uses the rest
+ * of htab, and contains characters.  There is plenty of room for any
+ * possible stack (stack used to be 8000 characters).
+ */
+
+#define tab_prefixof(i) CodeTabOf(i)
+#define tab_suffixof(i)        ((char_type *)(htab))[i]
+#define de_stack               ((char_type *)&tab_suffixof(1<<GP_BITS))
+
+static int free_ent = 0;                  /* first unused entry */
+
+/*
+ * block compression parameters -- after all codes are used up,
+ * and compression rate changes, start over.
+ */
+static int clear_flg = 0;
+
+static long int out_count = 0;           /* # of codes output (for debugging) */
+
+/*
+ * compress stdin to stdout
+ *
+ * Algorithm:  use open addressing double hashing (no chaining) on the 
+ * prefix code / next character combination.  We do a variant of Knuth's
+ * algorithm D (vol. 3, sec. 6.4) along with G. Knott's relatively-prime
+ * secondary probe.  Here, the modular division first probe is gives way
+ * to a faster exclusive-or manipulation.  Also do block compression with
+ * an adaptive reset, whereby the code table is cleared when the compression
+ * ratio decreases, but after the table fills.  The variable-length output
+ * codes are re-sized at this point, and a special CLEAR code is generated
+ * for the decompressor.  Late addition:  construct the table according to
+ * file size for noticeable speed improvement on small files.  Please direct
+ * questions about this implementation to ames!jaw.
+ */
+
+static int g_init_bits;
+
+static int ClearCode;
+static int EOFCode;
+
+
+/********************************************************/
+static void compress(init_bits, data, len)
+         int   init_bits;
+         unsigned char **data;
+         int   len;
+{
+  register long fcode;
+  register int i = 0;
+  register int c;
+  register int ent;
+  register int disp;
+  register int hsize_reg;
+  register int hshift;
+  int      code_count = 0;
+
+  /* Use the output buffer as temporary storage for GIF data */
+
+  if (OutBufSize < HSIZE*(sizeof(count_int) + sizeof(unsigned short))) {
+    gif_error =1;
+    return;
+  }
+
+  /* Put htab and codetab arrays into the output buffer */
+
+  GIFOutBufSize = OutBufSize - HSIZE*(sizeof(count_int) + sizeof(unsigned short)) - 16;
+  GIFOutBufSize = GIFOutBufSize & (~0x3);       /* Make sure it's double word alligned */
+
+  htab = (count_int *) ((char *) OutBuffer + GIFOutBufSize);
+  codetab = (unsigned short *) ((char *) OutBuffer + GIFOutBufSize + HSIZE*sizeof(count_int));
+
+  /*
+   * Set up the globals:  g_init_bits - initial number of bits
+   *                      g_outfile   - pointer to output file
+   */
+  g_init_bits = init_bits;
+
+
+  /* initialize 'compress' globals */
+  maxbits = GP_BITS;
+  maxmaxcode = 1<<GP_BITS;
+  memset(htab,0,sizeof(htab));
+  memset(codetab,0,sizeof(codetab));
+  hsize = HSIZE;
+  free_ent = 0;
+  clear_flg = 0;
+  out_count = 0;
+  cur_accum = 0;
+  cur_bits = 0;
+
+  /*
+   * Set up the necessary values
+   */
+  out_count = 0;
+  clear_flg = 0;
+  maxcode = MAXCODE(n_bits = g_init_bits);
+
+  ClearCode = (1 << (init_bits - 1));
+  EOFCode = ClearCode + 1;
+  free_ent = ClearCode + 2;
+
+  char_init();
+  ent = pc2nc[data[cury][curx]];
+  curx++;
+  if (curx >= GIF_frame->width) {
+    curx = 0;
+    cury--;
+  }
+  len--;
+
+  hshift = 0;
+  for ( fcode = (long) hsize;  fcode < 65536L; fcode *= 2L )
+    hshift++;
+  hshift = 8 - hshift;                /* set hash code range bound */
+
+  hsize_reg = hsize;
+  cl_hash( (count_int) hsize_reg);            /* clear hash table */
+
+  output_GIF(ClearCode);
+  while (len) {
+    c = pc2nc[data[cury][curx]];
+    curx++;
+    if (curx >= GIF_frame->width) {
+      curx = 0;
+      cury--;
+    }
+    len--;
+
+    fcode = (long) ( ( (long) c << maxbits) + ent);
+       i = (((int) c << hshift) ^ ent);    /* xor hashing */
+
+    if ( HashTabOf (i) == fcode ) {
+      ent = CodeTabOf (i);
+      continue;
+    }
+
+    if ( (long)HashTabOf (i) < 0 )      /* empty slot */
+      goto nomatch;
+
+    disp = hsize_reg - i;           /* secondary hash (after G. Knott) */
+    if ( i == 0 )
+      disp = 1;
+
+probe:
+    if ( (i -= disp) < 0 )
+      i += hsize_reg;
+
+    if ( HashTabOf (i) == fcode ) {
+      ent = CodeTabOf (i);
+      continue;
+    }
+
+    if ( (long)HashTabOf (i) >= 0 ) 
+      goto probe;
+
+nomatch:
+    output_GIF(ent);
+    out_count++;
+    ent = c;
+
+    if ( free_ent < maxmaxcode ) {
+      CodeTabOf (i) = free_ent++; /* code -> hashtable */
+      HashTabOf (i) = fcode;
+    }
+    else
+      cl_block();
+
+  }
+  /* Put out the final code */
+  output_GIF(ent);
+  output_GIF(EOFCode);
+}
+
+
+/*****************************************************************
+ * TAG( output_GIF )
+ *
+ * Output the given code.
+ * Inputs:
+ *      code:   A n_bits-bit integer.  If == -1, then EOF.  This assumes
+ *              that n_bits =< (long)wordsize - 1.
+ * Outputs:
+ *      Outputs code to the file.
+ * Assumptions:
+ *      Chars are 8 bits long.
+ * Algorithm:
+ *      Maintain a BITS character long buffer (so that 8 codes will
+ * fit in it exactly).  Use the VAX insv instruction to insert each
+ * code in turn.  When the buffer fills up empty it and start over.
+ */
+
+static
+unsigned long masks[] = { 0x0000, 0x0001, 0x0003, 0x0007, 0x000F,
+                                  0x001F, 0x003F, 0x007F, 0x00FF,
+                                  0x01FF, 0x03FF, 0x07FF, 0x0FFF,
+                                  0x1FFF, 0x3FFF, 0x7FFF, 0xFFFF };
+
+static void output_GIF(code)
+int code;
+{
+  cur_accum &= masks[cur_bits];
+
+  if (cur_bits > 0)
+    cur_accum |= ((long)code << cur_bits);
+  else
+    cur_accum = code;
+	
+  cur_bits += n_bits;
+
+  while( cur_bits >= 8 ) {
+    char_out( (int) (cur_accum & 0xff) );
+    cur_accum >>= 8;
+    cur_bits -= 8;
+  }
+
+  /*
+   * If the next entry is going to be too big for the code size,
+   * then increase it, if possible.
+   */
+
+  if (free_ent > maxcode || clear_flg) {
+
+    if( clear_flg ) {
+      maxcode = MAXCODE (n_bits = g_init_bits);
+      clear_flg = 0;
+    }
+    else {
+      n_bits++;
+      if ( n_bits == maxbits )
+	maxcode = maxmaxcode;
+      else
+	maxcode = MAXCODE(n_bits);
+    }
+  }
+	
+  if( code == EOFCode ) {
+    /* At EOF, write the rest of the buffer */
+    while( cur_bits > 0 ) {
+      char_out( (int)(cur_accum & 0xff) );
+      cur_accum >>= 8;
+      cur_bits -= 8;
+    }
+
+    flush_char();
+  }
+}
+
+
+/********************************/
+static void cl_block ()             /* table clear for block compress */
+{
+  /* Clear out the hash table */
+
+  cl_hash ( (count_int) hsize );
+  free_ent = ClearCode + 2;
+  clear_flg = 1;
+
+  output_GIF(ClearCode);
+}
+
+
+/********************************/
+static void cl_hash(hsize)          /* reset code table */
+register count_int hsize;
+{
+  register count_int *htab_p = htab+hsize;
+  register long i;
+  register long m1 = -1;
+
+  i = hsize - 16;
+  do {                            /* might use Sys V memset(3) here */
+    *(htab_p-16) = m1;
+    *(htab_p-15) = m1;
+    *(htab_p-14) = m1;
+    *(htab_p-13) = m1;
+    *(htab_p-12) = m1;
+    *(htab_p-11) = m1;
+    *(htab_p-10) = m1;
+    *(htab_p-9) = m1;
+    *(htab_p-8) = m1;
+    *(htab_p-7) = m1;
+    *(htab_p-6) = m1;
+    *(htab_p-5) = m1;
+    *(htab_p-4) = m1;
+    *(htab_p-3) = m1;
+    *(htab_p-2) = m1;
+    *(htab_p-1) = m1;
+    htab_p -= 16;
+  } while ((i -= 16) >= 0);
+
+  for ( i += 16; i > 0; i-- )
+    *--htab_p = m1;
+}
+
+
+/******************************************************************************
+ *
+ * GIF Specific routines
+ *
+ ******************************************************************************/
+
+/*
+ * Number of characters so far in this 'packet'
+ */
+static int a_count;
+
+/*
+ * Set up the 'byte output' routine
+ */
+static void char_init()
+{
+	a_count = 0;
+}
+
+/*
+ * Define the storage for the packet accumulator
+ */
+static char accum[ 256 ];
+
+/*
+ * Add a character to the end of the current packet, and if it is 254
+ * characters, flush the packet to disk.
+ */
+static void char_out(c)
+int c;
+{
+  accum[ a_count++ ] = c;
+  if( a_count >= 254 ) 
+    flush_char();
+}
+
+/*
+ * Flush the packet to disk, and reset the accumulator
+ */
+static void flush_char()
+{
+  if (gif_error) return;
+  if( a_count > 0 ) {
+    *(op++) = a_count;
+    memcpy(op,accum,a_count);
+    op+=a_count;
+    a_count = 0;
+    
+    if (op > (unsigned char *) ((char *) OutBuffer + (GIFOutBufSize - 2048))) {
+      gif_error = 1;
+    }
+  } 
+}	
+
+
+/* ----------------------------------------------------------------------
+   int FrameBuffer_writeGIF(char *filename)
+   
+   Write a GIF file to filename
+   ----------------------------------------------------------------------- */
+
+int FrameBuffer_writeGIF(FrameBuffer *f, ColorMap *c, char *filename) {
+
+  FILE *file;
+  void *buffer;
+  int nbytes;
+  int bufsize;
+
+  file = fopen(filename,"wb");
+  if (file == NULL) return -1;
+
+  bufsize = (f->width*f->height*3)/2;
+  buffer = (void *) malloc(bufsize);
+  nbytes = FrameBuffer_makeGIF(f,c,buffer,bufsize);
+  if (nbytes == -1) { 
+    free(buffer);
+    fclose(file);
+    return -1;
+  }
+  if (fwrite(buffer,nbytes,1,file) != 1) {
+    free(buffer);
+    fclose(file);
+    return -1;
+  }
+  fclose(file);
+  free(buffer);
+  return 0;
+}
+
+    
+  
+  
+  
diff --git a/Examples/GIFPlot/Lib/matrix.c b/Examples/GIFPlot/Lib/matrix.c
new file mode 100644
index 0000000..ef0cf3a
--- /dev/null
+++ b/Examples/GIFPlot/Lib/matrix.c
@@ -0,0 +1,343 @@
+/* ----------------------------------------------------------------------------- 
+ * matrix.c
+ *
+ *     Some 4x4 matrix operations
+ * 
+ * Author(s) : David Beazley (beazley@cs.uchicago.edu)
+ * Copyright (C) 1995-1996
+ *
+ * See the file LICENSE for information on usage and redistribution.	
+ * ----------------------------------------------------------------------------- */
+
+#define MATRIX
+#include "gifplot.h"
+#include <math.h>
+
+/* ------------------------------------------------------------------------
+   Matrix new_Matrix()
+
+   Create a new 4x4 matrix.
+   ------------------------------------------------------------------------ */
+Matrix
+new_Matrix() {
+  Matrix m;
+  m = (Matrix) malloc(16*sizeof(double));
+  return m;
+}
+
+/* ------------------------------------------------------------------------
+   delete_Matrix(Matrix *m);
+
+   Destroy a matrix
+   ------------------------------------------------------------------------ */
+
+void
+delete_Matrix(Matrix m) {
+  if (m)
+    free((char *) m);
+}
+
+/* ------------------------------------------------------------------------
+   Matrix Matrix_copy(Matrix a)
+
+   Makes a copy of matrix a and returns it.
+   ------------------------------------------------------------------------ */
+
+Matrix Matrix_copy(Matrix a) {
+  int i;
+  Matrix r = 0;
+  if (a) {
+    r = new_Matrix();
+    if (r) {
+      for (i = 0; i < 16; i++)
+	r[i] = a[i];
+    }
+  }
+  return r;
+}
+    
+/* ------------------------------------------------------------------------
+   Matrix_multiply(Matrix a, Matrix b, Matrix c)
+
+   Multiplies a*b = c
+   c may be one of the source matrices
+   ------------------------------------------------------------------------ */
+void
+Matrix_multiply(Matrix a, Matrix b, Matrix c) {
+  double temp[16];
+  int    i,j,k;
+
+  for (i =0; i < 4; i++)
+    for (j = 0; j < 4; j++) {
+      temp[i*4+j] = 0.0;
+      for (k = 0; k < 4; k++) 
+	temp[i*4+j] += a[i*4+k]*b[k*4+j];
+    }
+  for (i = 0; i < 16; i++)
+    c[i] = temp[i];
+}
+  
+/* ------------------------------------------------------------------------
+   Matrix_identity(Matrix a)
+
+   Puts an identity matrix in matrix a
+   ------------------------------------------------------------------------ */
+
+void
+Matrix_identity(Matrix a) {
+  int i;
+  for (i = 0; i < 16; i++) a[i] = 0;
+  a[0] = 1;
+  a[5] = 1;
+  a[10] = 1;
+  a[15] = 1;
+}
+
+/* ------------------------------------------------------------------------
+   Matrix_zero(Matrix a)
+   
+   Puts a zero matrix in matrix a
+   ------------------------------------------------------------------------ */
+void
+Matrix_zero(Matrix a) {
+  int i;
+  for (i = 0; i < 16; i++) a[i] = 0;
+}
+
+/* ------------------------------------------------------------------------
+   Matrix_transpose(Matrix a, Matrix result)
+
+   Transposes matrix a and puts it in result.
+   ------------------------------------------------------------------------ */
+void
+Matrix_transpose(Matrix a, Matrix result) {
+  double temp[16];
+  int i,j;
+
+  for (i = 0; i < 4; i++)
+    for (j = 0; j < 4; j++) 
+      temp[4*i+j] = a[4*j+i];
+
+  for (i = 0; i < 16; i++)
+    result[i] = temp[i];
+}
+
+
+/* ------------------------------------------------------------------------
+   Matrix_gauss(Matrix a, Matrix b)
+
+   Solves ax=b for x, using Gaussian elimination. Destroys a.
+   Really only used for calculating inverses of 4x4 transformation
+   matrices.
+   ------------------------------------------------------------------------ */
+
+void Matrix_gauss(Matrix a, Matrix b) {
+  int ipiv[4], indxr[4], indxc[4];
+  int i,j,k,l,ll;
+  int irow=0, icol=0;
+  double big, pivinv;
+  double dum;
+  for (j = 0; j < 4; j++)
+    ipiv[j] = 0;
+  for (i = 0; i < 4; i++) {
+    big = 0;
+    for (j = 0; j < 4; j++) {
+      if (ipiv[j] != 1) {
+	for (k = 0; k < 4; k++) {
+	  if (ipiv[k] == 0) {
+	    if (fabs(a[4*j+k]) >= big) {
+	      big = fabs(a[4*j+k]);
+	      irow = j;
+	      icol = k;
+	    }
+	  } else if (ipiv[k] > 1)
+	    return;  /* Singular matrix */
+	}
+      }
+    }
+    ipiv[icol] = ipiv[icol]+1;
+    if (irow != icol) {
+      for (l = 0; l < 4; l++) {
+	dum = a[4*irow+l];
+	a[4*irow+l] = a[4*icol+l];
+	a[4*icol+l] = dum;
+      }
+      for (l = 0; l < 4; l++) {
+	dum = b[4*irow+l];
+	b[4*irow+l] = b[4*icol+l];
+	b[4*icol+l] = dum;
+      }
+    }
+    indxr[i] = irow;
+    indxc[i] = icol;
+    if (a[4*icol+icol] == 0) return;
+    pivinv = 1.0/a[4*icol+icol];
+    a[4*icol+icol] = 1.0;
+    for (l = 0; l < 4; l++)
+      a[4*icol+l] = a[4*icol+l]*pivinv;
+    for (l = 0; l < 4; l++)
+      b[4*icol+l] = b[4*icol+l]*pivinv;
+    for (ll = 0; ll < 4; ll++) {
+      if (ll != icol) {
+	dum = a[4*ll+icol];
+	a[4*ll+icol] = 0;
+	for (l = 0; l < 4; l++)
+	  a[4*ll+l] = a[4*ll+l] - a[4*icol+l]*dum;
+	for (l = 0; l < 4; l++)
+	  b[4*ll+l] = b[4*ll+l] - b[4*icol+l]*dum;
+      }
+    }
+  }
+  for (l = 3; l >= 0; l--) {
+    if (indxr[l] != indxc[l]) {
+      for (k = 0; k < 4; k++) {
+	dum = a[4*k+indxr[l]];
+	a[4*k+indxr[l]] = a[4*k+indxc[l]];
+	a[4*k+indxc[l]] = dum;
+      }
+    }
+  }
+}
+
+/* ------------------------------------------------------------------------
+   Matrix_invert(Matrix a, Matrix inva)
+
+   Inverts Matrix a and places the result in inva.
+   Relies on the Gaussian Elimination code above. (See Numerical recipes).
+   ------------------------------------------------------------------------ */
+void
+Matrix_invert(Matrix a, Matrix inva) {
+
+  double  temp[16];
+  int     i;
+
+  for (i = 0; i < 16; i++)
+    temp[i] = a[i];
+  Matrix_identity(inva);
+  Matrix_gauss(temp,inva);
+}
+    
+/* ------------------------------------------------------------------------
+   Matrix_transform(Matrix a, GL_Vector *r, GL_Vector *t)
+
+   Transform a vector.    a*r ----> t
+   ------------------------------------------------------------------------ */
+
+void   Matrix_transform(Matrix a, GL_Vector *r, GL_Vector *t) {
+
+  double  rx, ry, rz, rw;
+
+  rx = r->x;
+  ry = r->y;
+  rz = r->z;
+  rw = r->w;
+  t->x = a[0]*rx + a[1]*ry + a[2]*rz + a[3]*rw;
+  t->y = a[4]*rx + a[5]*ry + a[6]*rz + a[7]*rw;
+  t->z = a[8]*rx + a[9]*ry + a[10]*rz + a[11]*rw;
+  t->w = a[12]*rx + a[13]*ry + a[14]*rz + a[15]*rw;
+}
+
+/* ------------------------------------------------------------------------
+   Matrix_transform4(Matrix a, double x, double y, double z, double w, GL_Vector *t)
+
+   Transform a vector from a point specified as 4 doubles
+   ------------------------------------------------------------------------ */
+
+void   Matrix_transform4(Matrix a, double rx, double ry, double rz, double rw,
+			 GL_Vector *t) {
+
+  t->x = a[0]*rx + a[1]*ry + a[2]*rz + a[3]*rw;
+  t->y = a[4]*rx + a[5]*ry + a[6]*rz + a[7]*rw;
+  t->z = a[8]*rx + a[9]*ry + a[10]*rz + a[11]*rw;
+  t->w = a[12]*rx + a[13]*ry + a[14]*rz + a[15]*rw;
+}
+
+/* ---------------------------------------------------------------------
+   Matrix_translate(Matrix a, double tx, double ty, double tz)
+
+   Put a translation matrix in Matrix a
+   ---------------------------------------------------------------------- */
+
+void Matrix_translate(Matrix a, double tx, double ty, double tz) {
+  Matrix_identity(a);
+  a[3] = tx;
+  a[7] = ty;
+  a[11] = tz;
+  a[15] = 1;
+}
+
+/* -----------------------------------------------------------------------
+   Matrix_rotatex(Matrix a, double deg)
+
+   Produce an x-rotation matrix for given angle in degrees.
+   ----------------------------------------------------------------------- */
+void
+Matrix_rotatex(Matrix a, double deg) {
+  double r;
+
+  r = 3.1415926*deg/180.0;
+  Matrix_zero(a);
+  a[0] = 1.0;
+  a[5] = cos(r);
+  a[6] = -sin(r);
+  a[9] = sin(r);
+  a[10] = cos(r);
+  a[15] = 1.0;
+}
+
+/* -----------------------------------------------------------------------
+   Matrix_rotatey(Matrix a, double deg)
+
+   Produce an y-rotation matrix for given angle in degrees.
+   ----------------------------------------------------------------------- */
+void
+Matrix_rotatey(Matrix a, double deg) {
+  double r;
+
+  r = 3.1415926*deg/180.0;
+  Matrix_zero(a);
+  a[0] = cos(r);
+  a[2] = sin(r);
+  a[5] = 1.0;
+  a[8] = -sin(r);
+  a[10] = cos(r);
+  a[15] = 1;
+  
+}
+/* -----------------------------------------------------------------------
+   Matrix_RotateZ(Matrix a, double deg)
+
+   Produce an z-rotation matrix for given angle in degrees.
+   ----------------------------------------------------------------------- */
+void
+Matrix_rotatez(Matrix a, double deg) {
+  double r;
+
+  r = 3.1415926*deg/180.0;
+  Matrix_zero(a);
+  a[0] = cos(r);
+  a[1] = -sin(r);
+  a[4] = sin(r);
+  a[5] = cos(r);
+  a[10] = 1.0;
+  a[15] = 1.0;
+}
+
+
+/* A debugging routine */
+
+void Matrix_set(Matrix a, int i, int j, double val) {
+  a[4*j+i] = val;
+}
+
+void Matrix_print(Matrix a) {
+  int i,j;
+  for (i = 0; i < 4; i++) {
+    for (j = 0; j < 4; j++) {
+      fprintf(stdout,"%10f ",a[4*i+j]);
+    }
+    fprintf(stdout,"\n");
+  }
+  fprintf(stdout,"\n");
+}
+
diff --git a/Examples/GIFPlot/Lib/pixmap.c b/Examples/GIFPlot/Lib/pixmap.c
new file mode 100644
index 0000000..a55cf04
--- /dev/null
+++ b/Examples/GIFPlot/Lib/pixmap.c
@@ -0,0 +1,159 @@
+/* ----------------------------------------------------------------------------- 
+ * pixmap.c
+ *
+ *     Pixel maps (i.e., bitmaps)
+ * 
+ * Author(s) : David Beazley (beazley@cs.uchicago.edu)
+ * Copyright (C) 1995-1996
+ *
+ * See the file LICENSE for information on usage and redistribution.	
+ * ----------------------------------------------------------------------------- */
+
+#define PIXMAP
+#include "gifplot.h"
+
+/* -----------------------------------------------------------------------
+   PixMap *new_PixMap(int width, int height, int centerx, int centery)
+
+   Create a new pixmap of given size
+   ----------------------------------------------------------------------- */
+PixMap *new_PixMap(int width, int height, int centerx, int centery) {
+  PixMap *pm;
+  if ((width > 0) && (height > 0)) {
+    pm = (PixMap *) malloc(sizeof(PixMap));
+    pm->width = width;
+    pm->height = height;
+    pm->centerx = centerx;
+    pm->centery = centery;
+    pm->map = (int *) malloc(height*width*sizeof(int));
+    return pm;
+  }
+  return (PixMap *) 0;
+}
+
+/* --------------------------------------------------------------------------
+   void delete_PixMap(PixMap *pm)
+
+   Destroy a pixmap
+   -------------------------------------------------------------------------- */
+
+void delete_PixMap(PixMap *pm) {
+  if (pm) {
+    free((char *) pm->map);
+    free((char *) pm);
+  }
+}
+
+/* ---------------------------------------------------------------------------
+   void PixMap_set(PixMap *pm, int x, int y, int pix)
+
+   Set a pixel in the bitmap
+   --------------------------------------------------------------------------- */
+void
+PixMap_set(PixMap *pm, int x, int y, int pix) {
+  if ((x < 0) || (x>=pm->width)) return;
+  if ((y < 0) || (y>=pm->height)) return;
+
+  pm->map[pm->width*y + x] = pix;
+}
+
+/* -----------------------------------------------------------------------------
+   void FrameBuffer_drawpixmap(FrameBuffer *f, PixMap *pm, int x, int y, int fgcolor, int bgcolor) 
+
+   Draw a pixmap onto the framebuffer.   This is somewhat optimized for speed.
+   ------------------------------------------------------------------------------ */
+
+void
+FrameBuffer_drawpixmap(FrameBuffer *f, PixMap *pm, int x, int y, int fgcolor, int bgcolor) {
+
+  int startx, starty;         /* Starting location on framebuffer */
+  int startpixx = 0, startpixy = 0;   /* Starting location in pixmap      */
+  int endx, endy;             /* Ending location on framebuffer   */
+  int i,j, px, py;
+  int c;
+  
+  startx = x - pm->centerx;
+  starty = y + pm->centery;
+  endx = startx + pm->width;
+  endy = starty - pm->height;
+
+  /* Figure out if we need to clip */
+
+  if (startx < f->xmin) {
+    startpixx = f->xmin - startx;
+    startx = f->xmin;
+  }
+  if (starty >= f->ymax) {
+    startpixy = starty - f->ymax;
+    starty = f->ymax-1;
+  }
+  if (endx >= f->xmax) {
+    endx = f->xmax-1;
+  }
+  if (endy < f->ymin) {
+    endy = f->ymin;
+  }
+  py = startpixy;
+  for (j = starty; j >= endy; j--) {
+    px = startpixx;
+    for (i = startx; i < endx; i++) {
+      c = pm->map[py*pm->width + px];
+      switch (c) {
+      case GIFPLOT_FOREGROUND:
+	f->pixels[j][i] = fgcolor;
+	break;
+      case GIFPLOT_BACKGROUND:
+	f->pixels[j][i] = bgcolor;
+	break;
+      default:
+	break;
+      }
+      px++;
+    }
+    py++;
+  }
+}
+
+/**************************************************************************
+ * Some common PixMaps  (for plotting)
+ *
+ **************************************************************************/
+
+int _SQUARE_MAP[] = {
+  0,1,1,1,1,1,1,1,
+  0,1,1,1,1,1,1,1,
+  0,1,1,1,1,1,1,1,
+  0,1,1,1,1,1,1,1,
+  0,1,1,1,1,1,1,1,
+  0,1,1,1,1,1,1,1,
+  0,1,1,1,1,1,1,1,
+  0,0,0,0,0,0,0,0 };
+
+PixMap  PixMap_SQUARE = { 8,8,4,4, _SQUARE_MAP};
+
+int _TRIANGLE_MAP[] = {
+  0,0,0,1,0,0,0,0,
+  0,0,0,1,0,0,0,0,
+  0,0,1,1,1,0,0,0,
+  0,0,1,1,1,0,0,0,
+  0,1,1,1,1,1,0,0,
+  0,1,1,1,1,1,0,0,
+  1,1,1,1,1,1,1,0,
+  0,0,0,0,0,0,0,0 };
+
+PixMap  PixMap_TRIANGLE = { 8,8,4,4,_TRIANGLE_MAP};
+
+int _CROSS_MAP[] = {
+  0,0,0,1,0,0,0,0,
+  0,0,0,1,0,0,0,0,
+  0,0,0,1,0,0,0,0,
+  1,1,1,1,1,1,1,0,
+  0,0,0,1,0,0,0,0,
+  0,0,0,1,0,0,0,0,
+  0,0,0,1,0,0,0,0,
+  0,0,0,0,0,0,0,0 };
+
+PixMap  PixMap_CROSS = { 8,8,4,4,_CROSS_MAP};
+
+  
+      
diff --git a/Examples/GIFPlot/Lib/plot2d.c b/Examples/GIFPlot/Lib/plot2d.c
new file mode 100644
index 0000000..e78107b
--- /dev/null
+++ b/Examples/GIFPlot/Lib/plot2d.c
@@ -0,0 +1,445 @@
+/* ----------------------------------------------------------------------------- 
+ * plot2d.c
+ *
+ *     2-Dimensional plotting
+ * 
+ * Author(s) : David Beazley (beazley@cs.uchicago.edu)
+ * Copyright (C) 1995-1996
+ *
+ * See the file LICENSE for information on usage and redistribution.	
+ * ----------------------------------------------------------------------------- */
+
+#define PLOT2D
+
+#include "gifplot.h"
+
+/* ------------------------------------------------------------------------
+   Plot2D *new_Plot2D(FrameBuffer *frame, xmin, ymin, xmax, ymax)
+
+   Create a new 2D plot with given minimum and maximum coordinates.
+   ------------------------------------------------------------------------ */
+Plot2D  *new_Plot2D(FrameBuffer *frame,double xmin,double ymin,double xmax,double ymax) {
+  Plot2D  *p2;
+  if (frame) {
+    if (xmax <= xmin) return (Plot2D *) 0;
+    if (ymax <= ymin) return (Plot2D *) 0;
+    p2 = (Plot2D *) malloc(sizeof(Plot2D));
+    p2->frame = frame;
+    p2->xmin = xmin;
+    p2->ymin = ymin;
+    p2->xmax = xmax;
+    p2->ymax = ymax;
+    p2->view_xmin = 0;
+    p2->view_xmax = frame->width;
+    p2->view_ymin = 0;
+    p2->view_ymax = frame->height;
+    p2->xscale = LINEAR;
+    p2->yscale = LINEAR;
+    p2->dx = (p2->view_xmax - p2->view_xmin)/(p2->xmax - p2->xmin);
+    p2->dy = (p2->view_ymax - p2->view_ymin)/(p2->ymax - p2->ymin);
+    return p2;
+  }
+  return (Plot2D *) 0;
+}
+
+/* ----------------------------------------------------------------------------
+   delete_Plot2D(Plot2D *p2)
+
+   Delete a 2D plot
+   ---------------------------------------------------------------------------- */
+void
+delete_Plot2D(Plot2D *p2) {
+  if (p2)
+    free((char *) p2);
+}
+
+/* -----------------------------------------------------------------------------
+   Plot2D *Plot2D_copy(Plot2D *p2)
+
+   Makes a copy of the Plot2D data structure.
+   ----------------------------------------------------------------------------- */
+
+Plot2D *Plot2D_copy(Plot2D *p2) {
+  Plot2D *c2;
+  if (p2) {
+    c2 = (Plot2D *) malloc(sizeof(Plot2D));
+    if (c2) {
+      c2->frame = p2->frame;
+      c2->view_xmin = p2->view_xmin;
+      c2->view_ymin = p2->view_ymin;
+      c2->view_xmax = p2->view_xmax;
+      c2->view_ymax = p2->view_ymax;
+      c2->xmin = p2->xmin;
+      c2->ymin = p2->ymin;
+      c2->xmax = p2->xmax;
+      c2->ymax = p2->ymax;
+      c2->xscale = p2->xscale;
+      c2->yscale = p2->yscale;
+      c2->dx = p2->dx;
+      c2->dy = p2->dy;
+    }
+    return c2;
+  } else {
+    return (Plot2D *) 0;
+  }
+}
+
+/* -----------------------------------------------------------------------------
+   Plot2D_clear(Plot2D *p2, Pixel c)
+
+   Clear the region assigned to this plot to the given color.
+   -------------------------------------------------------------------------- */
+
+void Plot2D_clear(Plot2D *p2, Pixel c) {
+  int i,j;
+  for (i = p2->view_xmin; i < p2->view_xmax; i++)
+    for (j = p2->view_ymin; j < p2->view_ymax; j++) {
+      p2->frame->pixels[j][i] = c;
+    }
+}
+
+/* ------------------------------------------------------------------------------
+   Plot2D_setview
+
+   Sets the plot region on the framebuffer
+   ------------------------------------------------------------------------------ */
+
+void
+Plot2D_setview(Plot2D *p2, int vxmin, int vymin, int vxmax, int vymax) {
+  if (p2) {
+    p2->view_xmin = vxmin;
+    p2->view_ymin = vymin;
+    p2->view_xmax = vxmax;
+    p2->view_ymax = vymax;
+    p2->dx = (p2->view_xmax - p2->view_xmin)/(p2->xmax - p2->xmin);
+    p2->dy = (p2->view_ymax - p2->view_ymin)/(p2->ymax - p2->ymin);
+    FrameBuffer_setclip(p2->frame,vxmin,vymin,vxmax,vymax);
+  }
+}
+
+/* -------------------------------------------------------------------------------
+   Plot2D_setrange(Plot2D *p2, double xmin, double ymin, double xmax, double ymax)
+
+   Sets the plotting range.
+   ------------------------------------------------------------------------------- */
+
+void
+Plot2D_setrange(Plot2D *p2, double xmin, double ymin, double xmax, double ymax) {
+  if (p2) {
+    p2->xmin = xmin;
+    p2->ymin = ymin;
+    p2->xmax = xmax;
+    p2->ymax = ymax;
+    p2->dx = (p2->view_xmax - p2->view_xmin)/(p2->xmax - p2->xmin);
+    p2->dy = (p2->view_ymax - p2->view_ymin)/(p2->ymax - p2->ymin);
+  }
+}
+
+/* -------------------------------------------------------------------------------
+   Plot2D_setscale(Plot2D *p2, int xscale, int yscale)
+
+   Sets the plotting scaling method
+   ------------------------------------------------------------------------------- */
+
+void
+Plot2D_setscale(Plot2D *p2, int xscale, int yscale) {
+  if (p2) {
+    p2->xscale = xscale;
+    p2->yscale = yscale;
+  }
+}
+
+/* ----------------------------------------------------------------------------
+   Plot2D_transform(Plot2D *p2, double x, double y, int *px, int *py)
+   
+   Transforms x,y into screen coordinates px and py.  Result is returned
+   in px and py.  Rounds to the nearest pixel instead of truncating.
+   ----------------------------------------------------------------------------- */
+
+void
+Plot2D_transform(Plot2D *p2, double x, double y, int *px, int *py) {
+  if (p2) {
+    *px = p2->view_xmin + (int) (p2->dx*(x-p2->xmin) + 0.5);
+    *py = p2->view_ymin + (int) (p2->dy*(y-p2->ymin) + 0.5);
+  }
+}
+
+/* -------------------------------------------------------------------------------
+   Plot2D_plot(Plot2D *p2, double x, double y, Pixel color)
+
+   Plot a 2D Point of a given color
+   ------------------------------------------------------------------------------- */
+void
+Plot2D_plot(Plot2D *p2, double x, double y, Pixel color) {
+  int px, py;
+
+  Plot2D_transform(p2,x,y,&px,&py);
+  FrameBuffer_plot(p2->frame, px, py, color);
+}
+
+/* -------------------------------------------------------------------------------
+   Plot2D_box(Plot2D *p2, double x1, double y1, double x2, double y2, Pixel Color)
+
+   Plot an outline box on the 2D plot
+   ------------------------------------------------------------------------------- */
+void
+Plot2D_box(Plot2D *p2, double x1, double y1,double x2, double y2, Pixel color) {
+  int ix1, ix2,iy1, iy2;
+
+  Plot2D_transform(p2,x1,y1,&ix1,&iy1);
+  Plot2D_transform(p2,x2,y2,&ix2,&iy2);
+  FrameBuffer_box(p2->frame,ix1,iy1,ix2,iy2,color);
+}
+
+/* -------------------------------------------------------------------------------
+   Plot2D_solidbox(Plot2D *p2, double x1, double y1, double x2, double y2, Pixel Color)
+
+   Plot a solid box box on the 2D plot
+   ------------------------------------------------------------------------------- */
+void
+Plot2D_solidbox(Plot2D *p2, double x1, double y1,double x2, double y2, Pixel color) {
+  int ix1, ix2,iy1, iy2;
+
+  Plot2D_transform(p2,x1,y1,&ix1,&iy1);
+  Plot2D_transform(p2,x2,y2,&ix2,&iy2);
+  FrameBuffer_solidbox(p2->frame,ix1,iy1,ix2,iy2,color);
+}
+
+/* -------------------------------------------------------------------------------
+   Plot2D_interpbox(Plot2D *p2, double x1, double y1, double x2, double y2,
+                    Pixel c1, Pixel c2, Pixel c3, Pixel c4)
+
+   Plot a color-interpolated box on the 2D plot
+   ------------------------------------------------------------------------------- */
+void
+Plot2D_interpbox(Plot2D *p2, double x1, double y1,double x2, double y2,
+		 Pixel c1, Pixel c2, Pixel c3, Pixel c4) {
+  int ix1, ix2,iy1, iy2;
+
+  Plot2D_transform(p2,x1,y1,&ix1,&iy1);
+  Plot2D_transform(p2,x2,y2,&ix2,&iy2);
+  FrameBuffer_interpbox(p2->frame,ix1,iy1,ix2,iy2,c1,c2,c3,c4);
+}
+
+/* -------------------------------------------------------------------------------
+   Plot2D_circle(Plot2D *p2, double x, double y, double radius, Pixel color)
+
+   Make an outline circle on the 2D plot.
+   ------------------------------------------------------------------------------- */
+void
+Plot2D_circle(Plot2D *p2, double x, double y, double radius, Pixel color) {
+  int ix, iy, ir;
+
+  Plot2D_transform(p2,x,y,&ix,&iy);
+  ir = p2->dx * radius;        /* This is really incorrect. Will need ellipse */
+  if (ir > 1) 
+    FrameBuffer_circle(p2->frame,ix,iy,ir,color);
+  else
+    FrameBuffer_plot(p2->frame,ix,iy,color);
+  
+}
+  
+/* -------------------------------------------------------------------------------
+   Plot2D_solidcircle(Plot2D *p2, double x, double y, double radius, Pixel color)
+
+   Make an solid circle on the 2D plot.
+   ------------------------------------------------------------------------------- */
+void
+Plot2D_solidcircle(Plot2D *p2, double x, double y, double radius, Pixel color) {
+  int ix, iy, ir;
+
+  Plot2D_transform(p2,x,y,&ix,&iy);
+  ir = p2->dx * radius;        /* This is really incorrect. Will need ellipse */
+  if (ir > 1) 
+    FrameBuffer_solidcircle(p2->frame,ix,iy,ir,color);
+  else
+    FrameBuffer_plot(p2->frame,ix,iy,color);
+}
+
+/* -------------------------------------------------------------------------------
+   Plot2D_line(Plot2D *p2, double x1, double y1, double x2, double y2, Pixel color)
+
+   Draw a line
+   ------------------------------------------------------------------------------- */
+
+void
+Plot2D_line(Plot2D *p2, double x1, double y1, double x2, double y2, Pixel color) {
+  int ix1, ix2, iy1, iy2;
+
+  Plot2D_transform(p2,x1,y1,&ix1,&iy1);
+  Plot2D_transform(p2,x2,y2,&ix2,&iy2);
+  FrameBuffer_line(p2->frame,ix1,iy1,ix2,iy2,color);
+}
+
+
+
+/* -------------------------------------------------------------------------------
+   Plot2D_start(Plot2D *p2)
+
+   This should be called before starting to make a 2D plot. It will change
+   the viewport coordinates for the framebuffer and do other stuff.
+   ------------------------------------------------------------------------------- */
+
+void Plot2D_start(Plot2D *p2) {
+  if (p2) {
+    FrameBuffer_setclip(p2->frame, p2->view_xmin,p2->view_ymin,p2->view_xmax, p2->view_ymax);
+    p2->dx = (p2->view_xmax - p2->view_xmin)/(p2->xmax - p2->xmin);
+    p2->dy = (p2->view_ymax - p2->view_ymin)/(p2->ymax - p2->ymin);
+  }
+}
+
+/* --------------------------------------------------------------------------
+   void Plot2D_drawpixmap(Plot2D *p2, PixMap *pm, double x, double y, Pixel color, Pixel bgcolor)
+
+   Draw a pixel map at the given coordinates.  (Used for putting symbols on 2D
+   plots).
+   -------------------------------------------------------------------------- */
+void
+Plot2D_drawpixmap(Plot2D *p2, PixMap *pm, double x, double y, Pixel color, Pixel bgcolor) {
+  int ix, iy;
+
+  Plot2D_transform(p2,x,y,&ix,&iy);
+  FrameBuffer_drawpixmap(p2->frame,pm,ix,iy,color,bgcolor);
+}
+
+/* ----------------------------------------------------------------------------
+   void Plot2D_xaxis(Plot2D *p2, double x, double y, double xtick, int ticklength, Pixel color)
+
+   Draw an X axis bar at location x,y with ticks spaced every xtick units.
+   Ticks are spaced starting at "x"
+   ----------------------------------------------------------------------------- */
+
+void Plot2D_xaxis(Plot2D *p2, double x, double y, double xtick, int ticklength, Pixel color) {
+  int ix, iy,iy2;
+  double xt;
+
+  /* Draw a line fox the axis */
+  
+  Plot2D_line(p2,p2->xmin,y,p2->xmax,y,color);
+  xt = x;
+  while (xt >= p2->xmin) {
+    Plot2D_transform(p2,xt,y,&ix,&iy);
+    iy2 = iy+ticklength;
+    iy = iy-ticklength;
+    FrameBuffer_line(p2->frame,ix,iy,ix,iy2,color);
+    xt = xt - xtick;
+  }
+  xt = x + xtick;
+  while (xt < p2->xmax) {
+    Plot2D_transform(p2,xt,y,&ix,&iy);
+    iy2 = iy+ticklength;
+    iy = iy-ticklength;
+    FrameBuffer_line(p2->frame,ix,iy,ix,iy2,color);
+    xt = xt + xtick;
+  }
+}
+
+
+/* ----------------------------------------------------------------------------
+   void Plot2D_yaxis(Plot2D *p2, double x, double y, double ytick, int ticklength, Pixel c)
+
+   Draw an Y axis bar at location x,y with ticks spaced every xtick units.
+   Ticks are spaced starting at "y"
+   ----------------------------------------------------------------------------- */
+
+void Plot2D_yaxis(Plot2D *p2, double x, double y, double ytick, int ticklength, Pixel color) {
+  int ix, iy, ix2;
+  double yt;
+
+  /* Draw a line fox the axis */
+  
+  Plot2D_line(p2,x,p2->ymin,x,p2->ymax,color);
+  yt = y;
+  while (yt >= p2->ymin) {
+    Plot2D_transform(p2,x,yt,&ix,&iy);
+    ix2 = ix+ticklength;
+    ix = ix-ticklength;
+    FrameBuffer_line(p2->frame,ix,iy,ix2,iy,color);
+    yt = yt - ytick;
+  }
+  yt = y + ytick;
+  while (yt < p2->ymax) {
+    Plot2D_transform(p2,x,yt,&ix,&iy);
+    ix2 = ix+ticklength;
+    ix = ix-ticklength;
+    FrameBuffer_line(p2->frame,ix,iy,ix2,iy,color);
+    yt = yt + ytick;
+  }
+}
+
+
+/* -------------------------------------------------------------------------
+   Plot2D_triangle(Plot2D *p2, double x1, double y1, 
+                               double x2, double y2, 
+			       double x3, double y3, 
+			       Pixel fillcolor)
+
+   This function draws a 2D outline triangle.
+   -------------------------------------------------------------------------- */
+
+void Plot2D_triangle(Plot2D *p2, double x1, double y1,
+		     double x2, double y2, 
+		     double x3, double y3, Pixel color) {
+
+  Plot2D_line(p2,x1,y1,x2,y2,color);
+  Plot2D_line(p2,x2,y2,x3,y3,color);
+  Plot2D_line(p2,x3,y3,x1,y1,color);
+
+}
+
+
+/* -------------------------------------------------------------------------
+   Plot2D_solidtriangle(Plot2D *p2, double x1, double y1, 
+                               double x2, double y2, 
+			       double x3, double y3, 
+			       Pixel color)
+
+   This function draws a 2D filled triangle.    Can be used to
+   draw other primitives such as quadralaterals, etc...
+
+   -------------------------------------------------------------------------- */
+
+void Plot2D_solidtriangle(Plot2D *p2, double x1, double y1, 
+
+			  double x2, double y2, 
+			  double x3, double y3, Pixel color) {
+
+  int        tx1, tx2, tx3, ty1, ty2, ty3;
+  
+  /* Transform the three points into screen coordinates */
+
+  Plot2D_transform(p2,x1,y1,&tx1,&ty1);
+  Plot2D_transform(p2,x2,y2,&tx2,&ty2);
+  Plot2D_transform(p2,x3,y3,&tx3,&ty3);
+
+  FrameBuffer_solidtriangle(p2->frame,tx1,ty1,tx2,ty2,tx3,ty3,color);
+  
+}    
+
+/* -------------------------------------------------------------------------
+   Plot2D_interptriangle(Plot2D *p2, double x1, double y1, Pixel c1,
+                               double x2, double y2, Pixel c2,
+			       double x3, double y3, Pixel c3);
+
+   This function draws a 2D filled triangle with color interpolation.
+   Can be used to draw other primitives such as quadralaterals, etc...
+   -------------------------------------------------------------------------- */
+
+void Plot2D_interptriangle(Plot2D *p2, double x1, double y1, Pixel c1,
+			   double x2, double y2, Pixel c2,
+			   double x3, double y3, Pixel c3) {
+
+  int        tx1, tx2, tx3, ty1, ty2, ty3;
+  
+  /* Transform the three points into screen coordinates */
+
+  Plot2D_transform(p2,x1,y1,&tx1,&ty1);
+  Plot2D_transform(p2,x2,y2,&tx2,&ty2);
+  Plot2D_transform(p2,x3,y3,&tx3,&ty3);
+
+  FrameBuffer_interptriangle(p2->frame,tx1,ty1,c1,tx2,ty2,c2,tx3,ty3,c3);
+  
+}    
+  
+  
+
diff --git a/Examples/GIFPlot/Lib/plot3d.c b/Examples/GIFPlot/Lib/plot3d.c
new file mode 100644
index 0000000..387e420
--- /dev/null
+++ b/Examples/GIFPlot/Lib/plot3d.c
@@ -0,0 +1,2181 @@
+/* ----------------------------------------------------------------------------- 
+ * plot3d.c
+ *
+ *     Three-dimensional plotting.
+ * 
+ * Author(s) : David Beazley (beazley@cs.uchicago.edu)
+ * Copyright (C) 1995-1996
+ *
+ * See the file LICENSE for information on usage and redistribution.	
+ * ----------------------------------------------------------------------------- */
+
+#define PLOT3D
+#include "gifplot.h"
+#include <math.h>
+#include <float.h>
+
+#define ORTHO        1
+#define PERSPECTIVE  2
+/* ------------------------------------------------------------------------
+   Plot3D *new_Plot3D(FrameBuffer *f, double xmin, double ymin, double zmin,
+   double xmax, double ymax, double zmax)
+
+   Creates a new 3D plot.  Min and max coordinates are primarily used to
+   pick some default parameters.  Returns NULL on failure
+   ------------------------------------------------------------------------- */
+
+Plot3D *new_Plot3D(FrameBuffer *f, double xmin, double ymin, double zmin,
+		   double xmax, double ymax, double zmax) {
+
+  Plot3D  *p3;
+  void Plot3D_maketransform(Plot3D *p3);
+  
+  /* Check to make sure the framebuffer and min/max parameters are valid */
+
+  if (!f) return (Plot3D *) 0;
+  if ((xmin > xmax) || (ymin > ymax) || (zmin > zmax)) return (Plot3D *) 0;
+
+  p3 = (Plot3D *) malloc(sizeof(Plot3D));
+  p3->frame = f;
+  p3->xmin = xmin;
+  p3->ymin = ymin;
+  p3->zmin = zmin;
+  p3->xmax = xmax;
+  p3->ymax = ymax;
+  p3->zmax = zmax;
+
+  /* Set view region to the entire size of the framebuffer */
+
+  p3->view_xmin = 0;
+  p3->view_ymin = 0;
+  p3->view_xmax = f->width;
+  p3->view_ymax = f->height;
+  p3->width = f->width;
+  p3->height = f->height;
+  
+  /* Calculate a center point based off the min and max coordinates given */
+
+  p3->xcenter = (xmax - xmin)/2.0 + xmin;
+  p3->ycenter = (ymax - ymin)/2.0 + ymin;
+  p3->zcenter = (zmax - zmin)/2.0 + zmin;
+
+  /* Calculate the aspect ratio of the viewing region */
+
+  p3->aspect = (double) f->width/(double) f->height;
+
+  /* Set default view parameters */
+  p3->xshift = 1.0;
+  p3->yshift = 1.0;
+  p3->zoom   = 0.5;
+  p3->fovy   = 40.0;         /* 40 degree field of view */
+
+  /* Create matrices */
+
+  p3->model_mat = new_Matrix();
+  p3->view_mat = new_Matrix();
+  p3->center_mat = new_Matrix();
+  p3->fullmodel_mat = new_Matrix();
+  p3->trans_mat = new_Matrix();
+  p3->pers_mode = ORTHO;
+
+  FrameBuffer_zresize(p3->frame,p3->width, p3->height);
+  Matrix_identity(p3->view_mat);
+  Matrix_identity(p3->model_mat);
+  Matrix_translate(p3->center_mat, -p3->xcenter,-p3->ycenter,-p3->zcenter);
+  Plot3D_maketransform(p3);
+  return p3;
+}
+
+/* ---------------------------------------------------------------------
+   delete_Plot3D(Plot3D *p3)
+
+   Deletes a 3D plot
+   --------------------------------------------------------------------- */
+
+void delete_Plot3D(Plot3D *p3) {
+  if (p3) {
+    delete_Matrix(p3->view_mat);
+    delete_Matrix(p3->model_mat);
+    delete_Matrix(p3->trans_mat);
+    free((char *) p3);
+  }
+}
+
+/* ---------------------------------------------------------------------
+   Plot3D *Plot3D_copy(Plot3D *p3)
+
+   This makes a copy of the 3D plot structure and returns a pointer to it.
+   --------------------------------------------------------------------- */
+
+Plot3D *Plot3D_copy(Plot3D *p3) {
+  Plot3D *c3;
+  if (p3) {
+    c3 = (Plot3D *) malloc(sizeof(Plot3D));
+    if (c3) {
+      c3->frame = p3->frame;
+      c3->view_xmin = p3->view_xmin;
+      c3->view_ymin = p3->view_ymin;
+      c3->view_xmax = p3->view_xmax;
+      c3->view_ymax = p3->view_ymax;
+      c3->xmin = p3->xmin;
+      c3->ymin = p3->ymin;
+      c3->zmin = p3->zmin;
+      c3->xmax = p3->xmax;
+      c3->ymax = p3->ymax;
+      c3->zmax = p3->zmax;
+      c3->xcenter = p3->xcenter;
+      c3->ycenter = p3->ycenter;
+      c3->zcenter = p3->zcenter;
+      c3->fovy = p3->fovy;
+      c3->aspect = p3->aspect;
+      c3->znear = p3->znear;
+      c3->zfar = p3->zfar;
+      c3->center_mat = Matrix_copy(p3->center_mat);
+      c3->model_mat = Matrix_copy(p3->model_mat);
+      c3->view_mat = Matrix_copy(p3->view_mat);
+      c3->fullmodel_mat = Matrix_copy(p3->fullmodel_mat);
+      c3->trans_mat = Matrix_copy(p3->trans_mat);
+      c3->lookatz = p3->lookatz;
+      c3->xshift = p3->xshift;
+      c3->yshift = p3->yshift;
+      c3->zoom = p3->zoom;
+      c3->width = p3->width;
+      c3->height = p3->height;
+      c3->pers_mode = p3->pers_mode;
+    }
+    return c3;
+  } else {
+    return (Plot3D *) 0;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   Plot3D_clear(Plot3D *p3, Pixel bgcolor)
+
+   Clear the pixel and zbuffer only for the view region of this image.
+   ---------------------------------------------------------------------- */
+void
+Plot3D_clear(Plot3D *p3, Pixel bgcolor) {
+  int i,j;
+
+  for (i = p3->view_xmin; i < p3->view_xmax; i++)
+    for (j = p3->view_ymin; j < p3->view_ymax; j++) {
+      p3->frame->pixels[j][i] = bgcolor;
+      p3->frame->zbuffer[j][i] = ZMIN;
+    }
+}
+
+/* ---------------------------------------------------------------------
+   Plot3D_maketransform(Plot3D *p3)
+
+   This function builds the total 3D transformation matrix from a
+   collection of components.
+
+   Trans = View * Rotation * Center
+
+   Where View is the viewing transformation matrix, Rotation is the
+   model rotation matrix, Center is the translation matrix used to
+   center the Model at the origin.
+   --------------------------------------------------------------------- */
+
+void
+Plot3D_maketransform(Plot3D *p3) {
+
+  Matrix_multiply(p3->model_mat,p3->center_mat, p3->fullmodel_mat);
+  Matrix_multiply(p3->view_mat,p3->fullmodel_mat, p3->trans_mat);
+}
+   
+/* ---------------------------------------------------------------------
+   Plot3D_perspective(Plot3D *p3, double fovy, double znear, double zfar)
+
+   Sets up the perspective viewing transformation.   Assumes "lookat"
+   has already been called.
+   --------------------------------------------------------------------- */
+
+void
+Plot3D_perspective(Plot3D *p3, double fovy, double znear, double zfar) {
+  double theta;
+  double mat[16];
+
+  p3->fovy = fovy;
+  p3->znear = znear;
+  p3->zfar = zfar;
+
+  theta = 3.1415926*fovy/180.0;
+
+  Matrix_identity(mat);
+  mat[0] = cos(theta/2.0)/(sin(theta/2.0)*p3->aspect);
+  mat[5] = cos(theta/2.0)/(sin(theta/2.0));
+  mat[10] = -(zfar + znear)/(zfar-znear);
+  mat[14] = -1.0;
+  mat[11] = -(2*zfar*znear)/(zfar - znear);
+  mat[15] = 0.0;
+
+  /* Update the view transformation matrix */
+  
+  Matrix_multiply(mat,p3->view_mat,p3->view_mat);
+
+  /* Update the global transformation matrix */
+  
+  Plot3D_maketransform(p3);
+  p3->pers_mode = PERSPECTIVE;
+
+}
+
+/* ---------------------------------------------------------------------
+   Plot3D_lookat(Plot3D *p3, double z)
+
+   A greatly simplified version of (lookat).  Specifies the position
+   of the viewpoint. (can be used for moving image in or out).
+
+   Destroys the current viewing transformation matrix, so it will have
+   to be recalculated.
+   --------------------------------------------------------------------- */
+  
+void
+Plot3D_lookat(Plot3D *p3, double z) {
+  if (p3) {
+    Matrix_translate(p3->view_mat, 0,0,-z);
+    p3->lookatz = z;
+    Plot3D_maketransform(p3);
+  }
+}
+
+/* -------------------------------------------------------------------------
+   Plot3D_autoperspective(Plot3D *p3, double fovy)
+
+   Automatically figures out a semi-decent viewpoint given the
+   min,max parameters currently set for this image
+   ------------------------------------------------------------------------- */
+
+void
+Plot3D_autoperspective(Plot3D *p3, double fovy) {
+
+  /* Make a perspective transformation matrix for this system */
+
+  double  zfar;
+  double  znear;
+  double  d, dmax;
+  double  cx,cy,cz;
+  double  xmin,xmax,ymin,ymax,zmin,zmax;
+
+  xmin = p3->xmin;
+  ymin = p3->ymin;
+  zmin = p3->zmin;
+  xmax = p3->xmax;
+  ymax = p3->ymax;
+  zmax = p3->zmax;
+  cx = p3->xcenter;
+  cy = p3->ycenter;
+  cz = p3->zcenter;
+  
+  /* Calculate longest point from center point */
+
+  dmax = (xmin-cx)*(xmin-cx) + (ymin-cy)*(ymin-cy) + (zmin-cz)*(zmin-cz);
+  d = (xmax-cx)*(xmax-cx) + (ymin-cy)*(ymin-cy) + (zmin-cz)*(zmin-cz);
+  if (d > dmax) dmax = d;
+  d = (xmin-cx)*(xmin-cx) + (ymax-cy)*(ymax-cy) + (zmin-cz)*(zmin-cz);
+  if (d > dmax) dmax = d;
+  d = (xmax-cx)*(xmax-cx) + (ymax-cy)*(ymax-cy) + (zmin-cz)*(zmin-cz);
+  if (d > dmax) dmax = d;
+  d = (xmin-cx)*(xmin-cx) + (ymin-cy)*(ymin-cy) + (zmax-cz)*(zmax-cz);
+  if (d > dmax) dmax = d;
+  d = (xmax-cx)*(xmax-cx) + (ymin-cy)*(ymin-cy) + (zmax-cz)*(zmax-cz);
+  if (d > dmax) dmax = d;
+  d = (xmin-cx)*(xmin-cx) + (ymax-cy)*(ymax-cy) + (zmax-cz)*(zmax-cz);
+  if (d > dmax) dmax = d;
+  d = (xmax-cx)*(xmax-cx) + (ymax-cy)*(ymax-cy) + (zmax-cz)*(zmax-cz);
+  if (d > dmax) dmax = d;
+
+  dmax = sqrt(dmax);
+  d = p3->lookatz;
+  
+  znear = d - dmax;
+  zfar = znear+1.5*dmax;
+  Plot3D_perspective(p3, fovy,znear,zfar);
+  
+}
+
+   
+/* ---------------------------------------------------------------------
+   Plot3D_ortho(Plot3D *p3, double left, double right, double bottom, double top)
+
+   Sets up an orthographic viewing transformation.   
+   --------------------------------------------------------------------- */
+
+void
+Plot3D_ortho(Plot3D *p3, double left, double right, double bottom, double top) {
+
+
+  Matrix_identity(p3->view_mat);
+  p3->view_mat[0] = (2.0/(right - left))/p3->aspect;
+  p3->view_mat[5] = 2.0/(top - bottom);
+  p3->view_mat[10] = -1;
+  p3->view_mat[15] = 1.0;
+  p3->view_mat[3] = -(right+left)/(right-left);
+  p3->view_mat[7] = -(top+bottom)/(top-bottom);
+
+  /* Update the global transformation matrix */
+  
+  Plot3D_maketransform(p3);
+  p3->pers_mode = ORTHO;
+  p3->ortho_left = left;
+  p3->ortho_right = right;
+  p3->ortho_bottom = bottom;
+  p3->ortho_top = top;
+
+}
+
+/* ---------------------------------------------------------------------
+   Plot3D_autoortho(Plot3D *p3)
+
+   Automatically pick an orthographic projection that's probably
+   pretty good.
+   --------------------------------------------------------------------- */
+
+void
+Plot3D_autoortho(Plot3D *p3) {
+
+  /* Make a perspective transformation matrix for this system */
+
+  double  d, dmax;
+  double  cx,cy,cz;
+  double  xmin,xmax,ymin,ymax,zmin,zmax;
+
+  xmin = p3->xmin;
+  ymin = p3->ymin;
+  zmin = p3->zmin;
+  xmax = p3->xmax;
+  ymax = p3->ymax;
+  zmax = p3->zmax;
+  cx = p3->xcenter;
+  cy = p3->ycenter;
+  cz = p3->zcenter;
+  
+  /* Calculate longest point from center point */
+
+  dmax = (xmin-cx)*(xmin-cx) + (ymin-cy)*(ymin-cy) + (zmin-cz)*(zmin-cz);
+  d = (xmax-cx)*(xmax-cx) + (ymin-cy)*(ymin-cy) + (zmin-cz)*(zmin-cz);
+  if (d > dmax) dmax = d;
+  d = (xmin-cx)*(xmin-cx) + (ymax-cy)*(ymax-cy) + (zmin-cz)*(zmin-cz);
+  if (d > dmax) dmax = d;
+  d = (xmax-cx)*(xmax-cx) + (ymax-cy)*(ymax-cy) + (zmin-cz)*(zmin-cz);
+  if (d > dmax) dmax = d;
+  d = (xmin-cx)*(xmin-cx) + (ymin-cy)*(ymin-cy) + (zmax-cz)*(zmax-cz);
+  if (d > dmax) dmax = d;
+  d = (xmax-cx)*(xmax-cx) + (ymin-cy)*(ymin-cy) + (zmax-cz)*(zmax-cz);
+  if (d > dmax) dmax = d;
+  d = (xmin-cx)*(xmin-cx) + (ymax-cy)*(ymax-cy) + (zmax-cz)*(zmax-cz);
+  if (d > dmax) dmax = d;
+  d = (xmax-cx)*(xmax-cx) + (ymax-cy)*(ymax-cy) + (zmax-cz)*(zmax-cz);
+  if (d > dmax) dmax = d;
+
+  dmax = sqrt(dmax);
+
+  Plot3D_ortho(p3,-dmax,dmax,-dmax,dmax);
+  
+}
+
+
+
+/* -------------------------------------------------------------------------
+   Plot3D_setview(Plot3D *p3, int vxmin, int vymin, int vxmax, int vymax)
+
+   Sets the viewport for this 3D graph.  Will recalculate all of the
+   local viewing transformation matrices accordingly.
+   ------------------------------------------------------------------------- */
+void
+Plot3D_setview(Plot3D *p3, int vxmin, int vymin, int vxmax, int vymax) {
+  if (p3) {
+    if ((vxmin > vxmax) || (vymin >vymax)) return;
+    p3->view_xmin = vxmin;
+    p3->view_ymin = vymin;
+    p3->view_xmax = vxmax;
+    p3->view_ymax = vymax;
+    p3->width = (vxmax - vxmin);
+    p3->height = (vymax - vymin);
+    p3->aspect = (double) p3->width/(double) p3->height;
+
+    /* Fix up the viewing transformation matrix */
+
+    if (p3->pers_mode == PERSPECTIVE) {
+      Plot3D_lookat(p3,p3->lookatz);
+      Plot3D_perspective(p3,p3->fovy,p3->znear,p3->zfar);
+    } else {
+      Plot3D_ortho(p3,p3->ortho_left,p3->ortho_right,p3->ortho_bottom, p3->ortho_top);
+    }
+    FrameBuffer_setclip(p3->frame,vxmin,vymin,vxmax,vymax);
+  }
+}
+
+/* ---------------------------------------------------------------------------
+   Plot2D_start(Plot2D *p3)
+
+   Set up viewing region and other parameters for this image.
+   --------------------------------------------------------------------------- */
+
+void
+Plot3D_start(Plot3D *p3) {
+  if (p3)
+    FrameBuffer_setclip(p3->frame, p3->view_xmin,p3->view_ymin,p3->view_xmax, p3->view_ymax);
+
+}
+
+/* -------------------------------------------------------------------------
+   Plot3D_plot(Plot3D *p3, double x, double y, double z, Pixel Color)
+
+   Plot a 3D point
+   ------------------------------------------------------------------------- */
+   
+void
+Plot3D_plot(Plot3D *p3, double x, double y, double z, Pixel color) { 
+
+  GL_Vector t;
+  int      ix, iy;
+  double   invw;
+  FrameBuffer *f;
+  
+  /* Perform a transformation */
+
+  Matrix_transform4(p3->trans_mat,x,y,z,1,&t);
+
+  /* Scale the coordinates into unit cube */
+  
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+#ifdef GL_DEBUG
+  fprintf(stdout,"t.x = %g,   t.y = %g,    t.z = %g\n", t.x,t.y,t.z);
+#endif
+  /* Calculate the x and y coordinates */
+
+  ix = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5); 
+  iy = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5); 
+
+  if ((ix >= 0) && (ix < p3->width) &&
+      (iy >= 0) && (ix < p3->height)) {
+    ix += p3->view_xmin;
+    iy += p3->view_ymin;
+    f = p3->frame;
+    if (t.z <= f->zbuffer[iy][ix]) {
+      f->pixels[iy][ix] = color;
+      f->zbuffer[iy][ix] = t.z;
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   Plot3D_rotx(Plot3D *p3, double deg)
+
+   Rotate the model around its x axis.
+   ---------------------------------------------------------------------- */
+
+void
+Plot3D_rotx(Plot3D *p3, double deg) {
+  double temp[16];
+
+  Matrix_rotatex(temp,deg);         /* Construct a x rotation matrix */
+  Matrix_multiply(p3->model_mat,temp,p3->model_mat);
+  Plot3D_maketransform(p3);
+
+}
+
+/* ----------------------------------------------------------------------
+   Plot3D_roty(Plot3D *p3, double deg)
+
+   Rotate the model around its y axis.
+   ---------------------------------------------------------------------- */
+
+void
+Plot3D_roty(Plot3D *p3, double deg) {
+  double temp[16];
+
+  Matrix_rotatey(temp,deg);         /* Construct a y rotation matrix */
+  Matrix_multiply(p3->model_mat,temp,p3->model_mat);
+  Plot3D_maketransform(p3);
+
+}
+
+/* ----------------------------------------------------------------------
+   Plot3D_rotz(Plot3D *p3, double deg)
+
+   Rotate the model around its z axis.
+   ---------------------------------------------------------------------- */
+
+void
+Plot3D_rotz(Plot3D *p3, double deg) {
+  double temp[16];
+
+  Matrix_rotatez(temp,deg);         /* Construct a z rotation matrix */
+  Matrix_multiply(p3->model_mat,temp,p3->model_mat);
+  Plot3D_maketransform(p3);
+
+}
+
+
+/* ----------------------------------------------------------------------
+   Plot3D_rotd(Plot3D *p3, double deg)
+
+   Rotate the model down
+   ---------------------------------------------------------------------- */
+
+void
+Plot3D_rotd(Plot3D *p3, double deg) {
+  double temp[16];
+
+  Matrix_rotatex(temp,deg);         /* Construct a x rotation matrix */
+  Matrix_multiply(temp, p3->model_mat,p3->model_mat);
+  Plot3D_maketransform(p3);
+
+}
+    
+
+/* ----------------------------------------------------------------------
+   Plot3D_rotu(Plot3D *p3, double deg)
+
+   Rotate the model up
+   ---------------------------------------------------------------------- */
+
+void
+Plot3D_rotu(Plot3D *p3, double deg) {
+  double temp[16];
+
+  Matrix_rotatex(temp,-deg);         /* Construct a x rotation matrix */
+  Matrix_multiply(temp,p3->model_mat,p3->model_mat);
+  Plot3D_maketransform(p3);
+
+}
+
+
+/* ----------------------------------------------------------------------
+   Plot3D_rotr(Plot3D *p3, double deg)
+
+   Rotate the model down
+   ---------------------------------------------------------------------- */
+
+void
+Plot3D_rotr(Plot3D *p3, double deg) {
+  double temp[16];
+
+  Matrix_rotatey(temp,deg);         /* Construct a y rotation matrix */
+  Matrix_multiply(temp, p3->model_mat,p3->model_mat);
+  Plot3D_maketransform(p3);
+
+}
+    
+
+/* ----------------------------------------------------------------------
+   Plot3D_rotl(Plot3D *p3, double deg)
+
+   Rotate the model left
+   ---------------------------------------------------------------------- */
+
+void
+Plot3D_rotl(Plot3D *p3, double deg) {
+  double temp[16];
+
+  Matrix_rotatey(temp,-deg);         /* Construct a y rotation matrix */
+  Matrix_multiply(temp,p3->model_mat,p3->model_mat);
+  Plot3D_maketransform(p3);
+
+}
+
+
+/* ----------------------------------------------------------------------
+   Plot3D_rotc(Plot3D *p3, double deg)
+
+   Rotate the model around center point
+   ---------------------------------------------------------------------- */
+
+void
+Plot3D_rotc(Plot3D *p3, double deg) {
+  double temp[16];
+
+  Matrix_rotatez(temp,-deg);         /* Construct a z rotation matrix */
+  Matrix_multiply(temp,p3->model_mat,p3->model_mat);
+  Plot3D_maketransform(p3);
+}
+
+/* -------------------------------------------------------------------------
+   Plot3D_zoom(Plot3D *p3, double percent)
+
+   Zooms in or out the current image.   percent defines a percentage of
+   zoom.
+
+   Zooming is actually done by adjusting the perspective field of view
+   instead of scaling the model or moving in the viewpoint.    This
+   seems to work the best.
+   ------------------------------------------------------------------------- */
+
+void
+Plot3D_zoom(Plot3D *p3, double percent) {
+
+  double scale;
+  double dx;
+  if (percent <= 0) return;
+  scale = percent/100.0;
+
+  dx = (1.0/scale - 1.0)/(2*p3->zoom); /* Don't even ask where this came from */
+  p3->xshift += dx;
+  p3->yshift += dx;
+  p3->zoom = p3->zoom*scale;
+
+#ifdef OLD
+  p3->fovy = p3->fovy*scale;
+  if (p3->fovy > 170.0) p3->fovy = 170.0;
+  if (p3->fovy == 0) p3->fovy = 0.0001;
+  Plot3D_lookat(p3,p3->lookatz);
+  Plot3D_perspective(p3,p3->fovy,p3->znear,p3->zfar);
+#endif
+}
+
+/* --------------------------------------------------------------------------
+   Plot3D_left(Plot3D *p3, double s)
+
+   Shifts the image to the left by s units.   This is a little funky.
+
+   s is scaled so that s = 100 equals one full screen.
+   -------------------------------------------------------------------------- */
+void
+Plot3D_left(Plot3D *p3, double s) {
+  p3->xshift -= (s/100.0)/p3->zoom;
+}
+
+/* --------------------------------------------------------------------------
+   Plot3D_right(Plot3D *p3, double s)
+
+   Shifts the image to the right by s units. 
+
+   s is scaled so that s = 100 equals one full screen.
+   -------------------------------------------------------------------------- */
+void
+Plot3D_right(Plot3D *p3, double s) {
+  p3->xshift += (s/100.0)/p3->zoom;
+}
+
+/* --------------------------------------------------------------------------
+   Plot3D_up(Plot3D *p3, double s)
+
+   Shifts the image up left by s units.
+
+   s is scaled so that s = 100 equals one full screen.
+   -------------------------------------------------------------------------- */
+void
+Plot3D_up(Plot3D *p3, double s) {
+  p3->yshift += (s/100.0)/p3->zoom;
+}
+
+/* --------------------------------------------------------------------------
+   Plot3D_down(Plot3D *p3, double s)
+
+   Shifts the image down by s units.  
+
+   s is scaled so that s = 100 equals one full screen.
+   -------------------------------------------------------------------------- */
+void
+Plot3D_down(Plot3D *p3, double s) {
+  p3->yshift -= (s/100.0)/p3->zoom;
+}
+
+/* -------------------------------------------------------------------------
+   Plot3D_center(Plot3D *p3, double cx, double cy)
+
+   Centers the image on a point in the range (0,0) - (100,100)
+   ------------------------------------------------------------------------- */
+void
+Plot3D_center(Plot3D *p3, double cx, double cy) {
+  Plot3D_left(p3,cx-50);
+  Plot3D_down(p3,cy-50);
+}
+
+
+
+/***************************************************************************
+ *                       3d Primitives                                     *
+ ***************************************************************************/
+
+/* -------------------------------------------------------------------------
+   Plot3D_horizontal(Plot3D *p3, int xmin, int xmax, int y, double z1, double z2, Pixel color)
+
+   Draws a "Horizontal" line on the framebuffer between two screen coordinates,
+   but also supplies z-values and zbuffering.   This function probably isn't
+   too useful by itself, but will be used by a number of other primitives.
+   -------------------------------------------------------------------------- */
+
+void Plot3D_horizontal(Plot3D *p3, int xmin, int xmax, int y, Zvalue z1, Zvalue z2, Pixel color) {
+  Pixel  *p;
+  FrameBuffer *f;
+  int     i;
+  Zvalue  *zbuf,z,mz;
+  int     startx, endx;
+
+  f = p3->frame;
+  if ((y < f->ymin) || (y >= f->ymax)) return;
+  if (xmin > f->xmax) return;
+  if (xmin < f->xmin) startx = f->xmin;
+  else startx = xmin;
+  if (xmax < f->xmin) return;
+  if (xmax >= f->xmax) endx = f->xmax - 1;
+  else endx = xmax;
+
+  /* Calculate z slope */
+
+  if (xmax != xmin) {
+    mz = (Zvalue) ((double) (z2 - z1)/(double) (xmax - xmin));
+  } else {
+    mz = 0;
+  }
+
+  /* Draw it */
+  
+  p = &f->pixels[y][startx];
+  zbuf = &f->zbuffer[y][startx];
+  z = (Zvalue) (mz*(startx-xmin) + z1);
+  for (i = startx; i <= endx; i++, p++, zbuf++,z+=mz) {
+    if (z <= *zbuf) {
+      *p = color;
+      *zbuf = z;
+    }
+  }
+}
+
+
+/* -------------------------------------------------------------------------
+   Plot3D_vertical(Plot3D *p3, int ymin, int ymax, int x, double z1, double z2, Pixel color)
+
+   Draws a "Vertical" line on the framebuffer between two screen coordinates,
+   but also supplies z-values and zbuffering.   This function probably isn't
+   too useful by itself, but will be used by a number of other primitives.
+   -------------------------------------------------------------------------- */
+
+void Plot3D_vertical(Plot3D *p3, int ymin, int ymax, int x, Zvalue z1, Zvalue z2, Pixel color) {
+  Pixel  *p;
+  FrameBuffer *f;
+  int     i;
+  Zvalue  *zbuf,z,mz;
+  int     starty, endy;
+
+  f = p3->frame;
+  if ((x < f->xmin) || (x >= f->xmax)) return;
+  if (ymin >= f->ymax) return;
+  if (ymin < f->ymin) starty = f->ymin;
+  else starty = ymin;
+  if (ymax < f->ymin) return;
+  if (ymax >= f->ymax) endy = f->ymax - 1;
+  else endy = ymax;
+
+  /* Calculate z slope */
+
+  mz = (double) (z2 - z1)/(double) (ymax - ymin);
+
+  /* Draw it */
+  
+  p = &f->pixels[starty][x];
+  zbuf = &f->zbuffer[starty][x];
+  for (i = starty; i <= endy; i++, p+=f->width, zbuf+=f->width) {
+    z = (Zvalue) (mz*(i-ymin) + z1);
+    if (z <= *zbuf) {
+      *p = color;
+      *zbuf = z;
+    }
+  }
+}
+
+/* -------------------------------------------------------------------------------
+   Plot3D_linetransform(Plot3D *p3, int x1, int y1, Zvalue z1,
+                                    int x2, int y2, Zvalue z2, Pixel c)
+
+   Draw a 3D line between points that have already been transformed into
+   3D space.
+
+   Uses a Bresenham line algorithm, but with linear interpolation between
+   Zvalues.
+   ------------------------------------------------------------------------------- */
+   
+void
+Plot3D_linetransform(Plot3D *p3, int x1, int y1, Zvalue z1, int x2, int y2, Zvalue z2, Pixel c) {
+
+  int       orig_x1, orig_y1, orig_x2,orig_y2;
+  Zvalue    zt;
+  
+  /* Bresenham line drawing parameters */
+  FrameBuffer *f;
+  int  dx,dy,dxneg,dyneg, inc1,inc2,di;
+  int  x, y, xpixels, ypixels, xt, yt;
+  Pixel *p;
+  double m;
+  int  end1 = 0, end2 = 0;
+  Zvalue  *zbuf,mz,z;
+
+  f = p3->frame;
+    
+  /* Need to figure out where in the heck this line is */
+
+  dx = x2 - x1;
+  dy = y2 - y1;
+
+  if ((dx == 0) && (dy == 0)) {
+    if ((x1 < f->xmin) || (x1 >= f->xmax) ||
+	(y1 < f->ymin) || (y1 >= f->ymax)) return;
+    if (z1 <= f->zbuffer[y1][x1]) {
+      f->pixels[y1][x1] = c;
+    }
+    return;
+  }
+  if (dx == 0) {
+    /* Draw a Vertical Line */
+    if (y1 < y2)
+      Plot3D_vertical(p3,y1,y2,x1,z1,z2,c);
+    else
+      Plot3D_vertical(p3,y2,y1,x1,z2,z1,c);
+    return;
+  }
+  if (dy == 0) {
+    /* Draw a Horizontal Line */
+    if (x1 < x2)
+      Plot3D_horizontal(p3,x1,x2,y1,z1,z2,c);
+    else
+      Plot3D_horizontal(p3,x2,x1,y1,z2,z1,c);
+    return;
+  }
+
+  /* Figure out where in the heck these lines are using the
+     Cohen-Sutherland Line Clipping Scheme. */
+
+  end1 = ((x1 - f->xmin) < 0) |
+    (((f->xmax- 1 - x1) < 0) << 1) |
+    (((y1 - f->ymin) < 0) << 2) |
+    (((f->ymax-1 - y1) < 0) << 3);
+
+  end2 = ((x2 - f->xmin) < 0) |
+    (((f->xmax-1 - x2) < 0) << 1) |
+    (((y2 - f->ymin) < 0) << 2) |
+    (((f->ymax-1 - y2) < 0) << 3);
+
+  if (end1 & end2) return;      /* Nope : Not visible */
+  
+  /* Make sure points have a favorable orientation */
+  
+  if (x1 > x2) {
+    xt = x1;
+    x1 = x2;
+    x2 = xt;
+    yt = y1;
+    y1 = y2;
+    y2 = yt;
+    zt = z1;
+    z1 = z2;
+    z2 = zt;
+  }
+
+  /* Save original points before we clip them off */
+  orig_x1 = x1;
+  orig_y1 = y1;
+  orig_x2 = x2;
+  orig_y2 = y2;
+  
+  /* Clip against the boundaries */
+  m = (y2 - y1)/(double) (x2-x1);
+  if (x1 < f->xmin) {
+    y1 = (f->xmin - x1)*m + y1;
+    x1 = f->xmin;
+  }
+  if (x2 >= f->xmax) {
+    y2 = (f->xmax -1 -x1)*m + y1;
+    x2 = f->xmax - 1;
+  }
+
+  if (y1 > y2) {
+    xt = x1;
+    x1 = x2;
+    x2 = xt;
+    yt = y1;
+    y1 = y2;
+    y2 = yt;
+    zt = z1;
+    z1 = z2;
+    z2 = zt;
+
+    /* Swap original points */
+
+    xt = orig_x1;
+    orig_x1 = orig_x2;
+    orig_x2 = xt;
+    yt = orig_y1;
+    orig_y1 = orig_y2;
+    orig_y2 = yt;
+  }
+
+  m = 1/m;
+  if (y1 < f->ymin) {
+    x1 = (f->ymin - y1)*m + x1;
+    y1 = f->ymin;
+  }
+  if (y2 >= f->ymax) {
+    x2 = (f->ymax-1-y1)*m + x1;
+    y2 = f->ymax-1;
+  }
+
+  if ((x1 < f->xmin) || (x1 >= f->xmax) || (y1 < f->ymin) || (y1 >= f->ymax) ||
+      (x2 < f->xmin) || (x2 >= f->xmax) || (y2 < f->ymin) || (y2 >= f->ymax)) return;
+
+  dx = x2 - x1;
+  dy = y2 - y1;
+  xpixels = f->width;
+  ypixels = f->height;
+  
+  dxneg = (dx < 0) ? 1 : 0;
+  dyneg = (dy < 0) ? 1 : 0;
+  
+  dx = abs(dx);
+  dy = abs(dy);
+  if (dx >= dy) {
+    /* Slope between -1 and 1. */
+    mz = (z2 - z1)/(orig_x2 - orig_x1);    /* Z interpolation slope */
+    if (dxneg) {
+      x = x1;
+      y = y1;
+      x1 = x2;
+      y1 = y2;
+      x2 = x;
+      y2 = y;
+      dyneg = !dyneg;
+    }
+    inc1 = 2*dy;
+    inc2 = 2*(dy-dx);
+    di = 2*dy-dx;
+
+    /* Draw a line using x as independent variable */
+    
+    p = &f->pixels[y1][x1];
+    zbuf = &f->zbuffer[y1][x1];
+    x = x1;
+    while (x <= x2) {
+      /* Do a z-buffer check */
+      z = mz*(x-orig_x1)+z1;
+      if (z <= *zbuf){
+	*p = c;
+	*zbuf = z;
+      }
+      p++;
+      zbuf++;
+      if (di < 0) {
+	di = di + inc1;
+      } else {
+	if (dyneg) {
+	  p = p - xpixels;
+	  zbuf = zbuf - xpixels;
+	  di = di + inc2;
+	} else {
+	  p = p + xpixels;
+	  zbuf = zbuf + xpixels;
+	  di = di + inc2;
+	}
+      }
+      x++;
+    }
+  } else {
+    /* Slope < -1 or > 1 */
+    mz = (z2 - z1)/(double) (orig_y2 - orig_y1);
+    if (dyneg) {
+      x = x1;
+      y = y1;
+      x1 = x2;
+      y1 = y2;
+      x2 = x;
+      y2 = y;
+      dxneg = !dxneg;
+    }
+    inc1 = 2*dx;
+    inc2 = 2*(dx-dy);
+    di = 2*dx-dy;
+
+    /* Draw a line using y as independent variable */
+    
+    p = &f->pixels[y1][x1];
+    zbuf = &f->zbuffer[y1][x1];
+    y = y1;
+    while (y <= y2) {
+      /* Do a z-buffer check */
+      z = mz*(y-orig_y1)+z1;
+      if (z <= *zbuf) {
+	*p = c;
+	*zbuf = z;
+      }
+      p = p + xpixels;
+      zbuf = zbuf + xpixels;
+      if (di < 0) {
+	di = di + inc1;
+      } else {
+	if (dxneg) {
+	  p = p - 1;
+	  zbuf = zbuf - 1;
+	  di = di + inc2;
+	} else {
+	  p = p + 1;
+	  zbuf = zbuf + 1;
+	  di = di + inc2;
+	}
+      }
+      y++;
+    }
+  }    
+}
+
+/* ---------------------------------------------------------------------------
+   Plot3D_line(Plot3D *p3, double x1, double y1, double z1, double x2, double y2, double z2,int color)
+
+   Draws a line in 3D space.   This is done as follows (for lack of a better
+   method).
+
+   1.  The points (x1,y1,z1) and (x2,y2,z2) are transformed into screen coordinates
+   2.  We draw the line using a modified Bresenham line algorithm.
+   3.  Zbuffer values are linearly interpolated between the two points.
+   ---------------------------------------------------------------------------- */
+
+void
+Plot3D_line(Plot3D *p3, double fx1, double fy1, double fz1, double fx2, double fy2,
+	    double fz2, Pixel c) {
+
+  /* 3D Transformation parameters */
+  GL_Vector t;
+  double    invw;
+  int       x1,y1,x2,y2;
+  Zvalue    z1,z2;
+  
+  /* Transform the two points into screen coordinates */
+  
+  Matrix_transform4(p3->trans_mat,fx1,fy1,fz1,1,&t);  /* Point 1 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  x1 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  y1 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  z1 = t.z;
+  
+  Matrix_transform4(p3->trans_mat,fx2,fy2,fz2,1,&t);  /* Point 2 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  x2 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  y2 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  z2 = t.z;
+  Plot3D_linetransform(p3,x1,y1,z1,x2,y2,z2,c);
+}
+
+
+/* -------------------------------------------------------------------------
+   Plot3D_triangle(Plot3D *p3, double x1, double y1, double z1,
+                               double x2, double y2, double z2,
+			       double x3, double y3, double z3,
+			       Pixel fillcolor)
+
+   This function draws a 3D z-buffered outline triangle.
+   -------------------------------------------------------------------------- */
+
+void Plot3D_triangle(Plot3D *p3, double x1, double y1, double z1,
+		     double x2, double y2, double z2,
+		     double x3, double y3, double z3, Pixel color) {
+
+  int        tx1, tx2, tx3, ty1, ty2, ty3;
+  Zvalue     tz1, tz2, tz3;
+  GL_Vector  t;
+  double     invw;
+  
+  /* Transform the three points into screen coordinates */
+
+  Matrix_transform4(p3->trans_mat,x1,y1,z1,1,&t);  /* Point 1 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx1 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty1 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz1 = (Zvalue) t.z;
+  
+  Matrix_transform4(p3->trans_mat,x2,y2,z2,1,&t);  /* Point 2 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx2 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty2 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz2 = (Zvalue) t.z;
+
+  Matrix_transform4(p3->trans_mat,x3,y3,z3,1,&t);  /* Point 3 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx3 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty3 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz3 = (Zvalue) t.z;
+
+
+  Plot3D_linetransform(p3,tx1,ty1,tz1,tx2,ty2,tz2,color);
+  Plot3D_linetransform(p3,tx1,ty1,tz1,tx3,ty3,tz3,color);
+  Plot3D_linetransform(p3,tx2,ty2,tz2,tx3,ty3,tz3,color);
+}
+
+
+/* -------------------------------------------------------------------------
+   Plot3D_solidtriangletransform(Plot3D *p3, int tx1, int ty2, Zvalue tz1,
+                                 int tx2, int ty2, Zvalue tz2,
+				 int tx3, int ty3, Zvalue tz3, Pixel color) 
+
+   This function draws a 3D z-buffered filled triangle.   Assumes three
+   points have already been transformed into screen coordinates.
+
+   General idea :
+         1.   Transform the three points into screen coordinates
+	 2.   Order three points vertically on screen.
+	 3.   Check for degenerate cases (where 3 points are colinear).
+	 4.   Fill in the resulting triangle using horizontal lines.
+   -------------------------------------------------------------------------- */
+
+void Plot3D_solidtriangletransform(Plot3D *p3, int tx1, int ty1, Zvalue tz1,
+				   int tx2, int ty2, Zvalue tz2,
+				   int tx3, int ty3, Zvalue tz3, Pixel color) {
+  int        tempx, tempy;
+  Zvalue     tempz;
+  double     m1,m2,m3, mz1, mz2, mz3;
+  int        y;
+  int        ix1, ix2;
+  Zvalue     zz1, zz2;
+  FrameBuffer *f;
+  register   double fy1,fy2;
+  register   Zvalue fz1,fz2;
+
+  f = p3->frame;
+
+  /* Check for degenerate cases here */
+
+  if ((ty1 == ty2) && (ty2 == ty3)) {
+    if (tx2 < tx1) {   /* Swap points 1 and 2 if 2 is higher */
+      tempx = tx1;
+      tempz = tz1;
+      tx1 = tx2;
+      tz1 = tz2;
+      tx2 = tempx;
+      tz2 = tempz;
+    }
+    if (tx3 < tx1) {  /* Swap points 1 and 3 if 3 is higher */
+      tempx = tx1;
+      tempz = tz1;
+      tx1 = tx3;
+      tz1 = tz3;
+      tx3 = tempx;
+      tz3 = tempz;
+    }
+    if (tx3 < tx2) { /* Swap points 2 and 3 if 3 is higher */ 
+      tempx = tx2;
+      tempz = tz2;
+      tx2 = tx3;
+      tz2 = tz3;
+      tx3 = tempx;
+      tz3 = tempz;
+    }
+
+    /* Points are aligned horizontally.   Handle as a special case */
+    /* Just draw three lines using the outline color */
+
+    Plot3D_horizontal(p3,tx1,tx2,ty1,tz1,tz3,color);
+
+    /*    Plot3D_linetransform(p3,tx1,ty1,tz1,tx2,ty2,tz2,color);
+    Plot3D_linetransform(p3,tx1,ty1,tz1,tx3,ty3,tz3,color);
+    Plot3D_linetransform(p3,tx2,ty2,tz2,tx3,ty3,tz3,color);
+    */
+
+    return;
+  }
+
+  /* Figure out which point has the greatest "y" value */
+
+  if (ty2 > ty1) {   /* Swap points 1 and 2 if 2 is higher */
+    tempx = tx1;
+    tempy = ty1;
+    tempz = tz1;
+    tx1 = tx2;
+    ty1 = ty2;
+    tz1 = tz2;
+    tx2 = tempx;
+    ty2 = tempy;
+    tz2 = tempz;
+  }
+  if (ty3 > ty1) {  /* Swap points 1 and 3 if 3 is higher */
+    tempx = tx1;
+    tempy = ty1;
+    tempz = tz1;
+    tx1 = tx3;
+    ty1 = ty3;
+    tz1 = tz3;
+    tx3 = tempx;
+    ty3 = tempy;
+    tz3 = tempz;
+  }
+  if (ty3 > ty2) { /* Swap points 2 and 3 if 3 is higher */ 
+    tempx = tx2;
+    tempy = ty2;
+    tempz = tz2;
+    tx2 = tx3;
+    ty2 = ty3;
+    tz2 = tz3;
+    tx3 = tempx;
+    ty3 = tempy;
+    tz3 = tempz;
+  }
+
+  /* Points are now order so that t_1 is the highest point, t_2 is the
+     middle point, and t_3 is the lowest point */
+
+    if (ty2 < ty1) {
+      /* First process line segments between (x1,y1)-(x2,y2)
+	 And between (x1,y1),(x3,y3) */
+
+      m1 = (double) (tx2 - tx1)/(double) (ty2 - ty1);
+      m2 = (double) (tx3 - tx1)/(double) (ty3 - ty1);
+      mz1 = (tz2 - tz1)/(double) (ty2 - ty1);
+      mz2 = (tz3 - tz1)/(double) (ty3 - ty1);
+
+      y = ty1;
+      fy1 = m1*(y-ty1)+0.5 + tx1;
+      fy2 = m2*(y-ty1)+0.5 + tx1;
+      fz1 = mz1*(y-ty1) + tz1;
+      fz2 = mz2*(y-ty1) + tz1;
+      while (y >= ty2) {
+	/* Replace with bresenham scheme */
+	/* Calculate x values from slope */
+	ix1 = (int) fy1;
+	ix2 = (int) fy2;
+	zz1 = fz1;
+	zz2 = fz2;
+	fy1-= m1;
+	fy2-= m2;
+	fz1-= mz1;
+	fz2-= mz2;
+	if (ix1 > ix2) 
+	  Plot3D_horizontal(p3,ix2,ix1,y,zz2,zz1,color);
+	else
+	  Plot3D_horizontal(p3,ix1,ix2,y,zz1,zz2,color);
+	y--;
+      }
+    }
+    if (ty3 < ty2) {
+      /* Draw lower half of the triangle */
+      m2 = (double) (tx3 - tx1)/(double) (ty3 - ty1);
+      m3 = (double) (tx3 - tx2)/(double)(ty3 - ty2);
+      mz2 = (tz3 - tz1)/(double) (ty3 - ty1);
+      mz3 = (tz3 - tz2)/(double) (ty3 - ty2);
+      y = ty2;
+      while (y >= ty3) {
+	ix1 = (int) (m3*(y-ty2)+0.5)+tx2;
+	ix2 = (int) (m2*(y-ty1)+0.5)+tx1;
+	zz1 = mz3*(y-ty2)+tz2;
+	zz2 = mz2*(y-ty1)+tz1;
+	if (ix1 > ix2)
+	  Plot3D_horizontal(p3,ix2,ix1,y,zz2,zz1,color);
+	else
+	  Plot3D_horizontal(p3,ix1,ix2,y,zz1,zz2,color);
+	y--;
+      }
+    }
+}
+
+/* -------------------------------------------------------------------------
+   Plot3D_solidtriangle(Plot3D *p3, double x1, double y1, double z1,
+                               double x2, double y2, double z2,
+			       double x3, double y3, double z3,
+			       Pixel color)
+
+   This function draws a 3D z-buffered filled triangle.    Can be used to
+   draw other primitives such as quadralaterals, etc...
+
+   This function simply transforms the given points and calls
+   Plot3D_SolidTriangleTransform().
+   -------------------------------------------------------------------------- */
+
+void Plot3D_solidtriangle(Plot3D *p3, double x1, double y1, double z1,
+			  double x2, double y2, double z2,
+			  double x3, double y3, double z3, Pixel color) {
+
+  int        tx1, tx2, tx3, ty1, ty2, ty3;
+  Zvalue     tz1, tz2, tz3;
+  GL_Vector  t;
+  double     invw;
+  Matrix     a;
+  register double     xshift, yshift, zoom, width, height, view_xmin, view_ymin;
+
+  a = p3->trans_mat;
+  xshift = p3->xshift;
+  yshift = p3->yshift;
+  zoom = p3->zoom;
+  height = p3->height;
+  width = p3->width;
+  view_xmin = p3->view_xmin;
+  view_ymin = p3->view_ymin;
+
+  /* Transform the three points into screen coordinates */
+
+  t.w = a[12]*x1 + a[13]*y1 + a[14]*z1 + a[15];
+  invw = 1.0/t.w;
+  t.x = (a[0]*x1 + a[1]*y1 + a[2]*z1 + a[3])*invw;
+  t.y = (a[4]*x1 + a[5]*y1 + a[6]*z1 + a[7])*invw;
+  t.z = (a[8]*x1 + a[9]*y1 + a[10]*z1 + a[11])*invw;
+
+  tx1 = (int) ((t.x +xshift)*zoom*width + 0.5) + view_xmin; 
+  ty1 = (int) ((t.y +yshift)*zoom*height + 0.5) + view_ymin;
+  tz1 = (Zvalue) t.z;
+
+
+  t.w = a[12]*x2 + a[13]*y2 + a[14]*z2 + a[15];
+  invw = 1.0/t.w;
+  t.x = (a[0]*x2 + a[1]*y2 + a[2]*z2 + a[3])*invw;
+  t.y = (a[4]*x2 + a[5]*y2 + a[6]*z2 + a[7])*invw;
+  t.z = (a[8]*x2 + a[9]*y2 + a[10]*z2 + a[11])*invw;
+  tx2 = (int) ((t.x +xshift)*zoom*width + 0.5) + view_xmin; 
+  ty2 = (int) ((t.y +yshift)*zoom*height + 0.5) + view_ymin;
+  tz2 = (Zvalue) t.z;
+
+  t.w = a[12]*x3 + a[13]*y3 + a[14]*z3 + a[15];
+  invw = 1.0/t.w;
+  t.x = (a[0]*x3 + a[1]*y3 + a[2]*z3 + a[3])*invw;
+  t.y = (a[4]*x3 + a[5]*y3 + a[6]*z3 + a[7])*invw;
+  t.z = (a[8]*x3 + a[9]*y3 + a[10]*z3 + a[11])*invw;
+  tx3 = (int) ((t.x +xshift)*zoom*width + 0.5) + view_xmin; 
+  ty3 = (int) ((t.y +yshift)*zoom*height + 0.5) + view_ymin;
+  tz3 = (Zvalue) t.z;
+
+  Plot3D_solidtriangletransform(p3,tx1,ty1,tz1,tx2,ty2,tz2,tx3,ty3,tz3,color); 
+  
+}    
+
+
+/* -------------------------------------------------------------------------
+   Plot3D_horizontalinterp(Plot3D *p3, int xmin, int xmax, int y,
+                           double z1, double z2, Pixel c1, Pixel c2)
+
+   Draws a "Horizontal" line on the framebuffer between two screen coordinates,
+   but also supplies z-values and zbuffering.  Performs a color interpolation
+   between c1 and c2.  This is primarily used by the SolidTriangleInterp()
+   function to give the illusion of smooth surfaces.
+   -------------------------------------------------------------------------- */
+
+void Plot3D_horizontalinterp(Plot3D *p3, int xmin, int xmax, int y,
+			     Zvalue z1, Zvalue z2, Pixel c1, Pixel c2) {
+  Pixel  *p;
+  FrameBuffer *f;
+  int     i;
+  Zvalue  *zbuf,z,mz;
+  double   mc;
+  int     startx, endx;
+  double   invdx;
+  
+  f = p3->frame;
+  if ((y < f->ymin) || (y >= f->ymax)) return;
+  if (xmin >= f->xmax) return;
+  if (xmin < f->xmin) startx = f->xmin;
+  else startx = xmin;
+  if (xmax < f->xmin) return;
+  if (xmax >= f->xmax) endx = f->xmax - 1;
+  else endx = xmax;
+
+  /* Calculate z slope */
+  if (xmax != xmin) {
+    invdx = 1.0/(double) (xmax-xmin);
+  } else {
+    invdx = 0;
+  }
+
+  mz = (Zvalue) (z2 - z1)*invdx;
+
+  /* Calculate c slope */
+
+  mc = (double) (c2 - c1)*invdx;
+    
+  /* Draw it */
+  
+  p = &f->pixels[y][startx];
+  zbuf = &f->zbuffer[y][startx];
+  for (i = startx; i <= endx; i++, p++, zbuf++) {
+    z = (Zvalue) (mz*(i-xmin) + z1);
+    if (z <= *zbuf) {
+      *p = (Pixel) (mc*(i-xmin)+c1);
+      *zbuf = z;
+    }
+  }
+}
+
+/* -------------------------------------------------------------------------
+   Plot3D_interptriangletransform(Plot3D *p3,
+                                  int tx1, int ty2, Zvalue tz1, Pixel c1,
+                                  int tx2, int ty2, Zvalue tz2, Pixel c2,
+				  int tx3, int ty3, Zvalue tz3, Pixel c3)
+
+   This function draws a 3D z-buffered filled triangle with color
+   interpolation.   Assumes three points have already been transformed
+   into screen coordinates.
+
+   General idea :
+         1.   Transform the three points into screen coordinates
+	 2.   Order three points vertically on screen.
+	 3.   Check for degenerate cases (where 3 points are colinear).
+	 4.   Fill in the resulting triangle using horizontal lines.
+	 5.   Colors are interpolated between end points
+   -------------------------------------------------------------------------- */
+
+void Plot3D_interptriangletransform(Plot3D *p3,
+				    int tx1, int ty1, Zvalue tz1, Pixel c1,
+				    int tx2, int ty2, Zvalue tz2, Pixel c2,
+				    int tx3, int ty3, Zvalue tz3, Pixel c3) {
+  int        tempx, tempy;
+  Zvalue     tempz;
+  double     m1,m2,m3, mz1, mz2, mz3;
+  double     mc1,mc2,mc3;
+  Pixel      ic1,ic2,tempc;
+  int        y;
+  int        ix1, ix2;
+  Zvalue     zz1, zz2;
+  FrameBuffer *f;
+  
+  f = p3->frame;
+
+  /* Figure out which point has the greatest "y" value */
+
+  if (ty2 > ty1) {   /* Swap points 1 and 2 if 2 is higher */
+    tempx = tx1;
+    tempy = ty1;
+    tempz = tz1;
+    tempc = c1;
+    tx1 = tx2;
+    ty1 = ty2;
+    tz1 = tz2;
+    c1 = c2;
+    tx2 = tempx;
+    ty2 = tempy;
+    tz2 = tempz;
+    c2 = tempc;
+  }
+  if (ty3 > ty1) {  /* Swap points 1 and 3 if 3 is higher */
+    tempx = tx1;
+    tempy = ty1;
+    tempz = tz1;
+    tempc = c1;
+    tx1 = tx3;
+    ty1 = ty3;
+    tz1 = tz3;
+    c1 = c3;
+    tx3 = tempx;
+    ty3 = tempy;
+    tz3 = tempz;
+    c3 = tempc;
+  }
+  if (ty3 > ty2) { /* Swap points 2 and 3 if 3 is higher */ 
+    tempx = tx2;
+    tempy = ty2;
+    tempz = tz2;
+    tempc = c2;
+    tx2 = tx3;
+    ty2 = ty3;
+    tz2 = tz3;
+    c2 = c3;
+    tx3 = tempx;
+    ty3 = tempy;
+    tz3 = tempz;
+    c3 = tempc;
+  }
+
+  /* Points are now order so that t_1 is the highest point, t_2 is the
+     middle point, and t_3 is the lowest point */
+
+  /* Check for degenerate cases here */
+
+  if ((ty1 == ty2) && (ty2 == ty3)) {
+
+    /* Points are aligned horizontally.   Handle as a special case */
+    /* Just draw three lines using the outline color */
+
+    if (tx2 > tx1)
+      Plot3D_horizontalinterp(p3,tx1,tx2,ty1,tz1,tz2,c1,c2);
+    else
+      Plot3D_horizontalinterp(p3,tx2,tx1,ty1,tz2,tz1,c2,c1);
+    if (tx3 > tx1)
+      Plot3D_horizontalinterp(p3,tx1,tx3,ty1,tz1,tz3,c1,c3);
+    else
+      Plot3D_horizontalinterp(p3,tx3,tx1,ty1,tz3,tz1,c3,c1);
+    if (tx3 > tx2)
+      Plot3D_horizontalinterp(p3,tx2,tx3,ty2,tz2,tz3,c2,c3);
+    else
+      Plot3D_horizontalinterp(p3,tx3,tx2,ty2,tz3,tz2,c3,c2);
+    
+  } else {
+
+    /* First process line segments between (x1,y1)-(x2,y2)
+       And between (x1,y1),(x3,y3) */
+
+    if (ty2 < ty1) {
+      m1 = (double) (tx2 - tx1)/(double) (ty2 - ty1);
+      m2 = (double) (tx3 - tx1)/(double) (ty3 - ty1);
+      mz1 = (tz2 - tz1)/(double) (ty2 - ty1);
+      mz2 = (tz3 - tz1)/(double) (ty3 - ty1);
+      mc1 = (c2 - c1)/(double) (ty2 - ty1);
+      mc2 = (c3 - c1)/(double) (ty3 - ty1);
+
+      y = ty1;
+      while (y >= ty2) {
+	/* Calculate x values from slope */
+	ix1 = (int) (m1*(y-ty1)+0.5) + tx1;
+	ix2 = (int) (m2*(y-ty1)+0.5) + tx1;
+	zz1 = mz1*(y-ty1) + tz1;
+	zz2 = mz2*(y-ty1) + tz1;
+	ic1 = mc1*(y-ty1) + c1;
+	ic2 = mc2*(y-ty1) + c1;
+	if (ix1 > ix2) 
+	  Plot3D_horizontalinterp(p3,ix2,ix1,y,zz2,zz1,ic2,ic1);
+	else
+	  Plot3D_horizontalinterp(p3,ix1,ix2,y,zz1,zz2,ic1,ic2);
+	y--;
+      }
+    }
+    if (ty3 < ty2) {
+      /* Draw lower half of the triangle */
+      m2 = (double) (tx3 - tx1)/(double) (ty3 - ty1);
+      mz2 = (tz3 - tz1)/(double) (ty3 - ty1);
+      mc2 = (c3 - c1)/(double) (ty3 - ty1);
+      m3 = (double) (tx3 - tx2)/(double)(ty3 - ty2);
+      mz3 = (tz3 - tz2)/(double) (ty3 - ty2);
+      mc3 = (c3 - c2)/(double) (ty3 - ty2);
+      y = ty2;
+      while (y >= ty3) {
+	ix1 = (int) (m3*(y-ty2)+0.5)+tx2;
+	ix2 = (int) (m2*(y-ty1)+0.5)+tx1;
+	zz1 = mz3*(y-ty2)+tz2;
+	zz2 = mz2*(y-ty1)+tz1;
+	ic1 = mc3*(y-ty2)+c2;
+	ic2 = mc2*(y-ty1)+c1;
+	if (ix1 > ix2)
+	  Plot3D_horizontalinterp(p3,ix2,ix1,y,zz2,zz1,ic2,ic1);
+	else
+	  Plot3D_horizontalinterp(p3,ix1,ix2,y,zz1,zz2,ic1,ic2);
+	y--;
+      }
+    }
+  }
+}
+
+/* -------------------------------------------------------------------------
+   Plot3D_interptriangle(Plot3D *p3,
+                         double x1, double y1, double z1, Pixel c1,
+			 double x2, double y2, double z2, Pixel c2,
+			 double x3, double y3, double z3, Pixel c3)
+
+   This function draws a 3D z-buffered filled triangle with color
+   interpolation.
+
+   This function simply transforms the given points and calls
+   Plot3D_InterpTriangleTransform().
+   -------------------------------------------------------------------------- */
+
+void Plot3D_interptriangle(Plot3D *p3,
+			   double x1, double y1, double z1, Pixel c1,
+			   double x2, double y2, double z2, Pixel c2,
+			   double x3, double y3, double z3, Pixel c3) {
+
+  int        tx1, tx2, tx3, ty1, ty2, ty3;
+  Zvalue     tz1, tz2, tz3;
+  GL_Vector  t;
+  double     invw;
+  
+  /* Transform the three points into screen coordinates */
+
+  Matrix_transform4(p3->trans_mat,x1,y1,z1,1,&t);  /* Point 1 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx1 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty1 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz1 = (Zvalue) t.z;
+  
+  Matrix_transform4(p3->trans_mat,x2,y2,z2,1,&t);  /* Point 2 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx2 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty2 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz2 = (Zvalue) t.z;
+
+  Matrix_transform4(p3->trans_mat,x3,y3,z3,1,&t);  /* Point 3 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx3 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty3 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz3 = (Zvalue) t.z;
+
+  Plot3D_interptriangletransform(p3,tx1,ty1,tz1,c1,tx2,ty2,tz2,c2,tx3,ty3,tz3,c3);
+}    
+
+/* -------------------------------------------------------------------------
+   Plot3D_quad(Plot3D *p3, double x1, double y1, double z1,
+                           double x2, double y2, double z2,
+    		           double x3, double y3, double z3,
+			   double x4, double y4, double z4,
+ 		           Pixel fillcolor)
+
+   This function draws a 3D outlined Quadralateral.  Used primarily for
+   drawing meshes and other things.
+
+   Plotting is done in the following order :
+       (x1,y1,z1)  --> (x2,y2,z2)
+       (x2,y2,z2)  --> (x3,y3,z3)
+       (x3,y3,z3)  --> (x4,y4,z4)
+       (x4,y4,z4)  --> (x1,y1,z1)
+   -------------------------------------------------------------------------- */
+
+void Plot3D_quad(Plot3D *p3, double x1, double y1, double z1,
+		 double x2, double y2, double z2,
+		 double x3, double y3, double z3,
+		 double x4, double y4, double z4,
+		 Pixel color) {
+
+  int        tx1, tx2, tx3, tx4, ty1, ty2, ty3, ty4;
+  Zvalue     tz1, tz2, tz3, tz4;
+  GL_Vector  t;
+  double     invw;
+  
+  /* Transform the three points into screen coordinates */
+
+  Matrix_transform4(p3->trans_mat,x1,y1,z1,1,&t);  /* Point 1 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx1 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty1 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz1 = (Zvalue) t.z;
+  
+  Matrix_transform4(p3->trans_mat,x2,y2,z2,1,&t);  /* Point 2 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx2 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty2 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz2 = (Zvalue) t.z;
+
+  Matrix_transform4(p3->trans_mat,x3,y3,z3,1,&t);  /* Point 3 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx3 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty3 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz3 = (Zvalue) t.z;
+
+  Matrix_transform4(p3->trans_mat,x4,y4,z4,1,&t);  /* Point 3 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx4 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty4 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz4 = (Zvalue) t.z;
+
+  Plot3D_linetransform(p3,tx1,ty1,tz1,tx2,ty2,tz2,color);
+  Plot3D_linetransform(p3,tx2,ty2,tz2,tx3,ty3,tz3,color);
+  Plot3D_linetransform(p3,tx3,ty3,tz3,tx4,ty4,tz4,color);
+  Plot3D_linetransform(p3,tx4,ty4,tz4,tx1,ty1,tz1,color);
+  
+}
+
+
+/* -------------------------------------------------------------------------
+   Plot3D_solidquad(Plot3D *p3, double x1, double y1, double z1,
+                           double x2, double y2, double z2,
+    		           double x3, double y3, double z3,
+			   double x4, double y4, double z4,
+ 		           Pixel fillcolor)
+
+   This function draws a 3D solid Quadralateral. Uses the function
+   Plot3D_SolidTriangleTransform() to fill in the region.
+
+   Plotting is done in the following order :
+       (x1,y1,z1)  --> (x2,y2,z2)
+       (x2,y2,z2)  --> (x3,y3,z3)
+       (x3,y3,z3)  --> (x4,y4,z4)
+       (x4,y4,z4)  --> (x1,y1,z1)
+   -------------------------------------------------------------------------- */
+
+void Plot3D_solidquad(Plot3D *p3, double x1, double y1, double z1,
+		 double x2, double y2, double z2,
+		 double x3, double y3, double z3,
+		 double x4, double y4, double z4,
+		 Pixel color) {
+
+  int        tx1, tx2, tx3, tx4, ty1, ty2, ty3, ty4;
+  Zvalue     tz1, tz2, tz3, tz4;
+  GL_Vector  t;
+  double     invw;
+  
+  /* Transform the three points into screen coordinates */
+
+  Matrix_transform4(p3->trans_mat,x1,y1,z1,1,&t);  /* Point 1 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx1 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty1 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz1 = (Zvalue) t.z;
+  
+  Matrix_transform4(p3->trans_mat,x2,y2,z2,1,&t);  /* Point 2 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx2 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty2 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz2 = (Zvalue) t.z;
+
+  Matrix_transform4(p3->trans_mat,x3,y3,z3,1,&t);  /* Point 3 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx3 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty3 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz3 = (Zvalue) t.z;
+
+  Matrix_transform4(p3->trans_mat,x4,y4,z4,1,&t);  /* Point 3 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx4 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty4 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz4 = (Zvalue) t.z;
+
+  Plot3D_solidtriangletransform(p3,tx1,ty1,tz1,tx2,ty2,tz2,tx3,ty3,tz3,color);
+  Plot3D_solidtriangletransform(p3,tx1,ty1,tz1,tx4,ty4,tz4,tx3,ty3,tz3,color);
+}
+
+/* -------------------------------------------------------------------------
+   Plot3D_interpquad(Plot3D *p3, double x1, double y1, double z1, Pixel c1,
+                           double x2, double y2, double z2, Pixel c2,
+    		           double x3, double y3, double z3, Pixel c3,
+			   double x4, double y4, double z4, Pixel c4)
+
+   This function draws a 3D color-interpolated Quadralateral. Uses the function
+   Plot3D_InterpTriangleTransform() to fill in the region.
+
+   Plotting is done in the following order :
+       (x1,y1,z1)  --> (x2,y2,z2)
+       (x2,y2,z2)  --> (x3,y3,z3)
+       (x3,y3,z3)  --> (x4,y4,z4)
+       (x4,y4,z4)  --> (x1,y1,z1)
+   -------------------------------------------------------------------------- */
+
+void Plot3D_interpquad(Plot3D *p3, double x1, double y1, double z1, Pixel c1,
+		       double x2, double y2, double z2, Pixel c2,
+		       double x3, double y3, double z3, Pixel c3,
+		       double x4, double y4, double z4, Pixel c4) {
+
+
+  int        tx1, tx2, tx3, tx4, ty1, ty2, ty3, ty4;
+  Zvalue     tz1, tz2, tz3, tz4;
+  GL_Vector  t;
+  double     invw;
+  
+  /* Transform the three points into screen coordinates */
+
+  Matrix_transform4(p3->trans_mat,x1,y1,z1,1,&t);  /* Point 1 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx1 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty1 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz1 = (Zvalue) t.z;
+  
+  Matrix_transform4(p3->trans_mat,x2,y2,z2,1,&t);  /* Point 2 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx2 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty2 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz2 = (Zvalue) t.z;
+
+  Matrix_transform4(p3->trans_mat,x3,y3,z3,1,&t);  /* Point 3 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx3 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty3 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz3 = (Zvalue) t.z;
+
+  Matrix_transform4(p3->trans_mat,x4,y4,z4,1,&t);  /* Point 3 */
+  invw = 1.0/t.w;
+  t.x = t.x *invw;
+  t.y = t.y *invw;
+  t.z = t.z *invw;
+  tx4 = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty4 = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz4 = (Zvalue) t.z;
+
+  Plot3D_interptriangletransform(p3,tx1,ty1,tz1,c1,tx2,ty2,tz2,c2,tx3,ty3,tz3,c3);
+  Plot3D_interptriangletransform(p3,tx1,ty1,tz1,c1,tx4,ty4,tz4,c4,tx3,ty3,tz3,c3);
+  
+}
+
+/* --------------------------------------------------------------------------
+   Plot3D_solidsphere(Plot3 *p3, double x, double y, double z, double radius,
+                      Pixel c)
+
+   Makes a 3D sphere at x,y,z with given radius and color. 
+
+   Basic strategy :
+      1.  Transform point to screen coordinates
+      2.  Figure out what the radius is in screen coordinates
+      3.  Use bresenham algorithm for large spheres
+      4.  Use bitmaps for small spheres
+   -------------------------------------------------------------------------- */
+
+/* This is used to fill in spheres */
+static int s_xmin;
+static int s_ymin;
+static int s_xmax;
+static int s_ymax;
+static Pixel **s_pixels;
+static Zvalue **s_zbuffer;
+
+void Plot3D_spherehorizontal(int xmin, int xmax, int y, Zvalue z, Pixel color) {
+  int     i;
+  int     startx, endx;
+  Pixel   *p;
+  Zvalue  *zbuf;
+  
+  if ((y < s_ymin) || (y >= s_ymax)) return;
+  if (xmin < s_xmin) startx = s_xmin;
+  else startx = xmin;
+  if (xmax >= s_xmax) endx = s_xmax - 1;
+  else endx = xmax;
+
+  /* Draw it */
+  
+  p = &s_pixels[y][xmin];
+  zbuf = &s_zbuffer[y][xmin];
+  for (i = startx; i <= endx; i++, p++, zbuf++) {
+    if (z <= *zbuf) {
+      *p = color;
+      *zbuf = z;
+    }
+  }
+}
+
+void Plot3D_solidsphere(Plot3D *p3, double x, double y, double z, double radius,
+			Pixel c) {
+
+  GL_Vector  t,r;
+  double     rad;
+  int        tx,ty, irad;
+  Zvalue     tz;
+  double     invw;
+  int        ix, iy, ix1,ix2,p;
+  FrameBuffer *f;
+
+  /* First transform the point into model coordinates */
+
+  Matrix_transform4(p3->fullmodel_mat,x,y,z,1,&t);  
+
+  /* Now transform two points in order to find proper sphere radius */
+
+  Matrix_transform4(p3->view_mat,t.x+radius,t.y,t.z,t.w,&r);  /* transform radius */
+  Matrix_transform4(p3->view_mat,t.x,t.y,t.z,t.w,&t);
+
+  invw = 1.0/t.w;
+  t.x = t.x*invw;
+  t.y = t.y*invw;
+  t.z = t.z*invw;
+  invw = 1.0/r.w;
+  r.x = r.x*invw;
+  r.y = r.y*invw;
+  r.z = r.z*invw;
+  invw = 1.0/r.w;
+  
+  rad = fabs(t.x - r.x);
+  
+  /* Transform everything into screen coordinates */
+  
+  tx = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz = (Zvalue) t.z;
+  irad = (int) (p3->zoom*(rad*p3->width + 0.5));
+
+  /* This is only a temporary solution (maybe). */
+         
+#define fill_zcircle(x,y,c) \
+  ix1 = tx - x; \
+  ix2 = tx + x; \
+  if (ix1 < s_xmin) ix1 = s_xmin; \
+  if (ix2 >= s_xmax) ix2 = s_xmax; \
+  Plot3D_spherehorizontal(ix1,ix2,y,tz,c);
+  
+  f = p3->frame;
+  s_xmin = f->xmin;
+  s_ymin = f->ymin;
+  s_xmax = f->xmax;
+  s_ymax = f->ymax;
+  s_pixels = f->pixels;
+  s_zbuffer = f->zbuffer;
+  if (irad <= 1) {
+    /* Plot a single pixel */
+    if ((tx >= f->xmin) && (tx < f->xmax)) {
+      if ((ty >= f->ymin) && (ty <f->ymax)) {
+        if (tz <= f->zbuffer[ty][tx]) {
+	  f->pixels[ty][tx] = c;
+	  f->zbuffer[ty][tx] = tz;
+	}
+      }
+    }
+    return;
+  }
+  ix = 0;
+  iy = irad;
+  p = 3-2*irad;
+  while (ix <= iy) {
+    fill_zcircle(ix,ty+iy,c);
+    fill_zcircle(ix,ty-iy,c);
+    fill_zcircle(iy,ty+ix,c);
+    fill_zcircle(iy,ty-ix,c);
+    if (p < 0) p = p + 4*ix + 6;
+    else {
+      p = p + 4*(ix-iy) + 10;
+      iy = iy -1;
+    }
+    ix++;
+  }
+}
+
+
+/* --------------------------------------------------------------------
+   Plot3D_outlinesphere(Plot3D *p3, double x, double y, double z,
+                        double radius, Pixel color, Pixel bc)
+
+   Draws an outlined sphere.
+   -------------------------------------------------------------------- */
+
+void Plot3D_outlinesphere(Plot3D *p3, double x, double y, double z,
+			  double radius, Pixel c, Pixel bc)
+{
+  GL_Vector  t,r;
+  double     rad;
+  int        tx,ty, irad;
+  Zvalue     tz;
+  double     invw;
+  int        ix, iy, ix1,ix2,p;
+  
+  FrameBuffer *f;
+
+  /* First transform the point into model coordinates */
+
+  Matrix_transform4(p3->fullmodel_mat,x,y,z,1,&t);  
+
+  /* Now transform two points in order to find proper sphere radius */
+
+  Matrix_transform4(p3->view_mat,t.x+radius,t.y,t.z,t.w,&r);  /* transform radius */
+  Matrix_transform4(p3->view_mat,t.x,t.y,t.z,t.w,&t);
+
+  invw = 1.0/t.w;
+  t.x = t.x*invw;
+  t.y = t.y*invw;
+  t.z = t.z*invw;
+  invw = 1.0/r.w;
+  r.x = r.x*invw;
+  r.y = r.y*invw;
+  r.z = r.z*invw;
+  invw = 1.0/r.w;
+  
+  rad = fabs(t.x - r.x);
+  
+  /* Transform everything into screen coordinates */
+  
+  tx = (int) ((t.x +p3->xshift)*p3->zoom*p3->width + 0.5) + p3->view_xmin; 
+  ty = (int) ((t.y +p3->yshift)*p3->zoom*p3->height + 0.5) + p3->view_ymin;
+  tz = (Zvalue) t.z;
+  irad = (int) (p3->zoom*(rad*p3->width + 0.5));
+
+  /* This is only a temporary solution (maybe). */
+#define plot_zcircle(x,y,c) \
+   if ((x >= s_xmin) && (x < s_xmax) && \
+       (y >= s_ymin) && (y < s_ymax)) {\
+         if (tz <= s_zbuffer[y][x]) { \
+             s_pixels[y][x] = c; \
+	     s_zbuffer[y][x] = tz; } \
+   }
+     
+  f = p3->frame;
+  s_xmin = f->xmin;
+  s_ymin = f->ymin;
+  s_xmax = f->xmax;
+  s_ymax = f->ymax;
+  s_pixels = f->pixels;
+  s_zbuffer = f->zbuffer;
+  
+  if (irad <= 1) {
+    /* Plot a single pixel */
+    if ((tx >= f->xmin) && (tx < f->xmax)) {
+      if ((ty >= f->ymin) && (ty <f->ymax)) {
+        if (tz <= f->zbuffer[ty][tx]) {
+	  f->pixels[ty][tx] = c;
+	  f->zbuffer[ty][tx] = tz;
+	}
+      }
+    }
+    return;
+  }
+  ix = 0;
+  iy = irad;
+  p = 3-2*irad;
+  while (ix <= iy) {
+    fill_zcircle(ix,ty+iy,c);
+    fill_zcircle(ix,ty-iy,c);
+    fill_zcircle(iy,ty+ix,c);
+    fill_zcircle(iy,ty-ix,c);
+    
+    plot_zcircle(tx+ix,ty+iy,bc);
+    plot_zcircle(tx-ix,ty+iy,bc);
+    plot_zcircle(tx+ix,ty-iy,bc);
+    plot_zcircle(tx-ix,ty-iy,bc);
+    plot_zcircle(tx+iy,ty+ix,bc);
+    plot_zcircle(tx-iy,ty+ix,bc);
+    plot_zcircle(tx+iy,ty-ix,bc);
+    plot_zcircle(tx-iy,ty-ix,bc);
+    if (p < 0) p = p + 4*ix + 6;
+    else {
+      p = p + 4*(ix-iy) + 10;
+      iy = iy -1;
+    }
+    ix++;
+  }
+}
+
+/* QUAD Test
+   Test out quad functions for graphing */
+
+double zf(double x, double y) {
+  return cos(sqrt(x*x + y*y)*10.0)/(sqrt(x*x+y*y)+1);
+}
+
+void Quad_Test(Plot3D *p3, int npoints) {
+  int i,j;
+  double dx;
+  double x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,za;
+  int c;
+  dx = 2.0/npoints;
+
+  
+  for (i = 0; i < npoints; i++)
+    for (j = 0; j < npoints; j++) {
+      x1 = i*dx + -1.0;
+      y1 = j*dx + -1.0;
+      x2 = x1 + dx;
+      x3 = x1 + dx;
+      x4 = x1;
+      y2 = y1;
+      y3 = y1 + dx;
+      y4 = y1 + dx;
+      z1 = zf(x1,y1);
+      z2 = zf(x2,y2);
+      z3 = zf(x3,y3);
+      z4 = zf(x4,y4);
+      za = 0.25*(z1+z2+z3+z4);
+      c = 16+((za + 1)*120);
+      if (c > 254) c = 254;
+      Plot3D_quad(p3,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,(Pixel) c);
+    }
+}
+
+
+void Quad_SolidTest(Plot3D *p3, int npoints) {
+  int i,j;
+  double dx;
+  double x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,za;
+  int c;
+  dx = 2.0/npoints;
+
+  
+  for (i = 0; i < npoints; i++)
+    for (j = 0; j < npoints; j++) {
+      x1 = i*dx + -1.0;
+      y1 = j*dx + -1.0;
+      x2 = x1 + dx;
+      x3 = x1 + dx;
+      x4 = x1;
+      y2 = y1;
+      y3 = y1 + dx;
+      y4 = y1 + dx;
+      z1 = zf(x1,y1);
+      z2 = zf(x2,y2);
+      z3 = zf(x3,y3);
+      z4 = zf(x4,y4);
+      za = 0.25*(z1+z2+z3+z4);
+      c = 16+((za + 1)*120);
+      if (c > 254) c = 254;
+      Plot3D_solidquad(p3,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,(Pixel) c);
+    }
+}
+
+
+
+void Quad_InterpTest(Plot3D *p3, int npoints) {
+  int i,j;
+  double dx;
+  double x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4;
+  int c1,c2,c3,c4;
+  dx = 2.0/npoints;
+
+  
+  for (i = 0; i < npoints; i++)
+    for (j = 0; j < npoints; j++) {
+      x1 = i*dx + -1.0;
+      y1 = j*dx + -1.0;
+      x2 = x1 + dx;
+      x3 = x1 + dx;
+      x4 = x1;
+      y2 = y1;
+      y3 = y1 + dx;
+      y4 = y1 + dx;
+      z1 = zf(x1,y1);
+      z2 = zf(x2,y2);
+      z3 = zf(x3,y3);
+      z4 = zf(x4,y4);
+      c1 = 16+((z1 + 1)*120);
+      c2 = 16+((z2 + 1)*120);
+      c3 = 16+((z3 + 1)*120);
+      c4 = 16+((z4 + 1)*120);
+      if (c1 > 254) c1 = 254;
+      if (c2 > 254) c2 = 254;
+      if (c3 > 254) c3 = 254;
+      if (c4 > 254) c4 = 254;
+      Plot3D_interpquad(p3,x1,y1,z1,(Pixel) c1,x2,y2,z2,(Pixel) c2,x3,y3,z3,(Pixel) c3,x4,y4,z4,(Pixel) c4);
+    }
+}
+
+
+
+  
+  
+  
+   
+
+
+  
+  
+