float fd = 5.0;                        // focal length of the camera lens
float sc = 200.0;                      // scale to draw at
int windowHeight = 500;                // size of the window
int windowWidth = 750;
float windowCenterY = windowHeight/2;
float windowCenterX = windowWidth/2;

/* a linear transformation is represented as a 3x3 matrix 
   We use these to express rotations in three-dimensional space 
   
   A 3x3 matrix can be seen as a transformation from one 
   coordinate system into another--and in a very natural way.
   If I have three vectors (u,v,w) in space, they define a 
   coordinate system--they are the three axes against which to 
   measure other points. Any point in space can be seen as 
   so-and-so many times the length of u (and in the direction 
   of u), then so-and-so many times v, and so on.
   
   To construct the matrix that transforms a point from that 
   co-ordinate system into the "standard" one, we just fill in
   the first column of the matrix with the components of u, the
   second column with those of v, and the last with those of w.
   Then the standard operation of matrix multiplication will 
   do the transformation.
   */
class LinTrans {
  float ax, ay, az,
        bx, by, bz,
        cx, cy, cz;
  LinTrans(float _ax, float _ay, float _az, float _bx, float _by, float _bz, float _cx, float _cy, float _cz) {
    ax = _ax; ay = _ay; az = _az;
    bx = _bx; by = _by; bz = _bz;
    cx = _cx; cy = _cy; cz = _cz;
  }
}

class ZeroVectorException extends RuntimeException { };
class BadArgException extends RuntimeException { };

//*  Rats: enum doesn't seem to work in Processing.
// enum Axis { X, Y, Z };

// Point3 represents a point or vector in 3-space
class Point3 {
  
  final static int UNITX = 1;
  final static int UNITY = 2;
  final static int UNITZ = 3;
  public float x, y, z;

  // A const'r for the zero vector:   Point3 zeroVec = Point3(0);
  Point3(int zero) {
    if (zero == 0)
      x = y = z = 0;
    else throw new BadArgException();
  }
  // Construct a point from three coords.
  Point3(float ix, float iy, float iz) { x=ix;  y=iy;  z=iz;  };
  Point3 clone() {
    Point3 rslt = new Point3(x, y, z);
    return rslt;
  };

  Point3(Point3 copy) {
    x = copy.x; y = copy.y;  z = copy.z;
  };

  boolean isZero() {
    return (x == 0. && y == 0. && z == 0.);
  }

  float magnitude() {
    return sqrt(x * x + y * y + z * z);
  };

  void scale(float s) {
    x *= s; y *= s; z *= s;
  }

  void scaledown(float s) {
    x /= s; y /= s; z /= s;
  }

  // [invert]
  // Inverts the vector (acts in-place).
  void invert() {
    x = -x; y = -y; z = -z;
  }
  
  // Alter the target vector to make it unit length in the same direction
  void normalize() {
    float m = magnitude();
    if (m == 0.) throw new ZeroVectorException();
    scaledown(m);
  };

  // Return a clone of the target vector, having unit length in the same direction
  Point3 normalClone() {
    float m = magnitude();
    if (m == 0.) throw new ZeroVectorException();
    Point3 daClone = clone();
    daClone.scaledown(m);
    return daClone;
  };
  
  // If the target vector is longer than the given limit, scale it down to just that limit (same direction)
  void limitMag(float limit) {
    float mag = magnitude();
    if (mag > limit) {
      normalize();
      scale(limit);
    }
  }
 
  // rotateX, rotateY
  //   rotate the point by the given angle around the given axis through the given center
  void rotateY(Point3 center, float angle)
  {
    float nx, ny, nz;
    nx = cos(angle) * (x-center.x) + sin(angle) * (z - center.z) + center.x;
    ny = y;
    nz = -sin(angle) * (x-center.x) + cos(angle) * (z - center.z) + center.z;
    x = nx; y = ny; z = nz;
  };
  
