class Particle {
  Vec3 pos;
  Vec3 vel;
  Vec3 acc;
  boolean collided;
  Particle() {
    pos = new Vec3();
    vel = new Vec3();
    acc = new Vec3();
  }
  Particle(float x, float y, float z) {
    pos = new Vec3(x,y,z);
    vel = new Vec3();
    acc = new Vec3();
  }
  void reset(float x, float y, float z) {
    pos.set(x,y,z);
    vel.set(0f,0f,0f);
    acc.set(0f,0f,0f);
  }
  void integrate(float dt) {
    // apply acceleration to velocity
    vel.x += acc.x * dt;
    vel.y += acc.y * dt;
    vel.z += acc.z * dt;
    // apply velocity to position
    pos.x += vel.x * dt;
    pos.y += vel.y * dt;
    pos.z += vel.z * dt;
    // apply friction to velocity (cheaply)
    float fdt = 1.0 - (friction * dt);
    vel.x *= fdt;
    vel.y *= fdt;
    vel.z *= fdt;
  }
  void draw() {
    pushMatrix();
    translate(pos.x, pos.y, pos.z);
    box(2);
    popMatrix();
  }
}

//        * p1
//       / \
//    v /   \ u
//     /     \
// p3 *-------* p2
//
class Triangle {
  Vec3 p1, p2, p3; 
  Vec3 u, v, uw, vw;
  Vec3 norm;
  float d;
  Triangle(Vec3 p1, Vec3 p2, Vec3 p3) {
    this.p1 = p1.copy();
    this.p2 = p2.copy();
    this.p3 = p3.copy();
    u = p2.sub(p1);
    v = p3.sub(p1);
    this.norm = u.cross(v);
    this.norm.normalize();
    d = -p1.dot(norm);
    Vec3 w = norm;
    uw = u.cross(w);
    vw = v.cross(w);
    float det = 1f / ( w.z*u.x*v.y - w.z*u.y*v.x - u.z*w.x*v.y -
                       u.x*v.z*w.y + v.z*w.x*u.y + u.z*v.x*w.y );
    uw.mulby(-det);    
    vw.mulby(det);
  }
  void draw() {
    beginShape(TRIANGLES);
    vertex(p1.x, p1.y, p1.z);
    vertex(p2.x, p2.y, p2.z);
    vertex(p3.x, p3.y, p3.z);
    //vertex(p1.x, p1.y, p1.z);
    endShape();
  }
}


class Vec3 {
  float x, y, z;
  Vec3() { set(0f,0f,0f); }
  Vec3(float x, float y, float z) { set(x,y,z); }
  Vec3(Vec3 b) { set(b.x,b.y,b.z); }
  Vec3 copy() { return new Vec3(x,y,z); }
  void set(Vec3 b) { x=b.x;  y=b.y;  z=b.z; }
  void set(float x, float y, float z) {
    this.x = x;
    this.y = y;
    this.z = z;
  }
  float length() { return sqrt(x*x+y*y+z*z); }
  float lengthsq() { return x*x+y*y+z*z; }
  void normalize() {
    // assumes length != 0
    float oneoverl = 1.0f / sqrt(x*x+y*y+z*z);
    x *= oneoverl;
    y *= oneoverl;
    z *= oneoverl;
  }
  // the "by" variants perform their work in place
  // the others return new instances (avoid these where possible for performance reasons)
  float dot(Vec3 b) { return x*b.x + y*b.y + z*b.z; }
  Vec3 mul(float s) { return new Vec3(x*s, y*s, z*s); }
  void mulby(float s) { x *= s;  y *= s;  z *= s; }
  Vec3 div(float s) { return new Vec3(x/s, y/s, z/s); }
  void divby(float s) { x /= s;  y /= s;  z /= s; }
  Vec3 add(float s) { return new Vec3(x+s, y+s, z+s); }
  Vec3 add(Vec3 b) { return new Vec3(x+b.x, y+b.y, z+b.z); }
  void addby(float s) { x += s;  y += s;  z += s; }
  void addby(float xx, float yy, float zz) { x+=xx; y+=yy; z+=zz; }
  void addby(Vec3 b) { x+=b.x;  y+=b.y;  z+=b.z; }
  Vec3 sub(float s) { return new Vec3(x-s, y-s, z-s); }
  Vec3 sub(Vec3 b) { return new Vec3(x-b.x, y-b.y, z-b.z); }
  void subby(float s) { x -= s;  y -= s;  z -= s; }
  void subby(Vec3 b) { x-=b.x;  y-=b.y;  z-=b.z; }
  Vec3 cross(Vec3 b) {
    return new Vec3(y*b.z-z*b.y, z*b.x-x*b.z, x*b.y-y*b.x);
  }
  void cross(Vec3 a, Vec3 b) {
    x = a.y*b.z - a.z*b.y;
    y = a.z*b.x - a.x*b.z;
    z = a.x*b.y - a.y*b.x;
  }
}