  void rotateX(Point3 center, float angle)
  {
    float nx, ny, nz;
    nx = x;
    ny = cos(angle) * (y-center.y) + sin(angle) * (y - center.y) + center.y;
    nz = -sin(angle) * (z-center.z) + cos(angle) * (z - center.z) + center.z;
    x = nx; y = ny; z = nz;
  }
  
  // a special way of constructing a vector perpendicular to a certain axis and \
  //  making a certain angle with the other two.
  Point3(int selector, float angle) {
    switch (selector) {
      case UNITX: { x = 0.; y = sin(angle); z = cos(angle); break; }
      case UNITY: { x = sin(angle); y = 0.; z = cos(angle); break; }
      case UNITZ: { x = sin(angle); y = cos(angle); z = 0.; break; }
    }
  }

  void translate(float dx, float dy, float dz) {
     x += dx; y += dy; z += dz; };

  void translate(Point3 offset) {
     x += offset.x; y += offset.y; z += offset.z; };
  
  Point3 translateClone(Point3 offset) {
    Point3 result = clone();
    result.x += offset.x; result.y += offset.y; result.z += offset.z; 
    return result;
  };

  void subtract(Point3 p) {
    x -= p.x; y -= p.y; z -= p.z;
  }

  void linTransform(LinTrans tr) {
    float nx, ny, nz;
    nx = tr.ax * x  +  tr.ay * y  +  tr.az * z;
    ny = tr.bx * x  +  tr.by * y  +  tr.bz * z;
    nz = tr.cx * x  +  tr.cy * y  +  tr.cz * z;
    x = nx; y = ny; z = nz;
  };
}

// Return the vector which is the cross product of the two vectors.
// The cross product is always perpendicular to both arguments
// Its magnitude is equal to the product of the two arguments' 
// magnitudes times the sine of the angle between them.
Point3 cross(Point3 u, Point3 v) {
    return new Point3(u.y * v.z - u.z * v.y,
                      u.z * v.x - u.x * v.z,
                      u.x * v.y - u.y * v.x);
};

final Point3 origin = new Point3(0,0,0);

// Return the dot product of two vectors
// The dot product is defined as the sum of the products of the corresponding components
// It is always equal to:
//   * the product of the magnitudes of the two vectors, times the cosine of the angle between them
//   * the projection of either vector onto the line of the other, given as a multiple of that other's length
float dot(Point3 u, Point3 v) {
  return (u.x * v.x + u.y * v.y + u.z * v.z);
};

// [center]
// Given a set of points, return their arithmetic mean, or geographic center.
Point3 center(Point3[] pts) {
  Point3 result = new Point3(0);
  for (int i = 0; i < pts.length; i++) {
    result.translate(pts[i]);
  }
  result.scaledown(pts.length);
  return result;
}

// [uprightRotation]
// Given an orientation, uprightRotation constructs the rotation 
// that takes the z axis to that orientation. We have a choice of
// how to "roll" the object around that axis, so we try to keep it
// upright.
LinTrans uprightRotation(Point3 v) {
  Point3 zVec = v.normalClone();
  // The trick is as follows:
  // First we cross the orientation with the positive y axis;
  // this gives us a vector in the xz plane that's perpendicular
  // to the desired orientation--ultimately the "rightward" direction
  // of the rotation.
  Point3 up = new Point3(0, 1, 0);
  Point3 xVec = cross(zVec, up);
  xVec.normalize();
  // Next we cross the result with the original orientation to get
  // a new "upward" vector--but this one is perpendicular to the 
  // original v
  Point3 yVec = cross(xVec, zVec);

  // Now we just fill in the matrix of the linear transformation.
  return new LinTrans(xVec.x, yVec.x, zVec.x,
                      xVec.y, yVec.y, zVec.y,
                      xVec.z, yVec.z, zVec.z);
}

// projectX/projectY
//   Given a point p in space and a focal length fd, projects the point (through 
//   the origin) onto a plane perpendicular to the z axis (that plane where z = fd).
//   A larger fd is telephoto; smaller is wideangle.
//   The [sc] parameter chooses the scale of the picture plane: It is the number of
//   screen-distance units that correspond to one space-unit if it is placed at a 
//   distance of one focal length from the observer.
float projectX(Point3 p, float fd, float sc)
{
  return (p.x/p.z)*fd*sc;
}

float projectY(Point3 p, float fd, float sc)
{
  return (p.y/p.z)*fd*sc;
}

// [sline]
// Given two points in 3-space, draw the line between them, as projected into
// the picture plane. The focal length used is a global variable :-(
// Since the "view" transformation is a linear transformation, we know that
// a straight line in space projects to a straight line on the picture plane.
void sline(Point3 p1, Point3 p2)
{
  if (p1.z > 0 && p2.z > 0) 
    line(projectX(p1, fd, sc)+windowCenterX, windowHeight-(projectY(p1, fd, sc)+windowCenterY),
         projectX(p2, fd, sc)+windowCenterX, windowHeight-(projectY(p2, fd, sc)+windowCenterY));
}

Point3 rotateY(Point3 p, Point3 center, float currentTheta)
{
  Point3 rslt;
  rslt = new Point3(
    cos(currentTheta) * (p.x-center.x) + sin(currentTheta) * (p.z - center.z) + center.x,
    p.y,
    -sin(currentTheta) * (p.x-center.x) + cos(currentTheta) * (p.z - center.z) + center.z);
  return rslt;
}

// [direction]
// Given two points, returns a normal vector indicating the direction from 
//   the first to the second.
Point3 direction(Point3 a, Point3 b) {
  Point3 result = new Point3(b.x - a.x, b.y - a.y, b.z - a.z);
  result.normalize();
  return result;
}

// [invert]
// Constructs the additive inverse of the given vector.
Point3 invert(Point3 pt) {
  return new Point3(-pt.x, -pt.y, -pt.z);
}

/*
Point3 translate(Point3 p, float x, float y, float z) {
  Point3 rslt; 
  rslt = new Point3(p.x + x, p.y + y, p.z + z);
  return rslt;
}

Point3 translate(Point3 p, Point3 delta) { 
  return translate(p, p.x, p.y, p.z); 
}
*/

// Some useful constants.
float pi = 3.14;
float halfCircle = pi;
float fullCircle = 2 * pi;

/*
float pin(float x, float min, float max) {
  return (x > max ? max : (x < min ? min : x));
}

float sign(float x) {
  if (x > 0.0) return 1.0;
  if (x < 0.0) return -1.0;
  if (x == 0.0) return 0.0;
  return -33.73;
}
*/

class Face {
  public Point3[] vertices;
  Point3 normal() {
    if (vertices.length > 2) {
      Point3 a = direction(vertices[0], vertices[1]);
      Point3 b = direction(vertices[0], vertices[2]);
      return cross(a, b);
    } else throw new ZeroVectorException();
  }
  Point3 center() {
    Point3 result = new Point3(0);
    for (int i=0; i < vertices.length; i++) 
      result.translate(vertices[i]);
    result.scaledown(vertices.length);
    return result;
  }
  Face clone() {
    Face rslt = new Face();
    rslt.vertices = (Point3[])vertices.clone();
    return rslt;
  }
}

class Edge {
  public Point3 a, b;
}

class IntPair {
  public int a, b; 
  IntPair(int _a, int _b) { a = _a; b = _b; };
  public IntPair clone() { return new IntPair(a, b); } 
}

class Volume {
  public Point3[] vertices;
  public IntPair[] edges;
  public int[][] faces;

  Edge[] edges() {
    Edge[] result;
    result = new Edge[edges.length];
    for (int i = 0; i < edges.length; i++) {
      result[i] = new Edge();
      result[i].a = vertices[edges[i].a];
      result[i].b = vertices[edges[i].b];
    }
    return result;
  }
  
  Face[] faces() {
    Face[] rslt = new Face[faces.length];
    for (int i = 0; i < faces.length; i++) {
      rslt[i] = new Face();
      rslt[i].vertices = new Point3[faces[i].length];
      for (int j = 0; j < faces[i].length; j++) {
        rslt[i].vertices[j] = new Point3(vertices[faces[i][j]]);
      }
    }
    return rslt;
  }

  void rotateY(Point3 center, float angle) {
    for (int i = 0; i < vertices.length; i++) {
      vertices[i].rotateY(center, angle);
    }
  }
  
  void translate(float dx, float dy, float dz) {
    for (int i = 0; i < vertices.length; i++) {
      vertices[i].translate(dx, dy, dz);
    }
  }
  
  void translate(Point3 offset) {
    for (int i = 0; i < vertices.length; i++) {
      vertices[i].translate(offset);
    }
  }
  
  // [linTrans]
  // Apply a linear transformation to the volume as a whole.
  void linTrans(LinTrans tx) {
    for (int i = 0; i < vertices.length; i++) {
      vertices[i].linTransform(tx);
    }
  }
  
  Volume clone() {
    Volume rslt = new Volume();
    rslt.vertices = new Point3[vertices.length];
    for (int i = 0; i < vertices.length; i++) {
       rslt.vertices[i] = vertices[i].clone();
    }
    rslt.edges = new IntPair[edges.length];
    for (int i = 0; i < edges.length; i++) {
      rslt.edges[i] = edges[i].clone();
    }
//    rslt.faces = (int[][])faces.clone();
/*
    new Face[faces.length];
    for (int i = 0; i < faces.length; i++) {
      rslt.faces[i] = faces[i].clone();
    }
*/
    return rslt;
  }
}

class Light {
  public Point3 position;
  Light(Point3 _pos) { position = _pos; }
}

/*
// There is a bug with Processing and jikes concerning interfaces, so I'm not using this for now.
interface RenderMethod {
  void render(Volume vol);
}
*/
class WfRenderMethod /*implements RenderMethod*/ {
  void render(Volume vol) {
    Edge[] edges;
    edges = vol.edges();
    for (int i = 0; i < edges.length; i++) {
      sline(edges[i].a, edges[i].b);
    }
  }
}

class FlatRenderMethod /*implements RenderMethod*/ { 
  Light theLight;
  FlatRenderMethod(Light lt) { theLight = lt; }
  void render(Volume vol) {
    Face[] faces = vol.faces();
    for (int i = 0; i < faces.length; i++) {
      renderFace(faces[i]);
    }
  }
  void renderFace(Face f) {
    Point3 faceNormal = f.normal();
    Point3 toLight = direction(f.vertices[0], theLight.position);
    float lightIndex = dot(faceNormal, toLight);
    noStroke();
    fill((lightIndex + 0.15) * 255);
    beginShape();
    for (int i = 0; i < f.vertices.length; i++) {
      vertex(projectX(f.vertices[i], fd, sc)+windowCenterX,
             projectY(f.vertices[i], fd, sc)+windowCenterY);
    }
    vertex(projectX(f.vertices[0], fd, sc)+windowCenterX,
           projectY(f.vertices[0], fd, sc)+windowCenterY);
    endShape();
  }
}

class Camera {
    public Point3 offset, orient;
   // offset is the inverse of the camera's position--the displacement by which a point should be translated to get its position relative to the camera.
}

void renderWithCamera(WfRenderMethod rm, Volume vol, Camera cam) {
//void renderWithCamera(FlatRenderMethod rm, Volume vol, Camera cam) {
  Volume workVol = vol.clone();
  workVol.translate(cam.offset);
  workVol.linTrans(uprightRotation(cam.orient));
  rm.render(workVol);
}

class IsNoneExc extends RuntimeException { } 

class MaybeFloat {
  public boolean isSome;
  float value;
  MaybeFloat() { isSome = false; };
  MaybeFloat(float _value) { isSome = true; value = _value; };
  MaybeFloat(boolean _isSome, float _value) { isSome = _isSome; value = _value; };
  float valOf() { if (isSome) { return value; } else {throw new IsNoneExc(); } }
}

MaybeFloat isInCone(Point3 center, Point3 dir, float angle, Point3 target) {
  Point3 dispmnt = direction(center, target);
  Point3 offset = cross(dispmnt, dir);
  if (offset.magnitude() < abs(sin(angle))) {
    return new MaybeFloat(dot(dir, dispmnt));
  } else return new MaybeFloat();
}

boolean isInTruncCone(Point3 center, Point3 dir, float angle, float depth, Point3 target) {
  MaybeFloat d = isInCone(center, dir, angle, target);
  return (d.isSome && d.valOf() <= depth);
}

//interface PointPred { boolean test(Point3 p); }

boolean pointInSphere(Point3 center, float radius, Point3 p) {
  Point3 c = center.clone();
  c.subtract(p);
  if (c.isZero()) return true;
  return (c.magnitude() <= radius);
}

/* The initBLAH routines define models--that is, Volume */
Volume initCube() {
  Volume cube = new Volume();
  cube.vertices = new Point3[8];
  cube.vertices[0] = new Point3(-10.0, 10.0, -5.0);
  cube.vertices[1] = new Point3( 10.0, 10.0, -5.0);
  cube.vertices[2] = new Point3(-10.0,-10.0, -5.0);
  cube.vertices[3] = new Point3( 10.0,-10.0, -5.0);
  cube.vertices[4] = new Point3(-10.0, 10.0,  5.0);
  cube.vertices[5] = new Point3( 10.0, 10.0,  5.0);
  cube.vertices[6] = new Point3(-10.0,-10.0,  5.0);
  cube.vertices[7] = new Point3( 10.0,-10.0,  5.0);
  cube.edges = new IntPair[12]; 
  cube.edges[ 0] = new IntPair(0, 1);
  cube.edges[ 1] = new IntPair(0, 2);
  cube.edges[ 2] = new IntPair(0, 4);
  cube.edges[ 3] = new IntPair(1, 3);
  cube.edges[ 4] = new IntPair(1, 5);
  cube.edges[ 5] = new IntPair(2, 3);
  cube.edges[ 6] = new IntPair(2, 6);
  cube.edges[ 7] = new IntPair(3, 7);
  cube.edges[ 8] = new IntPair(4, 5);
  cube.edges[ 9] = new IntPair(4, 6);
  cube.edges[10] = new IntPair(5, 7);
  cube.edges[11] = new IntPair(6, 7);

  return cube;
}

Volume initWedge() {
  Volume o = new Volume();
  o.vertices = new Point3[4];
  o.vertices[0] = new Point3(  0.0, 0.0, 5.0);
  o.vertices[1] = new Point3(  1.5, 0.0, 0.0);
  o.vertices[2] = new Point3( -1.5, 0.0, 0.0);
  o.vertices[3] = new Point3(  0.0, 2.0, 0.0);
  o.edges = new IntPair[6]; 
  o.edges[ 0] = new IntPair(0, 1);
  o.edges[ 1] = new IntPair(0, 2);
  o.edges[ 2] = new IntPair(0, 3);
  o.edges[ 3] = new IntPair(1, 2);
  o.edges[ 4] = new IntPair(1, 3);
  o.edges[ 5] = new IntPair(2, 3);

  o.faces = new int[4][3];
  o.faces[0][0] = 0; o.faces[0][1] = 1; o.faces[0][2] = 2;
  o.faces[1][0] = 0; o.faces[1][1] = 1; o.faces[1][2] = 3;
  o.faces[2][0] = 0; o.faces[2][1] = 2; o.faces[2][2] = 3;
  o.faces[3][0] = 1; o.faces[3][1] = 2; o.faces[3][2] = 3;

/*
  o.faces = new Face[3];
  o.faces[0] = new Face();
  o.faces[0].vertices = new Point3[3];
  o.faces[0].vertices[0] = new Point3( 0.0, 0.0, 5.0);
  o.faces[0].vertices[1] = new Point3( 1.5, 0.0, 0.0);
  o.faces[0].vertices[2] = new Point3(-1.5, 0.0, 0.0);
  o.faces[1] = new Face();
  o.faces[1].vertices = new Point3[3];
  o.faces[1].vertices[0] = new Point3( 0.0, 0.0, 5.0);
  o.faces[1].vertices[1] = new Point3( 1.5, 0.0, 0.0);
  o.faces[1].vertices[2] = new Point3( 0.0, 2.0, 0.0);
  o.faces[2] = new Face();
  o.faces[2].vertices = new Point3[3];
  o.faces[2].vertices[0] = new Point3( 0.0, 0.0, 5.0);
  o.faces[2].vertices[1] = new Point3(-1.5, 0.0, 0.0);
  o.faces[2].vertices[2] = new Point3( 0.0, 2.0, 0.0);
*/
  return o;
}

/*
Volume initWedgeSolid() {
  o.vertices = new Point3[4];
  o.vertices[0] = new Point3(  0.0, 0.0, 5.0);
  o.vertices[1] = new Point3(  1.5, 0.0, 0.0);
  o.vertices[2] = new Point3( -1.5, 0.0, 0.0);
  o.vertices[3] = new Point3(  0.0, 2.6, 0.0);
  o.edges = new IntPair[6]; 
  o.edges[ 0] = new IntPair(0, 1);
  o.edges[ 1] = new IntPair(0, 2);
  o.edges[ 2] = new IntPair(0, 3);
  o.edges[ 3] = new IntPair(1, 2);
  o.edges[ 4] = new IntPair(1, 3);
  o.edges[ 5] = new IntPair(2, 3);

  return o;
}
*/

Volume initGroundX() {
  Volume o = new Volume();
  o.vertices = new Point3[4];
  o.vertices[0] = new Point3( 20.0, 0.0,  20.0);
  o.vertices[1] = new Point3( 20.0, 0.0, -20.0);
  o.vertices[2] = new Point3(-20.0, 0.0, -20.0);
  o.vertices[3] = new Point3(-20.0, 0.0,  20.0);
  o.edges = new IntPair[2];
  o.edges[ 0] = new IntPair(0, 2);
  o.edges[ 1] = new IntPair(3, 1);
  return o;
}

// The "woods circle" is a circle on the ground that shows the bounds of the 
// "woods" in which the drubbins live. They don't likie to go outside of it.
Volume initWoodsCircle() {
  Volume o = new Volume();
  o.vertices = new Point3[12];
  o.vertices[ 0] = new Point3(  0.0, 0.0, 100.0);
  o.vertices[ 1] = new Point3( 50.0, 0.0,  87.0);
  o.vertices[ 2] = new Point3( 87.0, 0.0,  50.0);
  o.vertices[ 3] = new Point3(100.0, 0.0,   0.0);
  o.vertices[ 4] = new Point3( 87.0, 0.0, -50.0);
  o.vertices[ 5] = new Point3( 50.0, 0.0, -87.0);
  o.vertices[ 6] = new Point3(  0.0, 0.0,-100.0);
  o.vertices[ 7] = new Point3(-50.0, 0.0, -87.0);
  o.vertices[ 8] = new Point3(-87.0, 0.0, -50.0);
  o.vertices[ 9] = new Point3(-100.0, 0.0,  0.0);
  o.vertices[10] = new Point3(-87.0, 0.0,  50.0);
  o.vertices[11] = new Point3(-50.0, 0.0,  87.0);
  o.edges = new IntPair[12]; 
  o.edges[ 0] = new IntPair(0, 1);
  o.edges[ 1] = new IntPair(1, 2);
  o.edges[ 2] = new IntPair(2, 3);
  o.edges[ 3] = new IntPair(3, 4);
  o.edges[ 4] = new IntPair(4, 5);
  o.edges[ 5] = new IntPair(5, 6);
  o.edges[ 6] = new IntPair(6, 7);
  o.edges[ 7] = new IntPair(7, 8);
  o.edges[ 8] = new IntPair(8, 9);
  o.edges[ 9] = new IntPair(9, 10);
  o.edges[10] = new IntPair(10, 11);
  o.edges[11] = new IntPair(11, 0);
  return o;
}

// [drubbinate]:
//   Defines the behavior of the drubbins.
Point3[] drubbinate(Point3[] drubbinPos, Point3[] drubbinVels, int drubbinCount) {
  Point3[] nextVel = (Point3[])drubbinVels.clone();
  
  Point3 drubbinCenter = center(drubbinPos);
//  print("drubbinCenter: (" + drubbinCenter.x + ", " + drubbinCenter.y + ", " + drubbinCenter.z + ")\n");
  
  for (int i = 0; i < drubbinCount; i++) {
    if (drubbinPos[i].magnitude() >= woodsRadius) {
      // If we're out of the woods, head toward home
      Point3 adjustment = direction(drubbinPos[i], origin);
      adjustment.scale(.1);
      nextVel[i].translate(adjustment);
    } else {
      // otherwise, head toward the center of the flock.
      Point3 idealHeading = direction(drubbinPos[i], drubbinCenter);
      Point3 adjIndex = cross(drubbinVels[i], idealHeading);
      Point3 adjustment = cross(adjIndex, drubbinVels[i]);
      adjustment.scale(.1);
      nextVel[i].translate(adjustment);
    }
    if (nextVel[i].isZero())  print("drubbin " + i + " has zero velocity");
    if (drubbinPos[i].y <= 7. && drubbinVels[i].y < 0.) {
      // close to ground? stop dropping!
      nextVel[i].y += 0.1;
    }
    drubbinVels[i].normalize();
  }
    
  Point3[] futurePos = (Point3[])drubbinPos.clone();
  for (int i = 0; i < drubbinCount; i++) {
    futurePos[i].translate(drubbinVels[i]);
  }
  for (int i = 0; i < drubbinCount; i++) {
    // Test for proximity of other drubbins & slow down if one is nearby.
    boolean slow = false;
    int collisionCount = 0;
    for (int j = 0; j < drubbinCount; j++) {
      if (i != j) {
        if (pointInSphere(drubbinPos[i], 5., drubbinPos[j])) {
          println("Drubbin " + i + " collided with drubbin " + j); 
          // We should have them fall out of the sky!
          collisionCount++;
        }
        if (pointInSphere(futurePos[i], 5., futurePos[j])) {
          println("Drubbin " + i + " is slowing to avoid hitting drubbin " + j);
          slow = true;
        }
      }
    }
    if (slow) {
      drubbinVels[i].limitMag(0.5);
    }
    if (collisionCount > 0)
      println(collisionCount + " collisions!"); 
  }
  return nextVel;
}

Volume drubbinModel;

Point3 camPos;
Point3 camOrient;

float drubbinTurnRate = fullCircle/80.;

Volume groundX, woodsCircleModel;

Camera cam;

// [setup]
// The method called by Processing to initialize me.
void setup()
{
  size(windowWidth, windowHeight);
  PFont helv = loadFont("HelveticaNeue-13.vlw");
  textFont(helv, 12);
  frameRate(15);
//  origin = new Point3(0,0,0);
  drubbinModel = initWedge();
  drubbinPos = new Point3[drubbinCount];
  for (int i = 0; i < drubbinCount; i++)
    drubbinPos[i] = new Point3((random(1.)-.5) * 90., 0., (random(1.)-.5) * 90.);
  drubbinVels = new Point3[drubbinCount];
  for (int i = 0; i < drubbinCount; i++)
    drubbinVels[i] = new Point3(random(0.2)-0.1, random(0.2)-0.1, random(0.2)-0.1);
  groundX = initGroundX();
  woodsCircleModel = initWoodsCircle();

  cam = new Camera();
}

final float woodsRadius = 100.;
float cameraRange = 600.;

float arrowheadAngle = halfCircle/5;

void drawArrowhead(int x1, int y1, int x2, int y2, int headSize)
{
  line(x1, y1, x2, y2);

  float angle = y1 != y2 ? atan((x2 - x1)/(y2 - y1)) : halfCircle/2;
  float slope1x = sin(angle+arrowheadAngle) * headSize;
  float slope1y = cos(angle+arrowheadAngle) * headSize;
  float slope2x = sin(angle-arrowheadAngle) * headSize;
  float slope2y = cos(angle-arrowheadAngle) * headSize;
  line(x1, y1, x1+slope1x, y1+slope1y);
  line(x1, y1, x1+slope2x, y1+slope2y);
  line(x2, y2, x2-slope1x, y2-slope1y);
  line(x2, y2, x2-slope2x, y2-slope2y);
}


final int drubbinCount = 3;
Point3[] drubbinPos;
Point3[] drubbinVels;



// [draw]
// The method called by Processing to draw one frame.
void draw()
{
  background(200);
  stroke(0);
  
  // Set the camera position/orientation based on mouse position.
  cam.offset = new Point3(Point3.UNITX, halfCircle / 2.0 - mouseY/2.*halfCircle/windowHeight);
  cam.offset.scale(1. + cameraRange *  mouseX/windowWidth);
  cam.orient = direction(cam.offset, origin);
  cam.offset.invert();
  
  // Record previous values of drubbin velocities for later reference
  Point3[] oldDrubbinVels = new Point3[drubbinCount];
  for (int i = 0; i < drubbinCount; i++)
    oldDrubbinVels[i] = drubbinVels[i].clone();
    
  // Determine the drubbins' behavior:
  drubbinVels = drubbinate(drubbinPos, drubbinVels, drubbinCount);

  // For curiosity's sake, measure the acceleration of each drubbin.
  float[] drubbinAccel = new float[drubbinCount];
  float maxDrubbinAccel = 0.;
  for (int i = 0; i < drubbinCount; i++) {
    oldDrubbinVels[i].subtract(drubbinVels[i]);
    drubbinAccel[i] = oldDrubbinVels[i].magnitude();
    if (drubbinAccel[i] > maxDrubbinAccel) maxDrubbinAccel = drubbinAccel[i];
  }
//  print("greatest drubbin acceleration: " + maxDrubbinAccel + "\n");


  // Move the drubbins as per each one's velocity:
  for (int i = 0; i < drubbinCount; i++) {
    drubbinPos[i].translate(drubbinVels[i]);
  }

  // Finally, draw the drubbins:
//  Light lt = new Light(new Point3(100, 100, -100));
//  FlatRenderMethod rm = new FlatRenderMethod(lt);

  WfRenderMethod rm = new WfRenderMethod();
  
  for (int i = 0; i < drubbinCount; i++) {
    Volume drub = drubbinModel.clone();
    drub.linTrans(uprightRotation(drubbinVels[i]));
    drub.translate(drubbinPos[i]);
    renderWithCamera(rm, drub, cam);
  }
  renderWithCamera(rm, groundX, cam);
  renderWithCamera(rm, woodsCircleModel, cam);
  
  // Some helpful text
  stroke(255);
  drawArrowhead(12, windowHeight-8, 25, windowHeight-8, 5);
  drawArrowhead(22, windowHeight-36, 22, windowHeight-22, 5);
  text("track in/out", 32, windowHeight-5);
  text("rotate vertically", 32, windowHeight-25);
}