Move source to src.

src/SDLMain.h

@@ -0,0 +1,11 @@
+/*   SDLMain.m - main entry point for our Cocoa-ized SDL app
+       Initial Version: Darrell Walisser <dwaliss1@purdue.edu>
+       Non-NIB-Code & other changes: Max Horn <max@quendi.de>
+
+    Feel free to customize this file to suit your needs
+*/
+
+#import <Cocoa/Cocoa.h>
+
+@interface SDLMain : NSObject
+@end

src/SDLMain.m

@@ -0,0 +1,382 @@
+/*   SDLMain.m - main entry point for our Cocoa-ized SDL app
+       Initial Version: Darrell Walisser <dwaliss1@purdue.edu>
+       Non-NIB-Code & other changes: Max Horn <max@quendi.de>
+
+    Feel free to customize this file to suit your needs
+*/
+
+#import <SDL/SDL.h>
+#import "SDLMain.h"
+#import <sys/param.h> /* for MAXPATHLEN */
+#import <unistd.h>
+
+/* For some reaon, Apple removed setAppleMenu from the headers in 10.4,
+ but the method still is there and works. To avoid warnings, we declare
+ it ourselves here. */
+@interface NSApplication(SDL_Missing_Methods)
+- (void)setAppleMenu:(NSMenu *)menu;
+@end
+
+/* Use this flag to determine whether we use SDLMain.nib or not */
+#define		SDL_USE_NIB_FILE	0
+
+/* Use this flag to determine whether we use CPS (docking) or not */
+#define		SDL_USE_CPS		1
+#ifdef SDL_USE_CPS
+/* Portions of CPS.h */
+typedef struct CPSProcessSerNum
+{
+	UInt32		lo;
+	UInt32		hi;
+} CPSProcessSerNum;
+
+extern OSErr	CPSGetCurrentProcess( CPSProcessSerNum *psn);
+extern OSErr 	CPSEnableForegroundOperation( CPSProcessSerNum *psn, UInt32 _arg2, UInt32 _arg3, UInt32 _arg4, UInt32 _arg5);
+extern OSErr	CPSSetFrontProcess( CPSProcessSerNum *psn);
+
+#endif /* SDL_USE_CPS */
+
+static int    gArgc;
+static char  **gArgv;
+static BOOL   gFinderLaunch;
+static BOOL   gCalledAppMainline = FALSE;
+
+static NSString *getApplicationName(void)
+{
+    NSDictionary *dict;
+    NSString *appName = 0;
+
+    /* Determine the application name */
+    dict = (NSDictionary *)CFBundleGetInfoDictionary(CFBundleGetMainBundle());
+    if (dict)
+        appName = [dict objectForKey: @"CFBundleName"];
+    
+    if (![appName length])
+        appName = [[NSProcessInfo processInfo] processName];
+
+    return appName;
+}
+
+#if SDL_USE_NIB_FILE
+/* A helper category for NSString */
+@interface NSString (ReplaceSubString)
+- (NSString *)stringByReplacingRange:(NSRange)aRange with:(NSString *)aString;
+@end
+#endif
+
+@interface SDLApplication : NSApplication
+@end
+
+@implementation SDLApplication
+/* Invoked from the Quit menu item */
+- (void)terminate:(id)sender
+{
+    /* Post a SDL_QUIT event */
+    SDL_Event event;
+    event.type = SDL_QUIT;
+    SDL_PushEvent(&event);
+}
+@end
+
+/* The main class of the application, the application's delegate */
+@implementation SDLMain
+
+/* Set the working directory to the .app's parent directory */
+- (void) setupWorkingDirectory:(BOOL)shouldChdir
+{
+    if (shouldChdir)
+    {
+        char appdir[MAXPATHLEN];
+		CFURLRef url = CFBundleCopyBundleURL(CFBundleGetMainBundle());
+		if (CFURLGetFileSystemRepresentation(url, true, (UInt8 *)appdir, MAXPATHLEN)) {
+	        assert ( chdir (appdir) == 0 );   /* chdir to the binary app */
+		}
+		CFRelease(url);
+	}
+
+}
+
+#if SDL_USE_NIB_FILE
+
+/* Fix menu to contain the real app name instead of "SDL App" */
+- (void)fixMenu:(NSMenu *)aMenu withAppName:(NSString *)appName
+{
+    NSRange aRange;
+    NSEnumerator *enumerator;
+    NSMenuItem *menuItem;
+
+    aRange = [[aMenu title] rangeOfString:@"SDL App"];
+    if (aRange.length != 0)
+        [aMenu setTitle: [[aMenu title] stringByReplacingRange:aRange with:appName]];
+
+    enumerator = [[aMenu itemArray] objectEnumerator];
+    while ((menuItem = [enumerator nextObject]))
+    {
+        aRange = [[menuItem title] rangeOfString:@"SDL App"];
+        if (aRange.length != 0)
+            [menuItem setTitle: [[menuItem title] stringByReplacingRange:aRange with:appName]];
+        if ([menuItem hasSubmenu])
+            [self fixMenu:[menuItem submenu] withAppName:appName];
+    }
+    [ aMenu sizeToFit ];
+}
+
+#else
+
+static void setApplicationMenu(void)
+{
+    /* warning: this code is very odd */
+    NSMenu *appleMenu;
+    NSMenuItem *menuItem;
+    NSString *title;
+    NSString *appName;
+    
+    appName = getApplicationName();
+    appleMenu = [[NSMenu alloc] initWithTitle:@""];
+    
+    /* Add menu items */
+    title = [@"About " stringByAppendingString:appName];
+    [appleMenu addItemWithTitle:title action:@selector(orderFrontStandardAboutPanel:) keyEquivalent:@""];
+
+    [appleMenu addItem:[NSMenuItem separatorItem]];
+
+    title = [@"Hide " stringByAppendingString:appName];
+    [appleMenu addItemWithTitle:title action:@selector(hide:) keyEquivalent:@"h"];
+
+    menuItem = (NSMenuItem *)[appleMenu addItemWithTitle:@"Hide Others" action:@selector(hideOtherApplications:) keyEquivalent:@"h"];
+    [menuItem setKeyEquivalentModifierMask:(NSAlternateKeyMask|NSCommandKeyMask)];
+
+    [appleMenu addItemWithTitle:@"Show All" action:@selector(unhideAllApplications:) keyEquivalent:@""];
+
+    [appleMenu addItem:[NSMenuItem separatorItem]];
+
+    title = [@"Quit " stringByAppendingString:appName];
+    [appleMenu addItemWithTitle:title action:@selector(terminate:) keyEquivalent:@"q"];
+
+    
+    /* Put menu into the menubar */
+    menuItem = [[NSMenuItem alloc] initWithTitle:@"" action:nil keyEquivalent:@""];
+    [menuItem setSubmenu:appleMenu];
+    [[NSApp mainMenu] addItem:menuItem];
+
+    /* Tell the application object that this is now the application menu */
+    [NSApp setAppleMenu:appleMenu];
+
+    /* Finally give up our references to the objects */
+    [appleMenu release];
+    [menuItem release];
+}
+
+/* Create a window menu */
+static void setupWindowMenu(void)
+{
+    NSMenu      *windowMenu;
+    NSMenuItem  *windowMenuItem;
+    NSMenuItem  *menuItem;
+
+    windowMenu = [[NSMenu alloc] initWithTitle:@"Window"];
+    
+    /* "Minimize" item */
+    menuItem = [[NSMenuItem alloc] initWithTitle:@"Minimize" action:@selector(performMiniaturize:) keyEquivalent:@"m"];
+    [windowMenu addItem:menuItem];
+    [menuItem release];
+    
+    /* Put menu into the menubar */
+    windowMenuItem = [[NSMenuItem alloc] initWithTitle:@"Window" action:nil keyEquivalent:@""];
+    [windowMenuItem setSubmenu:windowMenu];
+    [[NSApp mainMenu] addItem:windowMenuItem];
+    
+    /* Tell the application object that this is now the window menu */
+    [NSApp setWindowsMenu:windowMenu];
+
+    /* Finally give up our references to the objects */
+    [windowMenu release];
+    [windowMenuItem release];
+}
+
+/* Replacement for NSApplicationMain */
+static void CustomApplicationMain (int argc, char **argv)
+{
+    NSAutoreleasePool	*pool = [[NSAutoreleasePool alloc] init];
+    SDLMain				*sdlMain;
+
+    /* Ensure the application object is initialised */
+    [SDLApplication sharedApplication];
+    
+#ifdef SDL_USE_CPS
+    {
+        CPSProcessSerNum PSN;
+        /* Tell the dock about us */
+        if (!CPSGetCurrentProcess(&PSN))
+            if (!CPSEnableForegroundOperation(&PSN,0x03,0x3C,0x2C,0x1103))
+                if (!CPSSetFrontProcess(&PSN))
+                    [SDLApplication sharedApplication];
+    }
+#endif /* SDL_USE_CPS */
+
+    /* Set up the menubar */
+    [NSApp setMainMenu:[[NSMenu alloc] init]];
+    setApplicationMenu();
+    setupWindowMenu();
+
+    /* Create SDLMain and make it the app delegate */
+    sdlMain = [[SDLMain alloc] init];
+    [NSApp setDelegate:sdlMain];
+    
+    /* Start the main event loop */
+    [NSApp run];
+    
+    [sdlMain release];
+    [pool release];
+}
+
+#endif
+
+
+/*
+ * Catch document open requests...this lets us notice files when the app
+ *  was launched by double-clicking a document, or when a document was
+ *  dragged/dropped on the app's icon. You need to have a
+ *  CFBundleDocumentsType section in your Info.plist to get this message,
+ *  apparently.
+ *
+ * Files are added to gArgv, so to the app, they'll look like command line
+ *  arguments. Previously, apps launched from the finder had nothing but
+ *  an argv[0].
+ *
+ * This message may be received multiple times to open several docs on launch.
+ *
+ * This message is ignored once the app's mainline has been called.
+ */
+- (BOOL)application:(NSApplication *)theApplication openFile:(NSString *)filename
+{
+    const char *temparg;
+    size_t arglen;
+    char *arg;
+    char **newargv;
+
+    if (!gFinderLaunch)  /* MacOS is passing command line args. */
+        return FALSE;
+
+    if (gCalledAppMainline)  /* app has started, ignore this document. */
+        return FALSE;
+
+    temparg = [filename UTF8String];
+    arglen = SDL_strlen(temparg) + 1;
+    arg = (char *) SDL_malloc(arglen);
+    if (arg == NULL)
+        return FALSE;
+
+    newargv = (char **) realloc(gArgv, sizeof (char *) * (gArgc + 2));
+    if (newargv == NULL)
+    {
+        SDL_free(arg);
+        return FALSE;
+    }
+    gArgv = newargv;
+
+    SDL_strlcpy(arg, temparg, arglen);
+    gArgv[gArgc++] = arg;
+    gArgv[gArgc] = NULL;
+    return TRUE;
+}
+
+
+/* Called when the internal event loop has just started running */
+- (void) applicationDidFinishLaunching: (NSNotification *) note
+{
+    int status;
+
+    /* Set the working directory to the .app's parent directory */
+    [self setupWorkingDirectory:gFinderLaunch];
+
+#if SDL_USE_NIB_FILE
+    /* Set the main menu to contain the real app name instead of "SDL App" */
+    [self fixMenu:[NSApp mainMenu] withAppName:getApplicationName()];
+#endif
+
+    /* Hand off to main application code */
+    gCalledAppMainline = TRUE;
+    status = SDL_main (gArgc, gArgv);
+
+    /* We're done, thank you for playing */
+    exit(status);
+}
+@end
+
+
+@implementation NSString (ReplaceSubString)
+
+- (NSString *)stringByReplacingRange:(NSRange)aRange with:(NSString *)aString
+{
+    unsigned int bufferSize;
+    unsigned int selfLen = [self length];
+    unsigned int aStringLen = [aString length];
+    unichar *buffer;
+    NSRange localRange;
+    NSString *result;
+
+    bufferSize = selfLen + aStringLen - aRange.length;
+    buffer = NSAllocateMemoryPages(bufferSize*sizeof(unichar));
+    
+    /* Get first part into buffer */
+    localRange.location = 0;
+    localRange.length = aRange.location;
+    [self getCharacters:buffer range:localRange];
+    
+    /* Get middle part into buffer */
+    localRange.location = 0;
+    localRange.length = aStringLen;
+    [aString getCharacters:(buffer+aRange.location) range:localRange];
+     
+    /* Get last part into buffer */
+    localRange.location = aRange.location + aRange.length;
+    localRange.length = selfLen - localRange.location;
+    [self getCharacters:(buffer+aRange.location+aStringLen) range:localRange];
+    
+    /* Build output string */
+    result = [NSString stringWithCharacters:buffer length:bufferSize];
+    
+    NSDeallocateMemoryPages(buffer, bufferSize);
+    
+    return result;
+}
+
+@end
+
+
+
+#ifdef main
+#  undef main
+#endif
+
+
+/* Main entry point to executable - should *not* be SDL_main! */
+int main (int argc, char **argv)
+{
+    /* Copy the arguments into a global variable */
+    /* This is passed if we are launched by double-clicking */
+    if ( argc >= 2 && strncmp (argv[1], "-psn", 4) == 0 ) {
+        gArgv = (char **) SDL_malloc(sizeof (char *) * 2);
+        gArgv[0] = argv[0];
+        gArgv[1] = NULL;
+        gArgc = 1;
+        gFinderLaunch = YES;
+    } else {
+        int i;
+        gArgc = argc;
+        gArgv = (char **) SDL_malloc(sizeof (char *) * (argc+1));
+        for (i = 0; i <= argc; i++)
+            gArgv[i] = argv[i];
+        gFinderLaunch = NO;
+    }
+
+#if SDL_USE_NIB_FILE
+    [SDLApplication poseAsClass:[NSApplication class]];
+    NSApplicationMain (argc, argv);
+#else
+    CustomApplicationMain (argc, argv);
+#endif
+    return 0;
+}
+

src/ball.cpp

@@ -0,0 +1,26 @@
+// -*- C++ -*-
+
+#include "ball.h"
+#include "material.h"
+#include "model.h"
+#include "physics/shape.h"
+#include "physics/vector.h"
+#include "room.h"
+
+using namespace mbostock;
+
+Ball::Ball(const Vector& x, float radius)
+    : sphere_(x, radius), model_(sphere_) {
+}
+
+Model& Ball::model() {
+  return model_;
+}
+
+const Shape& Ball::shape() const {
+  return sphere_;
+}
+
+float Ball::slip() const {
+  return model_.material().slip();
+}

src/ball.h

@@ -0,0 +1,30 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_BALL_H
+#define MBOSTOCK_BALL_H
+
+#include "model.h"
+#include "physics/shape.h"
+#include "room_object.h"
+
+namespace mbostock {
+
+  class Ball : public RoomObject {
+  public:
+    Ball(const Vector& x, float radius);
+
+    virtual Model& model();
+    virtual const Shape& shape() const;
+    virtual float slip() const;
+
+    inline void setMaterial(const Material& m) { model_.setMaterial(m); }
+
+  private:
+    const Sphere sphere_;
+    SphereModel model_;
+  };
+
+}
+
+
+#endif

src/block.cpp

@@ -0,0 +1,58 @@
+// -*- C++ -*-
+
+#include "block.h"
+#include "material.h"
+#include "model.h"
+#include "physics/shape.h"
+#include "physics/vector.h"
+#include "room.h"
+
+using namespace mbostock;
+
+AxisAlignedBlock::AxisAlignedBlock(const Vector& min, const Vector& max)
+    : box_(min, max), model_(box_) {
+}
+
+Model& AxisAlignedBlock::model() {
+  return model_;
+}
+
+const Shape& AxisAlignedBlock::shape() const {
+  return box_;
+}
+
+void AxisAlignedBlock::setMaterial(const Material& m) {
+  model_.setMaterial(m);
+}
+
+void AxisAlignedBlock::setTopMaterial(const Material& m) {
+  model_.setTopMaterial(m);
+}
+
+float AxisAlignedBlock::slip() const {
+  return model_.material().slip();
+}
+
+Block::Block(const Vector& c, const Vector& x, const Vector& y, const Vector& z)
+    : box_(c, x, y, z), model_(box_) {
+}
+
+Model& Block::model() {
+  return model_;
+}
+
+const Shape& Block::shape() const {
+  return box_;
+}
+
+void Block::setMaterial(const Material& m) {
+  model_.setMaterial(m);
+}
+
+void Block::setTopMaterial(const Material& m) {
+  model_.setTopMaterial(m);
+}
+
+float Block::slip() const {
+  return model_.material().slip();
+}

src/block.h

@@ -0,0 +1,47 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_BLOCK_H
+#define MBOSTOCK_BLOCK_H
+
+#include "model.h"
+#include "physics/shape.h"
+#include "room_object.h"
+
+namespace mbostock {
+
+  class AxisAlignedBlock : public RoomObject {
+  public:
+    AxisAlignedBlock(const Vector& min, const Vector& max);
+
+    virtual Model& model();
+    virtual const Shape& shape() const;
+    virtual float slip() const;
+
+    void setMaterial(const Material& m);
+    void setTopMaterial(const Material& m);
+
+  protected:
+    const AxisAlignedBox box_;
+    AxisAlignedBoxModel model_;
+  };
+
+  class Block : public RoomObject {
+  public:
+    Block(const Vector& c, const Vector& x, const Vector& y, const Vector& z);
+
+    virtual Model& model();
+    virtual const Shape& shape() const;
+    virtual float slip() const;
+
+    void setMaterial(const Material& m);
+    void setTopMaterial(const Material& m);
+
+  private:
+    const Box box_;
+    BoxModel model_;
+  };
+
+}
+
+
+#endif

src/escalator.cpp

@@ -0,0 +1,50 @@
+// -*- C++ -*-
+
+#include "escalator.h"
+#include "material.h"
+#include "model.h"
+#include "physics/particle.h"
+#include "physics/shape.h"
+#include "physics/vector.h"
+
+using namespace mbostock;
+
+Escalator::Escalator(const Vector& min, const Vector& max, const Vector& v)
+    : box_(min, max), velocity_(v * ParticleSimulator::timeStep()),
+      model_(box_) {
+  model_.setTexOrientation((fabsf(v.x) > fabsf(v.z))
+      ? ((v.x > 0) ? AxisAlignedBoxModel::POSITIVE_X
+                   : AxisAlignedBoxModel::NEGATIVE_X)
+      : ((v.z > 0) ? AxisAlignedBoxModel::POSITIVE_Z
+                   : AxisAlignedBoxModel::NEGATIVE_Z));
+}
+
+Model& Escalator::model() {
+  return model_;
+}
+
+const Shape& Escalator::shape() const {
+  return box_;
+}
+
+void Escalator::step(const ParticleSimulator& s) {
+  offset_.x = fmodf(offset_.x + velocity_.x, 1.f);
+  offset_.z = fmodf(offset_.z + velocity_.z, 1.f);
+  model_.setTexOffset(-offset_.x, -offset_.z);
+}
+
+Vector Escalator::velocity(const Vector& x) const {
+  return velocity_;
+}
+
+void Escalator::setMaterial(const Material& m) {
+  model_.setMaterial(m);
+}
+
+void Escalator::setTopMaterial(const Material& m) {
+  model_.setTopMaterial(m);
+}
+
+float Escalator::slip() const {
+  return model_.material().slip();
+}

src/escalator.h

@@ -0,0 +1,39 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_ESCALATOR_H
+#define MBOSTOCK_ESCALATOR_H
+
+#include "model.h"
+#include "physics/shape.h"
+#include "physics/vector.h"
+#include "room_object.h"
+
+namespace mbostock {
+
+  class Material;
+  class ParticleSimulator;
+
+  class Escalator : public DynamicRoomObject {
+  public:
+    Escalator(const Vector& min, const Vector& max, const Vector& v);
+
+    virtual Model& model();
+    virtual const Shape& shape() const;
+    virtual void step(const ParticleSimulator& s);
+    virtual Vector velocity(const Vector& x) const;
+    virtual float slip() const;
+
+    void setMaterial(const Material& m);
+    void setTopMaterial(const Material& m);
+
+  private:
+    const AxisAlignedBox box_;
+    const Vector velocity_;
+    Vector offset_;
+
+    AxisAlignedBoxModel model_;
+  };
+
+}
+
+#endif

src/fan.cpp

@@ -0,0 +1,151 @@
+// -*- C++ -*-
+
+#include <iostream>
+#include <stdlib.h>
+
+#include "fan.h"
+#include "material.h"
+#include "physics/particle.h"
+#include "physics/shape.h"
+
+using namespace mbostock;
+
+namespace mbostock {
+
+  class StaticFanModel : public Model {
+  public:
+    StaticFanModel(float r);
+    virtual ~StaticFanModel();
+
+    virtual void initialize();
+    virtual void display();
+
+    void setMaterial(const Material& m);
+
+  private:
+    void displayBlades();
+
+    const Material* material_;
+    AxisAlignedBox bladeBox_;
+    AxisAlignedBoxModel bladeModel_;
+    GLUquadric* quadric_;
+    float r_;
+  };
+
+}
+
+StaticFanModel::StaticFanModel(float r)
+    : material_(&Materials::blank()),
+      bladeBox_(Vector(0.f, 0.f, 0.f),
+                Vector(r, r / 10.f, r / 20.f)),
+      bladeModel_(bladeBox_), quadric_(NULL), r_(r) {
+}
+
+StaticFanModel::~StaticFanModel() {
+  if (quadric_ != NULL) {
+    gluDeleteQuadric(quadric_);
+  }
+}
+
+void StaticFanModel::setMaterial(const Material& m) {
+  material_ = &m;
+  bladeModel_.setMaterial(m);
+}
+
+void StaticFanModel::initialize() {
+  if (quadric_ != NULL) {
+    gluDeleteQuadric(quadric_);
+  }
+  quadric_ = gluNewQuadric();
+  bladeModel_.initialize();
+}
+
+void StaticFanModel::display() {
+  float axleLength = r_ / 10.f * 1.5;
+  float axleRadius = r_ / 20.f;
+
+  material_->bind();
+  glPushMatrix();
+  glTranslatef(0.f, 0.f, axleLength);
+  gluSphere(quadric_, axleRadius, 16, 16);
+  glPopMatrix();
+  gluCylinder(quadric_, axleRadius, axleRadius, axleLength, 16, 4);
+  gluQuadricOrientation(quadric_, GLU_INSIDE);
+  gluDisk(quadric_, 0.f, axleRadius, 16, 4);
+  gluQuadricOrientation(quadric_, GLU_OUTSIDE);
+  displayBlades();
+}
+
+void StaticFanModel::displayBlades() {
+  int blades = 12;
+  for (int i = 0, n = blades; i < n; i++) {
+    glPushMatrix();
+    glRotatef(i * 360.f / n, 0.f, 0.f, 1.f);
+    glRotatef(60.f, 1.f, 0.f, 0.f);
+    bladeModel_.display();
+    glPopMatrix();
+  }
+}
+
+FanModel::FanModel(const Fan& fan)
+    : fan_(fan), staticModel_(new StaticFanModel(fan.cylinder_.radius())),
+      compiledModel_(Models::compile(staticModel_)) {
+  for (int i = 0; i < 15; i++) {
+    orientation_[i] = 0.f;
+  }
+  orientation_[15] = 1.f;
+}
+
+FanModel::~FanModel() {
+  delete compiledModel_;
+}
+
+float* FanModel::orientation() {
+  const Vector& z = fan_.cylinder_.z();
+  const Vector& n = (fabsf(z.z) > .5f) ? Vector::Y() : Vector::Z();
+  const Vector& x = n.cross(z);
+  const Vector& y = -x.cross(z);
+  orientation_[0] = x.x; orientation_[1] = x.y; orientation_[2] = x.z;
+  orientation_[4] = y.x; orientation_[5] = y.y; orientation_[6] = y.z;
+  orientation_[8] = z.x; orientation_[9] = z.y; orientation_[10] = z.z;
+  return orientation_;
+}
+
+void FanModel::setMaterial(const Material& m) {
+  staticModel_->setMaterial(m);
+}
+
+void FanModel::initialize() {
+  compiledModel_->initialize();
+}
+
+void FanModel::display() {
+  glPushMatrix();
+  glTranslatev(fan_.cylinder_.x0());
+  glMultMatrixf(orientation());
+  glRotatef(fan_.a_, 0.f, 0.f, 1.f);
+  compiledModel_->display();
+  glPopMatrix();
+}
+
+Fan::Fan(const Vector& x, const Vector& v, float r, float s)
+    : cylinder_(x, x + v * (r / 10.f), r),
+      s_(s * ParticleSimulator::timeStep()),
+      a_(0.f), model_(*this) {
+}
+
+Model& Fan::model() {
+  return model_;
+}
+
+const Shape& Fan::shape() const {
+  return cylinder_;
+}
+
+void Fan::step(const ParticleSimulator& s) {
+  a_ += s_;
+}
+
+void Fan::reset() {
+  a_ = 0.f;
+}

src/fan.h

@@ -0,0 +1,55 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_FAN_H
+#define MBOSTOCK_FAN_H
+
+#include "model.h"
+#include "physics/shape.h"
+#include "room_object.h"
+
+namespace mbostock {
+
+  class Fan;
+  class StaticFanModel;
+
+  class FanModel : public Model {
+  public:
+    FanModel(const Fan& fan);
+    virtual ~FanModel();
+
+    virtual void initialize();
+    virtual void display();
+
+    void setMaterial(const Material& m);
+
+  private:
+    float* orientation();
+
+    const Fan& fan_;
+    StaticFanModel* staticModel_;
+    Model* compiledModel_;
+    float orientation_[16];
+  };
+
+  class Fan : public DynamicRoomObject {
+  public:
+    Fan(const Vector& x, const Vector& v, float r, float s);
+
+    virtual Model& model();
+    virtual const Shape& shape() const;
+    virtual void step(const ParticleSimulator& s);
+    virtual void reset();
+
+    inline void setMaterial(const Material& m) { model_.setMaterial(m); }
+
+  private:
+    const Cylinder cylinder_;
+    const float s_;
+    float a_;
+    FanModel model_;
+
+    friend class FanModel;
+  };
+}
+
+#endif

src/lighting.cpp

@@ -0,0 +1,144 @@
+// -*- C++ -*-
+
+#include <stdlib.h>
+
+#include "lighting.h"
+
+using namespace mbostock;
+
+Light::Light()
+  : id_(-1), enabled_(false) {
+  setAmbient(0.f, 0.f, 0.f, 1.f);
+  setDiffuse(0.f, 0.f, 0.f, 1.f);
+  setSpecular(0.f, 0.f, 0.f, 1.f);
+  setPosition(0.f, 0.f, 1.f, 0.f);
+  setSpotDirection(0.f, 0.f, -1.f);
+  spotExponent_ = 0.f;
+  spotCutoff_ = 180.f;
+  constantAttenuation_ = 1.f;
+  linearAttenuation_ = 0.f;
+  quadraticAttenuation_ = 0.f;
+}
+
+void Light::setAmbient(float r, float g, float b, float a) {
+  ambient_[0] = r;
+  ambient_[1] = g;
+  ambient_[2] = b;
+  ambient_[3] = a;
+}
+
+void Light::setDiffuse(float r, float g, float b, float a) {
+  diffuse_[0] = r;
+  diffuse_[1] = g;
+  diffuse_[2] = b;
+  diffuse_[3] = a;
+}
+
+void Light::setSpecular(float r, float g, float b, float a) {
+  specular_[0] = r;
+  specular_[1] = g;
+  specular_[2] = b;
+  specular_[3] = a;
+}
+
+void Light::setPosition(float x, float y, float z, float w) {
+  position_[0] = x;
+  position_[1] = y;
+  position_[2] = z;
+  position_[3] = w;
+}
+
+void Light::setSpotDirection(float x, float y, float z) {
+  spotDirection_[0] = x;
+  spotDirection_[1] = y;
+  spotDirection_[2] = z;
+}
+
+void Light::setSpotExponent(float e) {
+  spotExponent_ = e;
+}
+
+void Light::setConstantAttenuation(float a) {
+  constantAttenuation_ = a;
+}
+
+void Light::setLinearAttenuation(float a) {
+  linearAttenuation_ = a;
+}
+
+void Light::setQuadraticAttenuation(float a) {
+  quadraticAttenuation_ = a;
+}
+
+void Light::initialize() const {
+  if (enabled_) {
+    glEnable(id_);
+    glLightfv(id_, GL_AMBIENT, ambient_);
+    glLightfv(id_, GL_DIFFUSE, diffuse_);
+    glLightfv(id_, GL_SPECULAR, specular_);
+    glLightf(id_, GL_SPOT_EXPONENT, spotExponent_);
+    glLightf(id_, GL_SPOT_CUTOFF, spotCutoff_);
+    glLightf(id_, GL_CONSTANT_ATTENUATION, constantAttenuation_);
+    glLightf(id_, GL_LINEAR_ATTENUATION, linearAttenuation_);
+    glLightf(id_, GL_QUADRATIC_ATTENUATION, quadraticAttenuation_);
+  } else {
+    glDisable(id_);
+  }
+}
+
+void Light::display() const {
+  if (enabled_) {
+    glLightfv(id_, GL_POSITION, position_);
+    glLightfv(id_, GL_SPOT_DIRECTION, spotDirection_);
+  }
+}
+
+void Light::enable() {
+  enabled_ = true;
+}
+
+void Light::disable() {
+  enabled_ = false;
+}
+
+Lighting::Lighting() {
+  for (int i = 0; i < lights(); i++) {
+    lights_[i].id_ = GL_LIGHT0 + i;
+  }
+  setGlobalAmbient(.2f, .2f, .2f, 1.f);
+  lights_[0].setDiffuse(1.f, 1.f, 1.f, 1.f);
+  lights_[0].setSpecular(1.f, 1.f, 1.f, 1.f);
+  lights_[0].enable();
+}
+
+void Lighting::setGlobalAmbient(float r, float g, float b, float a) {
+  globalAmbient_[0] = r;
+  globalAmbient_[1] = g;
+  globalAmbient_[2] = b;
+  globalAmbient_[3] = a;
+}
+
+static const Lighting* lastLighting_ = NULL;
+
+void Lighting::initialize() const {
+  lastLighting_ = this;
+  glLightModelfv(GL_LIGHT_MODEL_AMBIENT, globalAmbient_);
+  for (int i = 0; i < lights(); i++) {
+    lights_[i].initialize();
+  }
+}
+
+void Lighting::display() const {
+  if (lastLighting_ != this) {
+    lastLighting_ = this;
+    initialize();
+  }
+  for (int i = 0; i < lights(); i++) {
+    lights_[i].display();
+  }
+}
+
+const Lighting& Lightings::standard() {
+  static const Lighting l;
+  return l;
+}

src/lighting.h

@@ -0,0 +1,74 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_LIGHTING_H
+#define MBOSTOCK_LIGHTING_H
+
+#include <OpenGL/gl.h>
+
+namespace mbostock {
+
+  class Light {
+  public:
+    Light();
+
+    void initialize() const;
+    void display() const;
+
+    void enable();
+    void disable();
+    inline bool enabled() const { return enabled_; }
+
+    void setAmbient(float r, float g, float b, float a);
+    void setDiffuse(float r, float g, float b, float a);
+    void setSpecular(float r, float g, float b, float a);
+    void setPosition(float x, float y, float z, float w);
+    void setSpotDirection(float x, float y, float z);
+    void setSpotExponent(float e);
+    void setConstantAttenuation(float a);
+    void setLinearAttenuation(float a);
+    void setQuadraticAttenuation(float a);
+
+  private:
+    GLenum id_;
+    bool enabled_;
+    float ambient_[4];
+    float diffuse_[4];
+    float specular_[4];
+    float position_[4];
+    float spotDirection_[3];
+    float spotExponent_;
+    float spotCutoff_;
+    float constantAttenuation_;
+    float linearAttenuation_;
+    float quadraticAttenuation_;
+
+    friend class Lighting;
+  };
+
+  class Lighting {
+  public:
+    Lighting();
+
+    void initialize() const;
+    void display() const;
+
+    void setGlobalAmbient(float r, float g, float b, float a);
+    inline Light& light(int i) { return lights_[i]; }
+    inline int lights() const { return 8; }
+
+  private:
+    float globalAmbient_[4];
+    Light lights_[8];
+  };
+
+  class Lightings {
+  public:
+    static const Lighting& standard();
+
+  private:
+    Lightings();
+  };
+
+}
+
+#endif

src/main.cpp

@@ -0,0 +1,215 @@
+// -*- C++ -*-
+
+#include <OpenGL/gl.h>
+#include <OpenGL/glu.h>
+#include <SDL/SDL.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <TinyXML/tinyxml.h>
+
+#include "room.h"
+#include "shader.h"
+#include "sound.h"
+#include "texture.h"
+#include "world.h"
+#include "worlds.h"
+
+using namespace mbostock;
+
+static const int defaultWidth = 640;
+static const int defaultHeight = 480;
+static int screenWidth = 0;
+static int screenHeight = 0;
+static const int defaultX = 50;
+static const int defaultY = 50;
+static const float kd = .060f; // frame-rate dependent
+
+static bool run = true;
+static bool fullScreen = false;
+
+static World* world = NULL;
+static bool wireframe = false;
+
+static Shader* shaders[] = {
+  Shaders::defaultShader(),
+  Shaders::wireframeShader(),
+  Shaders::normalShader()
+};
+
+static int shaderi = 0;
+static const int shadern = 3;
+
+static Shader* shader() {
+  return shaders[shaderi];
+}
+
+static void resizeSurface(int width, int height) {
+  uint32_t flags = SDL_OPENGL | SDL_RESIZABLE;
+  if (fullScreen) {
+    flags |= SDL_FULLSCREEN;
+  }
+  SDL_SetVideoMode(width, height, 24, flags);
+
+  /* Store the screen resolution when going into full screen. */
+  if (width == 0) {
+    const SDL_VideoInfo* info = SDL_GetVideoInfo();
+    width = screenWidth = info->current_w;
+    height = screenHeight = info->current_h;
+  }
+
+  glViewport(0, 0, width, height);
+  glMatrixMode(GL_PROJECTION);
+  glLoadIdentity();
+  gluPerspective(45.f, width / (float) height, 1.0f, 100.f);
+  glMatrixMode(GL_MODELVIEW);
+  glClearColor(0.f, 0.f, 0.f, 0.f);
+
+  shader()->initialize();
+  Textures::initialize();
+  world->model().initialize();
+}
+
+static void handleDisplay() {
+  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+  glLoadIdentity();
+
+  world->simulate();
+
+  const Vector& p = world->player().origin();
+  const Vector& min = world->room().cameraBounds().min();
+  const Vector& max = world->room().cameraBounds().max();
+
+  /* Interpolate the eye location. */
+  Vector ee(p.x, p.y + 4.f, p.z + 6.f);
+  ee = Vector::min(Vector::max(min, ee), max);
+  static Vector e = ee;
+  e = e * (1.f - kd) + ee * kd;
+
+  /* Interpolate the camera direction towards the player. */
+  static Vector c = p;
+  c = c * (1.f - kd) + p * kd;
+
+  gluLookAt(e.x, e.y, e.z,
+            c.x, c.y, c.z,
+            0.f, 1.f, 0.f);
+
+  shader()->display(world->model());
+  SDL_GL_SwapBuffers();
+}
+
+static void toggleShader() {
+  shaderi = (shaderi + 1) % shadern;
+  shader()->initialize();
+}
+
+static void toggleFullScreen() {
+  fullScreen = !fullScreen;
+  if (fullScreen) {
+    resizeSurface(screenWidth, screenHeight);
+    SDL_ShowCursor(SDL_DISABLE);
+  } else {
+    resizeSurface(defaultWidth, defaultHeight);
+    SDL_ShowCursor(SDL_ENABLE);
+  }
+}
+
+static void handleKeyDown(SDL_Event* event) {
+  switch (event->key.keysym.sym) {
+    case SDLK_LEFT: {
+      if (event->key.keysym.mod & KMOD_META) {
+        world->previousRoom();
+      }
+      break;
+    }
+    case SDLK_DOWN: {
+      if (event->key.keysym.mod & KMOD_META) {
+        world->resetPlayer();
+      }
+      break;
+    }
+    case SDLK_RIGHT: {
+      if (event->key.keysym.mod & KMOD_META) {
+        world->nextRoom();
+      }
+      break;
+    }
+    case SDLK_a: world->player().move(Player::LEFT); break;
+    case SDLK_s: world->player().move(Player::BACKWARD); break;
+    case SDLK_d: world->player().move(Player::RIGHT); break;
+    case SDLK_w: world->player().move(Player::FORWARD); break;
+  }
+}
+
+static void handleKeyUp(SDL_Event* event) {
+  switch (event->key.keysym.sym) {
+    case SDLK_a: world->player().stop(Player::LEFT); break;
+    case SDLK_s: world->player().stop(Player::BACKWARD); break;
+    case SDLK_d: world->player().stop(Player::RIGHT); break;
+    case SDLK_w: world->player().stop(Player::FORWARD); break;
+    case SDLK_SPACE: world->togglePaused(); break;
+    case SDLK_q: if (!(event->key.keysym.mod & KMOD_META)) break;
+    case SDLK_ESCAPE: run = false; break;
+    case SDLK_F9: toggleShader(); break;
+    case SDLK_F10: world->toggleDebug(); break;
+    case SDLK_F11: toggleFullScreen(); break;
+  }
+}
+
+static void handleQuit() {
+  Sounds::dispose();
+  delete world;
+  SDL_Quit();
+}
+
+static void eventLoop() {
+  SDL_Event event;
+  while (run) {
+    handleDisplay();
+    if (SDL_PollEvent(&event)) {
+      switch (event.type) {
+        case SDL_VIDEORESIZE: {
+          resizeSurface(event.resize.w, event.resize.h);
+          break;
+        }
+        case SDL_KEYDOWN: {
+          handleKeyDown(&event);
+          break;
+        }
+        case SDL_KEYUP: {
+          handleKeyUp(&event);
+          break;
+        }
+        case SDL_QUIT: {
+          run = false;
+          break;
+        }
+      }
+    }
+    SDL_Delay(10);
+  }
+  handleQuit();
+  return;
+}
+
+int main(int argc, char** argv) {
+  SDL_Init(SDL_INIT_VIDEO | SDL_INIT_AUDIO);
+
+  SDL_GL_SetAttribute(SDL_GL_SWAP_CONTROL, 1);
+  SDL_GL_SetAttribute(SDL_GL_DOUBLEBUFFER, 1);
+  SDL_GL_SetAttribute(SDL_GL_DEPTH_SIZE, 24);
+  SDL_GL_SetAttribute(SDL_GL_RED_SIZE, 8);
+  SDL_GL_SetAttribute(SDL_GL_GREEN_SIZE, 8);
+  SDL_GL_SetAttribute(SDL_GL_BLUE_SIZE, 8);
+  SDL_GL_SetAttribute(SDL_GL_ALPHA_SIZE, 8);
+  SDL_GL_SetAttribute(SDL_GL_MULTISAMPLEBUFFERS, 1);
+  SDL_GL_SetAttribute(SDL_GL_MULTISAMPLESAMPLES, 4);
+  SDL_WM_SetCaption("POLLY-B-GONE", "POLLY-B-GONE");
+
+  Sounds::initialize();
+  world = Worlds::fromFile("world.xml");
+  // resizeSurface(defaultWidth, defaultHeight);
+  toggleFullScreen();
+  eventLoop();
+
+  return 0;
+}

src/material.cpp

@@ -0,0 +1,88 @@
+// -*- C++ -*-
+
+#include <OpenGL/gl.h>
+#include <iostream>
+#include <math.h>
+#include <stdlib.h>
+
+#include "material.h"
+#include "texture.h"
+
+using namespace mbostock;
+
+Material::Material()
+    : shininess_(0.f), texture_(NULL), slip_(0.f) {
+  for (int i = 0; i < 4; i++) {
+    ambient_[i] = 0.f;
+    diffuse_[i] = 0.f;
+    specular_[i] = 0.f;
+    emission_[i] = 0.f;
+  }
+}
+
+void Material::setAmbient(float r, float g, float b, float a) {
+  ambient_[0] = r;
+  ambient_[1] = g;
+  ambient_[2] = b;
+  ambient_[3] = a;
+}
+
+void Material::setDiffuse(float r, float g, float b, float a) {
+  diffuse_[0] = r;
+  diffuse_[1] = g;
+  diffuse_[2] = b;
+  diffuse_[3] = a;
+}
+
+void Material::setSpecular(float r, float g, float b, float a) {
+  specular_[0] = r;
+  specular_[1] = g;
+  specular_[2] = b;
+  specular_[3] = a;
+}
+
+void Material::setEmission(float r, float g, float b, float a) {
+  emission_[0] = r;
+  emission_[1] = g;
+  emission_[2] = b;
+  emission_[3] = a;
+}
+
+void Material::setShininess(float s) {
+  shininess_ = s;
+}
+
+void Material::setTexture(const char* path) {
+  texture_ = (path == NULL) ? NULL : &(Textures::fromFile(path));
+}
+
+void Material::setSlipAngle(float angle) {
+  slip_ = cosf(angle * (2.f * M_PI / 360.f));
+}
+
+void Material::bind() const {
+  glMaterialfv(GL_FRONT, GL_AMBIENT, ambient_);
+  glMaterialfv(GL_FRONT, GL_DIFFUSE, diffuse_);
+  glMaterialfv(GL_FRONT, GL_SPECULAR, specular_);
+  glMaterialfv(GL_FRONT, GL_EMISSION, emission_);
+  glMaterialf(GL_FRONT, GL_SHININESS, shininess_);
+  if (texture_ == NULL) {
+    glDisable(GL_BLEND);
+    glBindTexture(GL_TEXTURE_2D, GL_NONE);
+  } else {
+    texture_->bind();
+  }
+}
+
+const Material& Materials::blank() {
+  static Material m;
+  static bool init = false;
+  if (!init) {
+    init = true;
+    m.setAmbient(.2f, .2f, .1f, .1f);
+    m.setDiffuse(.4f, .1f, .1f, 1.f);
+    m.setSpecular(.05f, .05f, .05f, .1f);
+    m.setShininess(10.f);
+  }
+  return m;
+}

src/material.h

@@ -0,0 +1,47 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_MATERIAL_H
+#define MBOSTOCK_MATERIAL_H
+
+namespace mbostock {
+
+  class Texture;
+
+  class Material {
+  public:
+    Material();
+
+    void setAmbient(float r, float g, float b, float a);
+    void setDiffuse(float r, float g, float b, float a);
+    void setSpecular(float r, float g, float b, float a);
+    void setEmission(float r, float g, float b, float a);
+    void setShininess(float s);
+    void setTexture(const char* path);
+
+    void setSlipAngle(float angle);
+    inline float slip() const { return slip_; }
+
+    void bind() const;
+
+  private:
+    float ambient_[4];
+    float diffuse_[4];
+    float specular_[4];
+    float emission_[4];
+    float shininess_;
+    const Texture* texture_;
+
+    float slip_;
+  };
+
+  class Materials {
+  public:
+    static const Material& blank();
+
+  private:
+    Materials();
+  };
+
+}
+
+#endif

src/model.cpp

@@ -0,0 +1,607 @@
+// -*- C++ -*-
+
+#include <GLUT/glut.h>
+#include <OpenGL/glu.h>
+#include <algorithm>
+#include <iostream>
+#include <math.h>
+#include <stdlib.h>
+
+#include "material.h"
+#include "model.h"
+
+using namespace mbostock;
+
+WedgeModel::WedgeModel(const Wedge& wedge)
+    : wedge_(wedge), material_(&Materials::blank()), topMaterial_(NULL) {
+}
+
+void WedgeModel::setMaterial(const Material& m) {
+  material_ = &m;
+}
+
+void WedgeModel::setTopMaterial(const Material& m) {
+  topMaterial_ = &m;
+}
+
+const Material& WedgeModel::topMaterial() const {
+  return (topMaterial_ == NULL) ? *material_ : *topMaterial_;
+}
+
+void WedgeModel::display() {
+  Vector size = wedge_.x2() - wedge_.x0();
+
+  /* Top. */
+  ((topMaterial_ != NULL) ? topMaterial_ : material_)->bind();
+  glBegin(GL_QUADS);
+  glNormalv(wedge_.top().normal());
+  glTexCoord2f(0, 0);
+  glVertexv(wedge_.x0());
+  glTexCoord2f(sqrtf(size.x * size.x + size.y * size.y), 0);
+  glVertexv(wedge_.x1());
+  glTexCoord2f(sqrtf(size.x * size.x + size.y * size.y), size.z);
+  glVertexv(wedge_.x2());
+  glTexCoord2f(0, size.z);
+  glVertexv(wedge_.x3());
+  glEnd();
+
+  material_->bind();
+  glBegin(GL_QUADS);
+
+  /* Right. */
+  glNormalv(wedge_.right().normal());
+  glTexCoord2f(size.z, size.y);
+  glVertexv(wedge_.x2());
+  glTexCoord2f(0, size.y);
+  glVertexv(wedge_.x1());
+  glTexCoord2f(0, 0);
+  glVertexv(wedge_.x4());
+  glTexCoord2f(size.z, 0);
+  glVertexv(wedge_.x5());
+
+  /* Bottom. */
+  glNormalv(wedge_.bottom().normal());
+  glTexCoord2f(0, 0);
+  glVertexv(wedge_.x0());
+  glTexCoord2f(0, size.z);
+  glVertexv(wedge_.x3());
+  glTexCoord2f(size.x, size.z);
+  glVertexv(wedge_.x5());
+  glTexCoord2f(size.x, 0);
+  glVertexv(wedge_.x4());
+
+  glEnd();
+  glBegin(GL_TRIANGLES);
+
+  /* Front. */
+  glNormalv(wedge_.front().normal());
+  glTexCoord2f(size.x, size.y);
+  glVertexv(wedge_.x1());
+  glTexCoord2f(0, 0);
+  glVertexv(wedge_.x0());
+  glTexCoord2f(size.x, 0);
+  glVertexv(wedge_.x4());
+
+  /* Back. */
+  glNormalv(wedge_.back().normal());
+  glTexCoord2f(0, 0);
+  glVertexv(wedge_.x3());
+  glTexCoord2f(size.x, size.y);
+  glVertexv(wedge_.x2());
+  glTexCoord2f(size.x, 0);
+  glVertexv(wedge_.x5());
+
+  glEnd();
+}
+
+AxisAlignedBoxModel::AxisAlignedBoxModel(const AxisAlignedBox& box)
+    : box_(box), material_(&Materials::blank()), topMaterial_(NULL),
+      orientation_(POSITIVE_X), u_(0.f), v_(0.f) {
+}
+
+void AxisAlignedBoxModel::setMaterial(const Material& m) {
+  material_ = &m;
+}
+
+void AxisAlignedBoxModel::setTopMaterial(const Material& m) {
+  topMaterial_ = &m;
+}
+
+const Material& AxisAlignedBoxModel::topMaterial() const {
+  return (topMaterial_ == NULL) ? *material_ : *topMaterial_;
+}
+
+void AxisAlignedBoxModel::display() {
+  Vector size = box_.x7() - box_.x0();
+  Vector t4, t5, t6, t7;
+  switch (orientation_) {
+    case POSITIVE_X: {
+      t4 = Vector(size.x + u_, 0, v_);
+      t5 = Vector(u_, 0, v_);
+      t6 = Vector(u_, 0, size.z + v_);
+      t7 = Vector(size.x + u_, 0, size.z + v_);
+      break;
+    }
+    case POSITIVE_Z: {
+      t5 = Vector(v_, 0, u_);
+      t6 = Vector(size.z + v_, 0, u_);
+      t7 = Vector(size.z + v_, 0, size.x + u_);
+      t4 = Vector(v_, 0, size.x + u_);
+      break;
+    }
+    case NEGATIVE_X: {
+      t4 = Vector(-u_, 0, v_);
+      t5 = Vector(size.x - u_, 0, v_);
+      t6 = Vector(size.x - u_, 0, size.z + v_);
+      t7 = Vector(-u_, 0, size.z + v_);
+      break;
+    }
+    case NEGATIVE_Z: {
+      t5 = Vector(size.z - v_, 0, -u_);
+      t6 = Vector(-v_, 0, -u_);
+      t7 = Vector(-v_, 0, size.x - u_);
+      t4 = Vector(size.z - v_, 0, size.x - u_);
+      break;
+    }
+  }
+
+  Vector y(0, size.y, 0);
+  Vector t0 = t5 + y;
+  Vector t1 = t4 + y;
+  Vector t2 = t7 + y;
+  Vector t3 = t6 + y;
+
+  /* Top. */
+  ((topMaterial_ != NULL) ? topMaterial_ : material_)->bind();
+  glBegin(GL_QUADS);
+  glNormal3f(0.f, 1.f, 0.f);
+  glTexCoord2f(t4.x, t4.z);
+  glVertexv(box_.x4());
+  glTexCoord2f(t5.x, t5.z);
+  glVertexv(box_.x5());
+  glTexCoord2f(t6.x, t6.z);
+  glVertexv(box_.x6());
+  glTexCoord2f(t7.x, t7.z);
+  glVertexv(box_.x7());
+  glEnd();
+
+  material_->bind();
+  glBegin(GL_QUADS);
+
+  /* Bottom. */
+  glNormal3f(0.f, -1.f, 0.f);
+  glTexCoord2f(t0.x, t0.z);
+  glVertexv(box_.x0());
+  glTexCoord2f(t1.x, t1.z);
+  glVertexv(box_.x1());
+  glTexCoord2f(t2.x, t2.z);
+  glVertexv(box_.x2());
+  glTexCoord2f(t3.x, t3.z);
+  glVertexv(box_.x3());
+
+  /* Front. */
+  glNormal3f(0.f, 0.f, 1.f);
+  glTexCoord2f(t7.x, t7.y);
+  glVertexv(box_.x7());
+  glTexCoord2f(t6.x, t6.y);
+  glVertexv(box_.x6());
+  glTexCoord2f(t3.x, t3.y);
+  glVertexv(box_.x3());
+  glTexCoord2f(t2.x, t2.y);
+  glVertexv(box_.x2());
+
+  /* Back. */
+  glNormal3f(0.f, 0.f, -1.f);
+  glTexCoord2f(t5.x, t5.y);
+  glVertexv(box_.x5());
+  glTexCoord2f(t4.x, t4.y);
+  glVertexv(box_.x4());
+  glTexCoord2f(t1.x, t1.y);
+  glVertexv(box_.x1());
+  glTexCoord2f(t0.x, t0.y);
+  glVertexv(box_.x0());
+
+  /* Left. */
+  glNormal3f(-1.f, 0.f, 0.f);
+  glTexCoord2f(t6.z, t6.y);
+  glVertexv(box_.x6());
+  glTexCoord2f(t5.z, t5.y);
+  glVertexv(box_.x5());
+  glTexCoord2f(t0.z, t0.y);
+  glVertexv(box_.x0());
+  glTexCoord2f(t3.z, t3.y);
+  glVertexv(box_.x3());
+
+  /* Right. */
+  glNormal3f(1.f, 0.f, 0.f);
+  glTexCoord2f(t4.z, t4.y);
+  glVertexv(box_.x4());
+  glTexCoord2f(t7.z, t7.y);
+  glVertexv(box_.x7());
+  glTexCoord2f(t2.z, t2.y);
+  glVertexv(box_.x2());
+  glTexCoord2f(t1.z, t1.y);
+  glVertexv(box_.x1());
+
+  glEnd();
+}
+
+void AxisAlignedBoxModel::setTexOffset(float u, float v) {
+  u_ = u;
+  v_ = v;
+}
+
+void AxisAlignedBoxModel::setTexOrientation(Orientation orientation) {
+  orientation_ = orientation;
+}
+
+BoxModel::BoxModel(const Box& box)
+    : box_(box), material_(&Materials::blank()), topMaterial_(NULL) {
+}
+
+void BoxModel::setMaterial(const Material& m) {
+  material_ = &m;
+}
+
+void BoxModel::setTopMaterial(const Material& m) {
+  topMaterial_ = &m;
+}
+
+const Material& BoxModel::topMaterial() const {
+  return (topMaterial_ == NULL) ? *material_ : *topMaterial_;
+}
+
+void BoxModel::display() {
+  Vector size(
+      (box_.x1() - box_.x0()).length(),
+      (box_.x5() - box_.x0()).length(),
+      (box_.x3() - box_.x0()).length());
+
+  Vector t4(size.x, 0, 0);
+  Vector t5(0, 0, 0);
+  Vector t6(0, 0, size.z);
+  Vector t7(size.x, 0, size.z);
+
+  Vector y(0, size.y, 0);
+  Vector t0 = t5 + y;
+  Vector t1 = t4 + y;
+  Vector t2 = t7 + y;
+  Vector t3 = t6 + y;
+
+  /* Top. */
+  ((topMaterial_ != NULL) ? topMaterial_ : material_)->bind();
+  glBegin(GL_QUADS);
+  glNormalv(box_.top().normal());
+  glTexCoord2f(t4.x, t4.z);
+  glVertexv(box_.x4());
+  glTexCoord2f(t5.x, t5.z);
+  glVertexv(box_.x5());
+  glTexCoord2f(t6.x, t6.z);
+  glVertexv(box_.x6());
+  glTexCoord2f(t7.x, t7.z);
+  glVertexv(box_.x7());
+  glEnd();
+
+  material_->bind();
+  glBegin(GL_QUADS);
+
+  /* Bottom. */
+  glNormalv(box_.bottom().normal());
+  glTexCoord2f(t0.x, t0.z);
+  glVertexv(box_.x0());
+  glTexCoord2f(t1.x, t1.z);
+  glVertexv(box_.x1());
+  glTexCoord2f(t2.x, t2.z);
+  glVertexv(box_.x2());
+  glTexCoord2f(t3.x, t3.z);
+  glVertexv(box_.x3());
+
+  /* Front. */
+  glNormalv(box_.front().normal());
+  glTexCoord2f(t7.x, t7.y);
+  glVertexv(box_.x7());
+  glTexCoord2f(t6.x, t6.y);
+  glVertexv(box_.x6());
+  glTexCoord2f(t3.x, t3.y);
+  glVertexv(box_.x3());
+  glTexCoord2f(t2.x, t2.y);
+  glVertexv(box_.x2());
+
+  /* Back. */
+  glNormalv(box_.back().normal());
+  glTexCoord2f(t5.x, t5.y);
+  glVertexv(box_.x5());
+  glTexCoord2f(t4.x, t4.y);
+  glVertexv(box_.x4());
+  glTexCoord2f(t1.x, t1.y);
+  glVertexv(box_.x1());
+  glTexCoord2f(t0.x, t0.y);
+  glVertexv(box_.x0());
+
+  /* Left. */
+  glNormalv(box_.left().normal());
+  glTexCoord2f(t6.z, t6.y);
+  glVertexv(box_.x6());
+  glTexCoord2f(t5.z, t5.y);
+  glVertexv(box_.x5());
+  glTexCoord2f(t0.z, t0.y);
+  glVertexv(box_.x0());
+  glTexCoord2f(t3.z, t3.y);
+  glVertexv(box_.x3());
+
+  /* Right. */
+  glNormalv(box_.right().normal());
+  glTexCoord2f(t4.z, t4.y);
+  glVertexv(box_.x4());
+  glTexCoord2f(t7.z, t7.y);
+  glVertexv(box_.x7());
+  glTexCoord2f(t2.z, t2.y);
+  glVertexv(box_.x2());
+  glTexCoord2f(t1.z, t1.y);
+  glVertexv(box_.x1());
+
+  glEnd();
+}
+
+QuadModel::QuadModel(const Quad& quad)
+    : quad_(quad), material_(&Materials::blank()) {
+  float x = (quad_.x1() - quad_.x0()).length();
+  float y = (quad_.x3() - quad_.x0()).length();
+  if (x < y) {
+    texCoords_[0] = Vector(0, 0);
+    texCoords_[1] = Vector(0, x);
+    texCoords_[2] = Vector(y, x);
+    texCoords_[3] = Vector(y, 0);
+  } else {
+    texCoords_[0] = Vector(0, 0);
+    texCoords_[1] = Vector(x, 0);
+    texCoords_[2] = Vector(x, y);
+    texCoords_[3] = Vector(0, y);
+  }
+}
+
+void QuadModel::setMaterial(const Material& m) {
+  material_ = &m;
+}
+
+void QuadModel::setTexCoords(const Vector& t0, const Vector& t1,
+                             const Vector& t2, const Vector& t3) {
+  texCoords_[0] = t0;
+  texCoords_[1] = t1;
+  texCoords_[2] = t2;
+  texCoords_[3] = t3;
+}
+
+void QuadModel::display() {
+  material_->bind();
+  displaySide(true);
+  glFrontFace(GL_CW);
+  displaySide(false);
+  glFrontFace(GL_CCW);
+}
+
+void QuadModel::displaySide(bool front) {
+  glBegin(GL_QUADS);
+  glNormalv(front ? quad_.normal() : -quad_.normal());
+  glTexCoordv(texCoords_[0]);
+  glVertexv(quad_.x0());
+  glTexCoordv(texCoords_[1]);
+  glVertexv(quad_.x1());
+  glTexCoordv(texCoords_[2]);
+  glVertexv(quad_.x2());
+  glTexCoordv(texCoords_[3]);
+  glVertexv(quad_.x3());
+  glEnd();
+}
+
+TriangleModel::TriangleModel(const Triangle& triangle)
+    : triangle_(triangle), material_(&Materials::blank()) {
+  float x = (triangle_.x1() - triangle_.x0()).length();
+  float y = (triangle_.x2() - triangle_.x0()).length();
+  texCoords_[0] = Vector(0, 0);
+  texCoords_[1] = Vector(x, 0);
+  texCoords_[2] = Vector(x, y);
+}
+
+void TriangleModel::setMaterial(const Material& m) {
+  material_ = &m;
+}
+
+void TriangleModel::setTexCoords(const Vector& t0, const Vector& t1,
+                                 const Vector& t2) {
+  texCoords_[0] = t0;
+  texCoords_[1] = t1;
+  texCoords_[2] = t2;
+}
+
+void TriangleModel::display() {
+  material_->bind();
+  displaySide(true);
+  glFrontFace(GL_CW);
+  displaySide(false);
+  glFrontFace(GL_CCW);
+}
+
+void TriangleModel::displaySide(bool front) {
+  glBegin(GL_TRIANGLES);
+  glNormalv(front ? triangle_.normal() : -triangle_.normal());
+  glTexCoordv(texCoords_[0]);
+  glVertexv(triangle_.x0());
+  glTexCoordv(texCoords_[1]);
+  glVertexv(triangle_.x1());
+  glTexCoordv(texCoords_[2]);
+  glVertexv(triangle_.x2());
+  glEnd();
+}
+
+CylinderModel::CylinderModel(const Cylinder& cylinder, const Vector& y)
+    : cylinder_(cylinder), y_(y),
+      material_(&Materials::blank()), capMaterial_(NULL),
+      quadric_(NULL) {
+  for (int i = 0; i < 15; i++) {
+    orientation_[i] = 0.f;
+  }
+  orientation_[15] = 1.f;
+}
+
+CylinderModel::~CylinderModel() {
+  if (quadric_ != NULL) {
+    gluDeleteQuadric(quadric_);
+  }
+}
+
+void CylinderModel::initialize() {
+  if (quadric_ != NULL) {
+    gluDeleteQuadric(quadric_);
+  }
+  quadric_ = gluNewQuadric();
+}
+
+float* CylinderModel::orientation() {
+  const Vector& z = cylinder_.z();
+  const Vector& y = y_;
+  const Vector& x = y.cross(z);
+  orientation_[0] = x.x; orientation_[1] = x.y; orientation_[2] = x.z;
+  orientation_[4] = y.x; orientation_[5] = y.y; orientation_[6] = y.z;
+  orientation_[8] = z.x; orientation_[9] = z.y; orientation_[10] = z.z;
+  return orientation_;
+}
+
+void CylinderModel::setMaterial(const Material& m) {
+  material_ = &m;
+}
+
+void CylinderModel::setCapMaterial(const Material& m) {
+  capMaterial_ = &m;
+}
+
+const Material& CylinderModel::capMaterial() const {
+  return (capMaterial_ == NULL) ? *material_ : *capMaterial_;
+}
+
+void CylinderModel::display() {
+  float r = cylinder_.radius();
+  float l = cylinder_.length();
+  int slices = std::max(16, (int) roundf(64 * r * r));
+  int loops = 1;
+  int stacks = 1;
+  glPushMatrix();
+  glTranslatev(cylinder_.x0());
+  glMultMatrixf(orientation());
+  gluQuadricTexture(quadric_, GL_TRUE);
+  material_->bind();
+  gluCylinder(quadric_, r, r, l, slices, stacks);
+  glPushMatrix();
+  glRotatef(180.f, 0.f, 1.f, 0.f);
+  ((capMaterial_ != NULL) ? capMaterial_ : material_)->bind();
+  gluDisk(quadric_, 0.f, r, slices, loops);
+  glPopMatrix();
+  glTranslatef(0.f, 0.f, l);
+  gluDisk(quadric_, 0.f, r, slices, loops);
+  gluQuadricTexture(quadric_, GL_FALSE);
+  glPopMatrix();
+}
+
+SphereModel::SphereModel(const Sphere& sphere)
+    : sphere_(sphere), material_(&Materials::blank()), quadric_(NULL) {
+}
+
+SphereModel::~SphereModel() {
+  if (quadric_ != NULL) {
+    gluDeleteQuadric(quadric_);
+  }
+}
+
+void SphereModel::initialize() {
+  if (quadric_ != NULL) {
+    gluDeleteQuadric(quadric_);
+  }
+  quadric_ = gluNewQuadric();
+}
+
+void SphereModel::setMaterial(const Material& m) {
+  material_ = &m;
+}
+
+void SphereModel::display() {
+  float r2 = sphere_.radius() * sphere_.radius();
+  int slices = std::max(16, (int) roundf(32 * r2));
+  int stacks = std::max(4, (int) roundf(8 * r2));
+  material_->bind();
+  glPushMatrix();
+  glTranslatev(sphere_.x());
+  gluQuadricTexture(quadric_, GL_TRUE);
+  gluSphere(quadric_, sphere_.radius(), slices, stacks);
+  gluQuadricTexture(quadric_, GL_FALSE);
+  glPopMatrix();
+}
+
+static Model* teapotModel_ = NULL;
+static Model* icosahedronModel_ = NULL;
+
+class TeapotModel : public Model {
+public:
+  virtual void display() {
+    glFrontFace(GL_CW);
+    glutSolidTeapot(1.f);
+    glFrontFace(GL_CCW);
+  }
+};
+
+class IcosahedronModel : public Model {
+public:
+  virtual void display() {
+    glutSolidIcosahedron();
+  }
+};
+
+class DisplayListModel : public Model {
+public:
+  DisplayListModel(Model* model)
+      : model_(model), displayList_(GL_NONE) {
+  }
+
+  virtual ~DisplayListModel() {
+    if (displayList_ != GL_NONE) {
+      glDeleteLists(displayList_, 1);
+    }
+    delete model_;
+  }
+
+  virtual void initialize() {
+    if (displayList_ != GL_NONE) {
+      glDeleteLists(displayList_, 1);
+    }
+    displayList_ = glGenLists(1);
+    glNewList(displayList_, GL_COMPILE);
+    model_->initialize();
+    model_->display();
+    glEndList();
+  }
+
+  virtual void display() {
+    glCallList(displayList_);
+  }
+
+private:
+  Model* model_;
+  GLuint displayList_;
+};
+
+Model* Models::compile(Model* model) {
+  return new DisplayListModel(model);
+}
+
+Model* Models::teapot() {
+  if (teapotModel_ == NULL) {
+    teapotModel_ = new TeapotModel();
+  }
+  return teapotModel_;
+}
+
+Model* Models::icosahedron() {
+  if (icosahedronModel_ == NULL) {
+    icosahedronModel_ = new IcosahedronModel();
+  }
+  return icosahedronModel_;
+}

src/model.h

@@ -0,0 +1,257 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_MODEL_H
+#define MBOSTOCK_MODEL_H
+
+#include <OpenGL/glu.h>
+
+#include "physics/shape.h"
+#include "physics/vector.h"
+
+namespace mbostock {
+
+  class Material;
+
+  /**
+   * A generic base class for OpenGL models. Provides a set of convenience
+   * methods based on the Vector class. The constructor should not assume that
+   * OpenGL is initialized; instead, use the initialize method. Note that a
+   * model may be initialized multiple times if the OpenGL context is destroyed
+   * and recreated, as is the case with resizing the window or toggling
+   * fullscreen mode.
+   */
+  class Model {
+  public:
+    virtual ~Model() {}
+
+    /** Initializes any OpenGL state used by the model. */
+    virtual void initialize() {}
+
+    /** Displays the model. */
+    virtual void display() = 0;
+
+  protected:
+
+    /** A convenience wrapper for glTranslatef. */
+    inline void glTranslatev(const Vector& v) const {
+      glTranslatef(v.x, v.y, v.z);
+    }
+
+    /** A convenience wrapper for glRotatef. */
+    inline void glRotatev(float angle, const Vector& axis) const {
+      glRotatef(angle, axis.x, axis.y, axis.z);
+    }
+
+    /** A convenience wrapper for glVertex3f. */
+    inline void glVertexv(const Vector& v) const {
+      glVertex3f(v.x, v.y, v.z);
+    }
+
+    /** A convenience wrapper for glTexCoord2f. Uses the x and y coordinate. */
+    inline void glTexCoordv(const Vector& v) const {
+      glTexCoord2f(v.x, v.y);
+    }
+
+    /** A convenience wrapper for glNormal3f. */
+    inline void glNormalv(const Vector& v) const {
+      glNormal3f(v.x, v.y, v.z);
+    }
+
+    /** A convenience wrapper for glMaterialfv. */
+    inline void glMaterialv(GLenum face, GLenum pname, const Vector& v) const {
+      GLfloat params[] = { v.x, v.y, v.z, 1.f };
+      glMaterialfv(face, pname, params);
+    }
+
+    /** A convenience wrapper for glColor3f. */
+    inline void glColorv(const Vector& v) const {
+      glColor3f(v.x, v.y, v.z);
+    }
+  };
+
+  /**
+   * A model for Wedge. A separate material may be specified for the top face of
+   * the wedge.
+   */
+  class WedgeModel : public Model {
+  public:
+    WedgeModel(const Wedge& wedge);
+
+    virtual void display();
+
+    void setMaterial(const Material& m);
+    void setTopMaterial(const Material& m);
+
+    inline const Material& material() const { return *material_; }
+    const Material& topMaterial() const;
+
+  private:
+    const Wedge& wedge_;
+    const Material* material_;
+    const Material* topMaterial_;
+  };
+
+  /**
+   * A model for AxisAlignedBox. A separate material may be specified for the
+   * top face of the box. In addition, the texture may be oriented along either
+   * the X or Z axis, and the texture may be offset (e.g., for escalatars).
+   */
+  class AxisAlignedBoxModel : public Model {
+  public:
+    AxisAlignedBoxModel(const AxisAlignedBox& box);
+
+    virtual void display();
+
+    enum Orientation { POSITIVE_X, NEGATIVE_X, POSITIVE_Z, NEGATIVE_Z };
+
+    void setMaterial(const Material& m);
+    void setTopMaterial(const Material& m);
+    void setTexOffset(float u, float v);
+    void setTexOrientation(Orientation orientation);
+
+    inline const Material& material() const { return *material_; }
+    const Material& topMaterial() const;
+
+  private:
+    const AxisAlignedBox& box_;
+    const Material* material_;
+    const Material* topMaterial_;
+    Orientation orientation_;
+    float u_, v_;
+  };
+
+  /**
+   * A model for Box. A separate material may be specified for the top face of
+   * the box.
+   */
+  class BoxModel : public Model {
+  public:
+    BoxModel(const Box& box);
+
+    virtual void display();
+
+    void setMaterial(const Material& m);
+    void setTopMaterial(const Material& m);
+
+    inline const Material& material() const { return *material_; }
+    const Material& topMaterial() const;
+
+  private:
+    const Box& box_;
+    const Material* material_;
+    const Material* topMaterial_;
+  };
+
+  /**
+   * A model for Quad. The texture coordinates may be specified explicitly, or
+   * they will be inferred from the vertex coordinates.
+   */
+  class QuadModel : public Model {
+  public:
+    QuadModel(const Quad& quad);
+
+    virtual void display();
+
+    void setMaterial(const Material& m);
+    void setTexCoords(const Vector& t0, const Vector& t1,
+                      const Vector& t2, const Vector& t3);
+
+    inline const Material& material() const { return *material_; }
+
+  private:
+    void displaySide(bool front);
+
+    const Quad& quad_;
+    const Material* material_;
+    Vector texCoords_[4];
+  };
+
+  /**
+   * A model for Triangle. The texture coordinates may be specified explicitly,
+   * or they will be inferred from the vertex coordinates.
+   */
+  class TriangleModel : public Model {
+  public:
+    TriangleModel(const Triangle& triangle);
+
+    virtual void display();
+
+    void setMaterial(const Material& m);
+    void setTexCoords(const Vector& t0, const Vector& t1, const Vector& t2);
+
+    inline const Material& material() const { return *material_; }
+
+  private:
+    void displaySide(bool front);
+
+    const Triangle& triangle_;
+    const Material* material_;
+    Vector texCoords_[3];
+  };
+
+  /**
+   * A cylinder model. A separate material may be specified for the end caps. An
+   * additional Y-axis is required to orient the texture on the cylinder's
+   * surface; this axis must be orthogonal to the main axis of the cylinder. A
+   * reference to the Y-axis vector is retained so that changes to the Y-axis
+   * will be propagated to the displayed model.
+   */
+  class CylinderModel : public Model {
+  public:
+    CylinderModel(const Cylinder& cylinder, const Vector& y);
+    virtual ~CylinderModel();
+
+    virtual void initialize();
+    virtual void display();
+
+    void setMaterial(const Material& m);
+    void setCapMaterial(const Material& m);
+
+    inline const Material& material() const { return *material_; }
+    const Material& capMaterial() const;
+
+  private:
+    float* orientation();
+
+    const Cylinder& cylinder_;
+    const Vector& y_;
+    const Material* material_;
+    const Material* capMaterial_;
+    GLUquadric* quadric_;
+    float orientation_[16];
+  };
+
+  /** A model for Sphere. */
+  class SphereModel : public Model {
+  public:
+    SphereModel(const Sphere& sphere);
+    virtual ~SphereModel();
+
+    virtual void initialize();
+    virtual void display();
+
+    void setMaterial(const Material& m);
+    inline const Material& material() const { return *material_; }
+
+  private:
+    float* orientation();
+
+    const Sphere& sphere_;
+    const Material* material_;
+    GLUquadric* quadric_;
+    float orientation_[16];
+  };
+
+  class Models {
+  public:
+    static Model* teapot();
+    static Model* icosahedron();
+    static Model* compile(Model* model);
+
+  private:
+    Models();
+  };
+
+}
+
+#endif

src/physics/constraint.cpp

@@ -0,0 +1,143 @@
+// -*- C++ -*-
+
+#include <stdio.h>
+
+#include "constraint.h"
+#include "particle.h"
+#include "shape.h"
+
+using namespace mbostock;
+
+static Projection p;
+
+bool Constraints::minDistance(Particle& a, Vector p, float d) {
+  Vector v = p - a.position;
+  float l = v.length();
+  if (l < d) {
+    a.position += v * ((l - d) / l);
+    return true;
+  }
+  return false;
+}
+
+bool Constraints::minDistance(Particle& a, Particle& b, float d) {
+  Vector v = b.position - a.position;
+  float l = v.length();
+  if (l < d) {
+    v *= (l - d) / ((a.inverseMass + b.inverseMass) * l);
+    a.position += v * a.inverseMass;
+    b.position -= v * b.inverseMass;
+    return true;
+  }
+  return false;
+}
+
+bool Constraints::maxDistance(Particle& a, Vector p, float d) {
+  Vector v = p - a.position;
+  float l = v.length();
+  if (l > d) {
+    a.position += v * ((l - d) / l);
+    return true;
+  }
+  return false;
+}
+
+bool Constraints::maxDistance(Particle& a, Particle& b, float d) {
+  Vector v = b.position - a.position;
+  float l = v.length();
+  if (l > d) {
+    v *= (l - d) / ((a.inverseMass + b.inverseMass) * l);
+    a.position += v * a.inverseMass;
+    b.position -= v * b.inverseMass;
+    return true;
+  }
+  return false;
+}
+
+bool Constraints::distance(Particle& a, Vector p, float d) {
+  Vector v = p - a.position;
+  float l = v.length();
+  a.position += v * ((l - d) / l);
+  return l > 0.f;
+}
+
+bool Constraints::distance(Particle& a, Particle& b, float d) {
+  Vector v = b.position - a.position;
+  float l = v.length();
+  v *= (l - d) / ((a.inverseMass + b.inverseMass) * l);
+  a.position += v * a.inverseMass;
+  b.position -= v * b.inverseMass;
+  return l > 0.f;
+}
+
+bool Constraints::plane(Particle& a, Vector p, Vector n) {
+  float d = n.dot(a.position - p);
+  if (d < 0) {
+    a.position -= n * d;
+    return true;
+  }
+  return false;
+}
+
+bool Constraints::plane(Particle& a, Vector p, Vector n, float kr) {
+  float d = n.dot(a.position - p);
+  if (d < 0) {
+    Vector v = a.position - a.previousPosition;
+    Vector vNormal = n * n.dot(v);
+    Vector vTangent = v - vNormal;
+    a.position -= n * d;
+    a.previousPosition = a.position - vTangent + vNormal * kr;
+    return true;
+  }
+  return false;
+}
+
+/* XXX What if p.length is 0? Use the normal? */
+
+bool Constraints::inside(Particle& a, const Shape& s, float r) {
+  p = s.project(a.position);
+  if (p.length > -r) {
+    a.position = p.x + (p.x - a.position) * (r / p.length);
+    return true;
+  }
+  return false;
+}
+
+bool Constraints::inside(Particle& a, const Shape& s, float r, float kr) {
+  p = s.project(a.position);
+  if (p.length > -r) {
+    Vector v = a.position - a.previousPosition;
+    Vector vNormal = p.normal * p.normal.dot(v);
+    Vector vTangent = v - vNormal;
+    a.position = p.x + (p.x - a.position) * (r / p.length);
+    a.previousPosition = a.position - vTangent + vNormal * kr;
+    return true;
+  }
+  return false;
+}
+
+bool Constraints::outside(Particle& a, const Shape& s, float r) {
+  p = s.project(a.position);
+  if (p.length < r) {
+    a.position = p.x - (p.x - a.position) * (r / p.length);
+    return true;
+  }
+  return false;
+}
+
+bool Constraints::outside(Particle& a, const Shape& s, float r, float kr) {
+  p = s.project(a.position);
+  if (p.length < r) {
+    Vector v = a.position - a.previousPosition;
+    Vector vNormal = p.normal * p.normal.dot(v);
+    Vector vTangent = v - vNormal;
+    a.position = p.x - (p.x - a.position) * (r / p.length);
+    a.previousPosition = a.position - vTangent + vNormal * kr;
+    return true;
+  }
+  return false;
+}
+
+const Projection& Constraints::projection() {
+  return p;
+}

src/physics/constraint.h

@@ -0,0 +1,125 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_CONSTRAINT_H
+#define MBOSTOCK_CONSTRAINT_H
+
+#include "shape.h"
+#include "vector.h"
+
+namespace mbostock {
+
+  class Particle;
+  class Shape;
+
+  class Constraints {
+  public:
+
+    /**
+     * Constrains a particle a to be distance d from position p, independent of
+     * mass.
+     */
+    static bool distance(Particle& a, Vector p, float d);
+
+    /**
+     * Constrains two particles a and b to be distance d apart. The particles
+     * are moved in opposite directions according to their relative mass, either
+     * towards or away from each other, depending on whether their actual
+     * distance is greater or less than the specified distance d.
+     */
+    static bool distance(Particle& a, Particle& b, float d);
+
+    /**
+     * Constrains a particle a to be at least distance d from position p,
+     * independent of mass.
+     */
+    static bool minDistance(Particle& a, Vector p, float d);
+
+    /**
+     * Constrains two particles a and b to be at least distance d apart. The
+     * particles are moved in opposite directions according to their relative
+     * mass, away from each other, if their actual distance is less than the
+     * specified distance d.
+     */
+    static bool minDistance(Particle& a, Particle& b, float d);
+
+    /**
+     * Constrains a particle a to be at most distance d from position p,
+     * independent of mass.
+     */
+    static bool maxDistance(Particle& a, Vector p, float d);
+
+    /**
+     * Constrains two particles a and b to be at most distance d apart. The
+     * particles are moved in opposite directions according to their relative
+     * mass, towards each other, if their actual distance is greater than the
+     * specified distance d.
+     */
+    static bool maxDistance(Particle& a, Particle& b, float d);
+
+    /**
+     * Constrains the particle so that its y-value is at least min. This
+     * implicitly uses a coefficient of restitution of zero.
+     */
+    static bool minY(Particle& a, float min) {
+      return plane(a, Vector(0, min, 0), Vector(0, 1, 0));
+    }
+
+    /**
+     * Constrains the particle so that its y-value is at least min, simulating a
+     * collision with a plane of coefficient of restitution kr. Note that the
+     * bounce can only be simulated if the particle was previously above the
+     * plane.
+     */
+    static bool minY(Particle& a, float min, float kr) {
+      return plane(a, Vector(0, min, 0), Vector(0, 1, 0), kr);
+    }
+
+    /**
+     * Constrains the particle to be above the plane defined by the point p and
+     * the normal n. This implicitly uses a coefficient of restitution of zero.
+     */
+    static bool plane(Particle& a, Vector p, Vector n);
+
+    /**
+     * Constrains the particle to be above the plane defined by the point p and
+     * the normal n, using a coefficient of restitution kr.
+     */
+    static bool plane(Particle& a, Vector p, Vector n, float kr);
+
+    /**
+     * Constrains the particle to be inside the specified shape by at least
+     * radius r. This implicitly uses a coefficient of restitution of zero.
+     */
+    static bool inside(Particle& a, const Shape& s, float r);
+
+    /**
+     * Constrains the particle to be inside the specified shape by at least
+     * radius r, using a coefficient of restitution kr.
+     */
+    static bool inside(Particle& a, const Shape& s, float r, float kr);
+
+    /**
+     * Constrains the particle to be outside the specified shape by at least
+     * radius r. This implicitly uses a coefficient of restitution of zero.
+     */
+    static bool outside(Particle& a, const Shape& s, float r);
+
+    /**
+     * Constrains the particle to be outside the specified shape by at least
+     * radius r, using a coefficient of restitution kr.
+     */
+    static bool outside(Particle& a, const Shape& s, float r, float kr);
+
+    /**
+     * Returns the projection information for the last call to a shape-based
+     * constraint.
+     */
+    static const Projection& projection();
+
+  private:
+    Constraints();
+  };
+
+}
+
+#endif

src/physics/force.cpp

@@ -0,0 +1,94 @@
+// -*- C++ -*-
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "force.h"
+#include "particle.h"
+#include "vector.h"
+
+using namespace mbostock;
+
+static const float epsilon = 1E-20;
+
+void UnaryForce::apply(Particle& p) {
+  p.force += force(p);
+}
+
+void BinaryForce::apply(Particle& a, Particle& b) {
+  Vector f = force(a, b);
+  a.force += f;
+  b.force -= f;
+}
+
+RandomForce::RandomForce(float k)
+    : k_(k) {
+}
+
+Vector RandomForce::force(Particle& p) {
+  return Vector::randomVector(k_);
+}
+
+GravitationalForce::GravitationalForce(float g)
+    : g_(g) {
+}
+
+Vector GravitationalForce::force(const Particle& p) {
+  return Vector(0, -g_ / p.inverseMass, 0);
+}
+
+LinearDragForce::LinearDragForce(float kd)
+    : kd_(kd) {
+}
+
+Vector LinearDragForce::force(const Particle& p) {
+  return p.velocity() * -kd_;
+}
+
+QuadraticDragForce::QuadraticDragForce(float kd)
+    : kd_(kd) {
+}
+
+Vector QuadraticDragForce::force(const Particle& p) {
+  const Vector& v = p.velocity();
+  return v * v.length() * -kd_;
+}
+
+HookeForce::HookeForce(float r, float ks)
+    : r_(r), ks_(ks) {
+}
+
+Vector HookeForce::force(const Particle& a, const Particle& b) {
+  Vector l = a.position - b.position;
+  float ll = l.length();
+  if (ll == 0.f) {
+    l = Vector::randomVector(epsilon);
+    ll = epsilon;
+  }
+  return l * -ks_ * (ll - r_) / ll;
+}
+
+DampenedHookeForce::DampenedHookeForce(float r, float ks, float kd)
+    : r_(r), ks_(ks), kd_(kd) {
+}
+
+Vector DampenedHookeForce::force(const Particle& a, const Particle& b) {
+  Vector l = a.position - b.position;
+  Vector dl = a.velocity() - b.velocity();
+  float ll = l.length();
+  if (ll == 0.f) {
+    l = Vector::randomVector(epsilon);
+    ll = epsilon;
+  }
+  return l * -(ks_ * (ll - r_) + kd_ * dl.dot(l) / ll) / ll;
+}
+
+CoulombForce::CoulombForce(float e)
+    : e_(e) {
+}
+
+Vector CoulombForce::force(const Particle& a, const Particle& b) {
+  Vector l = a.position - b.position;
+  float ll = l.length();
+  return -l * e_ / (ll * ll * ll);
+}

src/physics/force.h

@@ -0,0 +1,105 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_FORCE_H
+#define MBOSTOCK_FORCE_H
+
+#include "vector.h"
+
+namespace mbostock {
+
+  class Particle;
+
+  class UnaryForce {
+  public:
+    virtual ~UnaryForce() {}
+
+    void apply(Particle& p);
+
+    virtual Vector force(const Particle& p)  = 0;
+  };
+
+  class BinaryForce {
+  public:
+    virtual ~BinaryForce() {}
+
+    void apply(Particle& a, Particle& b);
+
+    virtual Vector force(const Particle& a, const Particle& b) = 0;
+  };
+
+  class RandomForce : public UnaryForce {
+  public:
+    RandomForce(float k);
+
+    virtual Vector force(Particle& p);
+
+  private:
+    float k_;
+  };
+
+  class GravitationalForce : public UnaryForce {
+  public:
+    GravitationalForce(float g);
+
+    virtual Vector force(const Particle& p);
+
+  private:
+    float g_;
+  };
+
+  class LinearDragForce : public UnaryForce {
+  public:
+    LinearDragForce(float kd);
+
+    virtual Vector force(const Particle& p);
+
+  private:
+    float kd_;
+  };
+
+  class QuadraticDragForce : public UnaryForce {
+  public:
+    QuadraticDragForce(float kd);
+
+    virtual Vector force(const Particle& p);
+
+  private:
+    float kd_;
+  };
+
+  class HookeForce : public BinaryForce {
+  public:
+    HookeForce(float r, float ks);
+
+    virtual Vector force(const Particle& a, const Particle& b);
+
+  private:
+    float r_;
+    float ks_;
+  };
+
+  class DampenedHookeForce : public BinaryForce {
+  public:
+    DampenedHookeForce(float r, float ks, float kd);
+
+    virtual Vector force(const Particle& a, const Particle& b);
+
+  private:
+    float r_;
+    float ks_;
+    float kd_;
+  };
+
+  class CoulombForce : public BinaryForce {
+  public:
+    CoulombForce(float e);
+
+    virtual Vector force(const Particle& a, const Particle& b);
+
+  private:
+    float e_;
+  };
+
+}
+
+#endif

src/physics/particle.cpp

@@ -0,0 +1,37 @@
+// -*- C++ -*-
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "particle.h"
+
+using namespace mbostock;
+
+Particle::Particle() : inverseMass(1.f) {
+}
+
+Vector Particle::velocity() const {
+  /* Note: ignores force * (inverseMass * timeStep / 2). */
+  return (position - previousPosition) / ParticleSimulator::timeStep();
+}
+
+ParticleSimulator::ParticleSimulator()
+    : drag_(1.f) {
+}
+
+ParticleSimulator::ParticleSimulator(float drag) 
+    : drag_(drag) {
+}
+
+float ParticleSimulator::timeStep() {
+  return .003f; // TODO .005 and interpolate
+}
+
+void ParticleSimulator::step(Particle& p) const {
+  static const float timeStepSquared = timeStep() * timeStep();
+
+  Vector p0 = p.previousPosition;
+  p.previousPosition = p.position;
+  p.position += (p.position - p0) * drag_
+      + p.force * (p.inverseMass * timeStepSquared);
+}

src/physics/particle.h

@@ -0,0 +1,39 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_PARTICLE_H
+#define MBOSTOCK_PARTICLE_H
+
+#include "vector.h"
+
+namespace mbostock {
+
+  class Particle {
+  public:
+    Particle();
+
+    float inverseMass;
+    Vector position;
+    Vector previousPosition;
+    Vector force;
+
+    Vector velocity() const;
+  };
+
+  class ParticleSimulator {
+  public:
+    ParticleSimulator();
+    ParticleSimulator(float drag);
+
+    void step(Particle& p) const;
+
+    static float timeStep();
+
+  private:
+    static float timeStep_;
+    static float timeStepSquared_;
+    float drag_;
+  };
+
+}
+
+#endif

src/physics/rotation.cpp

@@ -0,0 +1,89 @@
+// -*- C++ -*-
+
+#include <math.h>
+
+#include "particle.h"
+#include "rotation.h"
+
+using namespace mbostock;
+
+Rotation::Rotation(const Vector& origin, const Vector& axis,
+                   float speed, float angle)
+    : origin_(origin), axis_(axis), startAngle_(angle),
+      speed_(speed * ParticleSimulator::timeStep()), angle_(angle) {
+  update();
+}
+
+void Rotation::reset() {
+  angle_ = startAngle_;
+}
+
+void Rotation::step() {
+  if (!enabled()) {
+    return;
+  }
+  angle_ += speed_;
+  angle_ = fmodf(angle_, 360.f);
+  update();
+}
+
+void Rotation::update() {
+  float r = angle_ * (2.f * M_PI / 360.f);
+  float c = cosf(r);
+  float s = sinf(r);
+
+  /* This look familiar to anyone? /wink */
+  matrix_[0] = axis_.x * axis_.x * (1.f - c) + c;
+  matrix_[1] = axis_.x * axis_.y * (1.f - c) - axis_.z * s;
+  matrix_[2] = axis_.x * axis_.z * (1.f - c) + axis_.y * s;
+  matrix_[3] = axis_.y * axis_.x * (1.f - c) + axis_.z * s;
+  matrix_[4] = axis_.y * axis_.y * (1.f - c) + c;
+  matrix_[5] = axis_.y * axis_.z * (1.f - c) - axis_.x * s;
+  matrix_[6] = axis_.x * axis_.z * (1.f - c) - axis_.y * s;
+  matrix_[7] = axis_.y * axis_.z * (1.f - c) + axis_.x * s;
+  matrix_[8] = axis_.z * axis_.z * (1.f - c) + c;
+}
+
+Vector Rotation::point(const Vector& x) const {
+  return origin_ + vector(x - origin_);
+}
+
+Vector Rotation::pointInverse(const Vector& x) const {
+  return origin_ + vectorInverse(x - origin_);
+}
+
+Vector Rotation::vector(const Vector& x) const {
+  return Vector(
+      matrix_[0] * x.x + matrix_[1] * x.y + matrix_[2] * x.z,
+      matrix_[3] * x.x + matrix_[4] * x.y + matrix_[5] * x.z,
+      matrix_[6] * x.x + matrix_[7] * x.y + matrix_[8] * x.z);
+}
+
+Vector Rotation::vectorInverse(const Vector& x) const {
+  return Vector(
+      matrix_[0] * x.x + matrix_[3] * x.y + matrix_[6] * x.z,
+      matrix_[1] * x.x + matrix_[4] * x.y + matrix_[7] * x.z,
+      matrix_[2] * x.x + matrix_[5] * x.y + matrix_[8] * x.z);
+}
+
+Vector Rotation::velocity(const Vector& x) const {
+  return enabled()
+      ? (origin_ - x).cross(axis_) * speed_ * (2.f * M_PI / 360.f)
+      : Vector::ZERO();
+}
+
+RotatingShape::RotatingShape(const Shape& s, const Rotation& r)
+    : shape_(s), rotation_(r) {
+}
+
+bool RotatingShape::intersects(const Sphere& s) const {
+  Sphere rs(rotation_.pointInverse(s.x()), s.radius());
+  return shape_.intersects(rs);
+}
+
+Projection RotatingShape::project(const Vector& x) const {
+  Projection p = shape_.project(rotation_.pointInverse(x));
+  p.x = rotation_.point(p.x);
+  p.normal = rotation_.vector(p.normal);
+  return p;
+}

src/physics/rotation.h

@@ -0,0 +1,79 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_ROTATION_H
+#define MBOSTOCK_ROTATION_H
+
+#include "shape.h"
+#include "vector.h"
+#include "transform.h"
+
+namespace mbostock {
+
+  /**
+   * Computes a rotation matrix for a given origin, axis and angle. This is
+   * roughly the equivalent to a glTranslatef and glRotatef call in OpenGL.
+   */
+  class Rotation : public Transform {
+  public:
+    Rotation(const Vector& origin, const Vector& axis,
+             float speed, float angle);
+
+    /** Rotates the specified point. */
+    Vector point(const Vector& x) const;
+
+    /** Inverse-rotates the specified point. */
+    Vector pointInverse(const Vector& x) const;
+
+    /** Rotates the specified vector. */
+    Vector vector(const Vector& x) const;
+
+    /** Inverse-rotates the specified vector. */
+    Vector vectorInverse(const Vector& x) const;
+
+    /** Returns the rotation's origin. */
+    inline const Vector& origin() const { return origin_; }
+
+    /** Returns the rotation's axis. */
+    inline const Vector& axis() const { return axis_; }
+
+    /** Returns the rotation's angle (in degrees, in the range [0, 360]). */
+    inline float angle() const { return angle_; }
+
+    /** Returns the speed of the rotation (in degrees per time step). */
+    inline float speed() const { return speed_; }
+
+    /** Advances the rotation by one time step. */
+    virtual void step();
+
+    /** Resets the rotation. */
+    virtual void reset();
+
+    /** Returns the velocity of the given point. */
+    Vector velocity(const Vector& x) const;
+
+  private:
+    void update();
+
+    Vector origin_;
+    Vector axis_;
+    float startAngle_;
+    float angle_;
+    float speed_;
+    float matrix_[9];
+  };
+
+  class RotatingShape : public Shape {
+  public:
+    RotatingShape(const Shape& s, const Rotation& r);
+
+    virtual bool intersects(const Sphere& s) const;
+    virtual Projection project(const Vector& x) const;
+
+  private:
+    const Shape& shape_;
+    const Rotation& rotation_;
+  };
+
+}
+
+#endif

src/physics/shape.cpp

@@ -0,0 +1,514 @@
+// -*- C++ -*-
+
+#include <algorithm>
+#include <math.h>
+
+#include "shape.h"
+
+using namespace mbostock;
+
+Projection::Projection() : length(0) {
+}
+
+Projection::Projection(const Vector& x, float length)
+    : x(x), length(length) {
+}
+
+Projection::Projection(const Vector& x, float length, const Vector& normal)
+    : x(x), length(length), normal(normal) {
+}
+
+bool Projection::operator<(const Projection& p) const {
+  return fabsf(length) < fabsf(p.length);
+}
+
+bool Projection::operator<=(const Projection& p) const {
+  return fabsf(length) <= fabsf(p.length);
+}
+
+bool Projection::operator>(const Projection& p) const {
+  return fabsf(length) > fabsf(p.length);
+}
+
+bool Projection::operator>=(const Projection& p) const {
+  return fabsf(length) >= fabsf(p.length);
+}
+
+Sphere::Sphere()
+    : r_(0.f), r2_(0.f) {
+}
+
+Sphere::Sphere(const Vector& x, float radius)
+    : x_(x), r_(radius), r2_(radius * radius) {
+}
+
+bool Sphere::above(const Plane& p) const {
+  /* Derived from Point-Plane Distance on MathWorld. */
+  return p.n_.dot(x_ - p.x_) > r_;
+}
+
+bool Sphere::below(const Plane& p) const {
+  /* Derived from Point-Plane Distance on MathWorld. */
+  return p.n_.dot(x_ - p.x_) < -r_;
+}
+
+bool Sphere::intersects(const Sphere& s) const {
+  /* Derived from Moller-Haines section 16.13.1. */
+  float r = r_ + s.r_;
+  return (s.x_ - x_).squared() < (r * r);
+}
+
+Projection Sphere::project(const Vector& x) const {
+  Vector v = x - x_;
+  float l = v.length();
+  float d = l - r_;
+  Vector n = v / l;
+  return Projection(x - n * d, fabsf(d), (d < 0.f) ? -n : n);
+}
+
+LineSegment::LineSegment() {
+}
+
+LineSegment::LineSegment(const Vector& x0, const Vector& x1)
+    : x0_(x0), x1_(x1), x01_(x1 - x0), l2_(x01_.dot(x01_)) {
+}
+
+bool LineSegment::intersects(const Sphere& s) const {
+  /* Derived from Point-Line Distance--3-Dimensional on MathWorld. */
+  return (x01_.cross(x0_ - s.x_).squared() / l2_) < s.r2_;
+}
+
+Projection LineSegment::project(const Vector& p) const {
+  /* Derived from Point-Line Distance--3-Dimensional on MathWorld. */
+  float t = (p - x0_).dot(x01_) / l2_;
+  Vector x = (t <= 0.f) ? x0_ : ((t >= 1.f) ? x1_ : (x0_ + x01_ * t));
+  return Projection(x, (p - x).length());
+}
+
+Plane::Plane() {
+}
+
+Plane::Plane(const Vector& x, const Vector& normal)
+    : x_(x), n_(normal) {
+}
+
+bool Plane::intersects(const Sphere& s) const {
+  /* Derived from Point-Plane Distance on MathWorld. */
+  return fabsf(n_.dot(s.x_ - x_)) < s.r_;
+}
+
+Projection Plane::project(const Vector& p) const {
+  /* Derived from Point-Plane Distance on MathWorld. */
+  float l = n_.dot(p - x_);
+  return Projection(p - n_ * l, fabsf(l), (l < 0.f) ? -n_ : n_);
+}
+
+Triangle::Triangle() {
+}
+
+Triangle::Triangle(const Vector& x0, const Vector& x1, const Vector& x2)
+    : x01_(x0, x1), x12_(x1, x2), x20_(x2, x0),
+      p_(x0, (x1 - x0).cross(x2 - x1).normalize()) {
+}
+
+bool Triangle::intersects(const Sphere& s) const {
+  /* Derived from ERIT section 2.3 (reordered). */
+  Projection p0 = p_.project(s.x_);
+  if (p0.length > s.r_) {
+    return false;
+  }
+  if (x01_.intersects(s) || x12_.intersects(s) || x20_.intersects(s)) {
+    return true;
+  }
+  return contains(p0.x);
+}
+
+Projection Triangle::project(const Vector& p) const {
+  Projection pp = p_.project(p);
+  if (contains(pp.x)) {
+    return pp;
+  }
+  Projection p01 = x01_.project(pp.x);
+  Projection p12 = x12_.project(pp.x);
+  Projection p20 = x20_.project(pp.x);
+  Projection pb = std::min(std::min(p01, p12), p20);
+  pb.length = sqrtf(pb.length * pb.length + pp.length * pp.length);
+  pb.normal = pp.normal;
+  return pb;
+}
+
+bool Triangle::contains(const Vector& p) const {
+  return ((x0() - x1()).cross(p - x0()).dot(p_.normal()) < 0.f)
+      && ((x1() - x2()).cross(p - x1()).dot(p_.normal()) < 0.f)
+      && ((x2() - x0()).cross(p - x2()).dot(p_.normal()) < 0.f);
+}
+
+Quad::Quad() {
+}
+
+Quad::Quad(const Vector& x0, const Vector& x1,
+           const Vector& x2, const Vector& x3)
+    : x01_(x0, x1), x12_(x1, x2), x23_(x2, x3), x30_(x3, x0),
+      p_(x0, (x1 - x0).cross(x2 - x1).normalize()) {
+}
+
+bool Quad::intersects(const Sphere& s) const {
+  /* Derived from Triangle-Sphere test above. */
+  Projection p0 = p_.project(s.x());
+  if (p0.length > s.radius()) {
+    return false;
+  }
+  if (x01_.intersects(s)
+      || x12_.intersects(s)
+      || x23_.intersects(s)
+      || x30_.intersects(s)) {
+    return true;
+  }
+  return contains(p0.x);
+}
+
+Projection Quad::project(const Vector& p) const {
+  Projection pp = p_.project(p);
+  if (contains(pp.x)) {
+    return pp;
+  }
+  Projection p01 = x01_.project(pp.x);
+  Projection p12 = x12_.project(pp.x);
+  Projection p23 = x23_.project(pp.x);
+  Projection p30 = x30_.project(pp.x);
+  Projection pb = std::min(std::min(std::min(p01, p12), p23), p30);
+  pb.length = sqrtf(pb.length * pb.length + pp.length * pp.length);
+  pb.normal = pp.normal;
+  return pb;
+}
+
+bool Quad::contains(const Vector& p) const {
+  return ((x0() - x1()).cross(p - x0()).dot(p_.normal()) < 0.f)
+      && ((x1() - x2()).cross(p - x1()).dot(p_.normal()) < 0.f)
+      && ((x2() - x3()).cross(p - x2()).dot(p_.normal()) < 0.f)
+      && ((x3() - x0()).cross(p - x3()).dot(p_.normal()) < 0.f);
+}
+
+Wedge::Wedge() {
+}
+
+Wedge::Wedge(const Vector& x0, const Vector& x1,
+             const Vector& x2, const Vector& x3)
+    : top_(x0, x1, x2, x3),
+      right_(x2, x1, Vector(x1.x, x0.y, x1.z), Vector(x2.x, x0.y, x2.z)),
+      bottom_(x0, x3, right_.x3(), right_.x2()),
+      front_(x1, x0, right_.x2()),
+      back_(x3, x2, right_.x3()) {
+}
+
+bool Wedge::intersects(const Sphere& s) const {
+  /* Derived (very approximately) from Moller-Haines section 16.14.2. */
+  if (s.above(top_.p_)
+      || s.above(right_.p_)
+      || s.above(front_.p_)
+      || s.above(back_.p_)
+      || s.above(bottom_.p_)) {
+    return false;
+  }
+  if (s.below(top_.p_)
+      && s.below(right_.p_)
+      && s.below(front_.p_)
+      && s.below(back_.p_)
+      && s.below(bottom_.p_)) {
+    return true;
+  }
+  return top_.intersects(s)
+    || right_.intersects(s)
+    || front_.intersects(s)
+    || back_.intersects(s)
+    || bottom_.intersects(s);
+}
+
+Projection Wedge::project(const Vector& p) const {
+  Projection pt = top_.project(p);
+  Projection pr = right_.project(p);
+  Projection pf = front_.project(p);
+  Projection pb = back_.project(p);
+  Projection pd = bottom_.project(p);
+  return std::min(std::min(std::min(std::min(pt, pr), pf), pb), pd);
+}
+
+AxisAlignedBox::AxisAlignedBox() {
+}
+
+AxisAlignedBox::AxisAlignedBox(const Vector& min, const Vector& max)
+    : min_(min), max_(max) {
+}
+
+bool AxisAlignedBox::intersects(const Sphere& s) const {
+  /* Derived from Moller-Haines section 16.13.2. */
+  Vector e = Vector::max(min_ - s.x_, Vector::ZERO())
+           + Vector::max(s.x_ - max_, Vector::ZERO());
+  return e.squared() < s.r2_;
+}
+
+bool AxisAlignedBox::contains(const Vector& p) const {
+  return ((p.x >= min_.x) && (p.x < max_.x)
+          && (p.y >= min_.y) && (p.y < max_.y)
+          && (p.z >= min_.z) && (p.z < max_.z));
+}
+
+Projection AxisAlignedBox::project(const Vector& p) const {
+  return contains(p) ? projectOut(p) : projectIn(p);
+}
+
+Projection AxisAlignedBox::projectOut(const Vector& p) const {
+  Vector x = p;
+  Vector n;
+  float l;
+
+  float dx1 = p.x - min_.x;
+  float dx2 = max_.x - p.x;
+  float dx = std::min(dx1, dx2);
+
+  float dy1 = p.y - min_.y;
+  float dy2 = max_.y - p.y;
+  float dy = std::min(dy1, dy2);
+
+  float dz1 = p.z - min_.z;
+  float dz2 = max_.z - p.z;
+  float dz = std::min(dz1, dz2);
+
+  if ((dx <= dy) && (dx <= dz)) {
+    if (dx1 <= dx2) {
+      l = dx1;
+      x.x = min_.x;
+      n.x = 1.f;
+    } else {
+      l = dx2;
+      x.x = max_.x;
+      n.x = -1.f;
+    }
+  } else if (dy <= dz) {
+    if (dy1 <= dy2) {
+      l = dy1;
+      x.y = min_.y;
+      n.y = 1.f;
+    } else {
+      l = dy2;
+      x.y = max_.y;
+      n.y = -1.f;
+    }
+  } else {
+    if (dz1 <= dz2) {
+      l = dz1;
+      x.z = min_.z;
+      n.z = 1.f;
+    } else {
+      l = dz2;
+      x.z = max_.z;
+      n.z = -1.f;
+    }
+  }
+
+  return Projection(x, -l, n);
+}
+
+Projection AxisAlignedBox::projectIn(const Vector& p) const {
+  Vector x = p;
+  Vector n;
+
+  if (p.x < min_.x) {
+    x.x = min_.x;
+    n = -Vector::X();
+  } else if (p.x >= max_.x) {
+    x.x = max_.x;
+    n = Vector::X();
+  }
+
+  if (p.z < min_.z) {
+    x.z = min_.z;
+    n = -Vector::Z();
+  } else if (p.z >= max_.z) {
+    x.z = max_.z;
+    n = Vector::Z();
+  }
+
+  if (p.y < min_.y) {
+    x.y = min_.y;
+    n = -Vector::Y();
+  } else if (p.y >= max_.y) {
+    x.y = max_.y;
+    n = Vector::Y();
+  }
+
+  /*
+   * TODO: If the particle does not project directly onto the triangle (and
+   * similarly for quads), we should interpolate the normal of the projection so
+   * that particles appear to bounce off the corner.
+   *
+   * On the other hand, this leads to some weird behavior when applying
+   * directional friction when Polly is on the edge of a platform. So for now,
+   * just use axis normals.
+   */
+
+  return Projection(x, (p - x).length(), n);
+}
+
+Box::Box() {
+}
+
+Box::Box(const AxisAlignedBox& box)
+    : top_(box.x4(), box.x5(), box.x6(), box.x7()),
+      bottom_(box.x0(), box.x1(), box.x2(), box.x3()),
+      left_(box.x0(), box.x3(), box.x6(), box.x5()),
+      right_(box.x1(), box.x4(), box.x7(), box.x2()),
+      front_(box.x2(), box.x7(), box.x6(), box.x3()),
+      back_(box.x0(), box.x5(), box.x4(), box.x1()) {
+}
+
+Box::Box(const Vector& c, const Vector& x, const Vector& y, const Vector& z) {
+  Vector x0 = c - x - y - z;
+  Vector x1 = c + x - y - z;
+  Vector x2 = c + x - y + z;
+  Vector x3 = c - x - y + z;
+  Vector x4 = c + x + y - z;
+  Vector x5 = c - x + y - z;
+  Vector x6 = c - x + y + z;
+  Vector x7 = c + x + y + z;
+  top_ = Quad(x4, x5, x6, x7);
+  bottom_ = Quad(x0, x1, x2, x3);
+  left_ = Quad(x0, x3, x6, x5);
+  right_ = Quad(x1, x4, x7, x2);
+  front_ = Quad(x2, x7, x6, x3);
+  back_ = Quad(x0, x5, x4, x1);
+}
+
+Box::Box(const Vector& x0, const Vector& x1, const Vector& x2, const Vector& x3,
+         const Vector& x4, const Vector& x5, const Vector& x6, const Vector& x7)
+    : top_(x4, x5, x6, x7),
+      bottom_(x0, x1, x2, x3),
+      left_(x0, x3, x6, x5),
+      right_(x1, x4, x7, x2),
+      front_(x2, x7, x6, x3),
+      back_(x0, x5, x4, x1) {
+}
+
+bool Box::intersects(const Sphere& s) const {
+  /* Derived (very approximately) from Moller-Haines section 16.14.2. */
+  if (s.above(top_.plane())
+      || s.above(left_.plane())
+      || s.above(right_.plane())
+      || s.above(front_.plane())
+      || s.above(back_.plane())
+      || s.above(bottom_.plane())) {
+    return false;
+  }
+  if (s.below(top_.plane())
+      && s.below(left_.plane())
+      && s.below(right_.plane())
+      && s.below(front_.plane())
+      && s.below(back_.plane())
+      && s.below(bottom_.plane())) {
+    return true;
+  }
+  return top_.intersects(s)
+    || left_.intersects(s)
+    || right_.intersects(s)
+    || front_.intersects(s)
+    || back_.intersects(s)
+    || bottom_.intersects(s);
+}
+
+Projection Box::project(const Vector& v) const {
+  Projection p = top_.project(v);
+  p = std::min(p, bottom_.project(v));
+  p = std::min(p, left_.project(v));
+  p = std::min(p, right_.project(v));
+  p = std::min(p, front_.project(v));
+  p = std::min(p, back_.project(v));
+  return p;
+}
+
+Cylinder::Cylinder() {
+}
+
+Cylinder::Cylinder(const Vector& x0, const Vector& x1, float radius)
+    : axis_(x0, x1), l_(sqrtf(axis_.l2_)), r_(radius), v_(axis_.x01_ / l_) {
+}
+
+bool Cylinder::intersects(const Sphere& s) const {
+  /* Derived from ERIT section 2.5. */
+  Vector pa = s.x_ - axis_.x0_;
+  float pqdotpa = axis_.x01_.dot(pa);
+  float rs1;
+
+  /* P -> Q -> B */
+  if (pqdotpa > axis_.l2_) {
+    float f = pqdotpa - axis_.l2_;
+    float d2 = (f * f) / axis_.l2_;
+    if (d2 > s.r2_) {
+      return false;
+    }
+    rs1 = sqrtf(s.r2_ - d2);
+  }
+
+  /* B -> P -> Q */
+  else if (pqdotpa < 0.f) {
+    float qpdotqa = -axis_.x01_.dot(s.x_ - axis_.x1_);
+    float f = qpdotqa - axis_.l2_;
+    float d2 = (f * f) / axis_.l2_;
+    if (d2 > s.r2_) {
+      return false;
+    }
+    rs1 = sqrtf(s.r2_ - d2);
+  }
+
+  /* P -> B -> Q */
+  else {
+    rs1 = s.r_;
+  }
+
+  float ld2 = axis_.x01_.cross(pa).squared() / axis_.l2_;
+  float r = r_ + rs1;
+  return ld2 <= (r * r);
+}
+
+Projection Cylinder::project(const Vector& p) const {
+  /* Derived from ERIT section 2.5. */
+  Vector pa = p - axis_.x0_;
+  float pqdotpa = axis_.x01_.dot(pa);
+
+  /* P -> Q -> B */
+  if (pqdotpa > axis_.l2_) {
+    Vector x = Plane(axis_.x1_, v_).project(p).x;
+
+    Vector v = x - axis_.x1_;
+    float l = v.length();
+    if (l > r_) {
+      float d = l - r_;
+      Vector n = v / l;
+      x -= n * d;
+    }
+
+    return Projection(x, (p - x).length(), v_);
+  }
+
+  /* B -> P -> Q */
+  else if (pqdotpa < 0.f) {
+    Vector x = Plane(axis_.x0_, v_).project(p).x;
+
+    Vector v = x - axis_.x0_;
+    float l = v.length();
+    if (l > r_) {
+      float d = l - r_;
+      Vector n = v / l;
+      x -= n * d;
+    }
+
+    return Projection(x, (p - x).length(), -v_);
+  }
+
+  /* P -> B -> Q */
+  Vector b = axis_.x0_ + axis_.x01_ * pqdotpa / axis_.l2_;
+  Vector v = p - b;
+  float l = v.length();
+  float d = l - r_;
+  Vector n = v / l;
+  return Projection(p - n * d, fabsf(d), (d < 0.f) ? -n : n);
+}

src/physics/shape.h

@@ -0,0 +1,431 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_SHAPE_H
+#define MBOSTOCK_SHAPE_H
+
+#include "vector.h"
+
+namespace mbostock {
+
+  class Plane;
+  class Sphere;
+
+  /**
+   * Represents the projection of a point onto the closest corresponding point
+   * on the surface of a shape. A projection is used when the bounding sphere of
+   * a particle interpenetrates a shape and, as this constitutes a collision,
+   * the sphere center is projected onto the surface of the shape to determine
+   * how to move the particle out in response to the collision.
+   *
+   * When a point is projected onto a surface, the projection includes the
+   * surface normal. However, if the point is projected onto a line, the
+   * corresponding normal is undefined.
+   *
+   * Projections define a natural order based on the length of the projection,
+   * so that it is convenient to find the shortest projection given a number of
+   * candidate projections.
+   */
+  class Projection {
+  public:
+    Projection();
+    Projection(const Vector& x, float length);
+    Projection(const Vector& x, float length, const Vector& normal);
+
+    /** The location of the projected point on the shape's surface. */
+    Vector x;
+
+    /** The distance from the original point to the projected point. */
+    float length;
+
+    /** The surface normal at the projected point. */
+    Vector normal;
+
+    bool operator<(const Projection& p) const;
+    bool operator<=(const Projection& p) const;
+    bool operator>(const Projection& p) const;
+    bool operator>=(const Projection& p) const;
+  };
+
+  /**
+   * Represents a shape in three dimensions. Since particles are represented
+   * solely as spheres, shapes need only define intersection tests with spheres,
+   * and how to project points onto the closest point on the surface of the
+   * shape.
+   */
+  class Shape {
+  public:
+    virtual ~Shape() {}
+
+    /** Returns true if the specified sphere intersects with this shape. */
+    virtual bool intersects(const Sphere& s) const = 0;
+
+    /**
+     * Projects the specified point onto the surface of the shape. The
+     * projection specifies the closest point on the surface of the shape to the
+     * specified point, the pre-computed length of the projection (the distance
+     * between two points), and the surface normal at the projected point, if
+     * defined.
+     */
+    virtual Projection project(const Vector& x) const = 0;
+  };
+
+  /** Represents a sphere using a center point and a radius. */
+  class Sphere : public Shape {
+  public:
+    Sphere();
+    Sphere(const Vector& x, float radius);
+
+    /** The center point of the sphere. */
+    inline Vector& x() { return x_; }
+    inline const Vector& x() const { return x_; }
+
+    /** The radius of the sphere. */
+    inline float radius() const { return r_; }
+
+    /** Returns true if the sphere is entirely above the specified plane. */
+    bool above(const Plane& p) const;
+
+    /** Returns true if the sphere is entirely below the specified plane. */
+    bool below(const Plane& p) const;
+
+    virtual bool intersects(const Sphere& s) const;
+    virtual Projection project(const Vector& p) const;
+
+  private:
+    Vector x_;
+    float r_;
+    float r2_;
+
+    friend class AxisAlignedBox;
+    friend class Cylinder;
+    friend class LineSegment;
+    friend class Plane;
+    friend class Triangle;
+  };
+
+  /** Represents a line segment using a start point and an end point. */
+  class LineSegment : public Shape {
+  public:
+    LineSegment();
+    LineSegment(const Vector& x0, const Vector& x1);
+
+    /** The start point of the line segment. */
+    inline const Vector& x0() const { return x0_; }
+
+    /** The end point of the line segment. */
+    inline const Vector& x1() const { return x1_; }
+
+    virtual bool intersects(const Sphere& s) const;
+    virtual Projection project(const Vector& p) const;
+
+  private:
+    Vector x0_;
+    Vector x1_;
+    Vector x01_;
+    float l2_;
+
+    friend class Cylinder;
+  };
+
+  /** Represents a plane using a point on the plane and a normal. */
+  class Plane : public Shape {
+  public:
+    Plane();
+    Plane(const Vector& x, const Vector& normal);
+
+    /** A point on the plane. */
+    inline const Vector& x() const { return x_; }
+
+    /** The normal. */
+    inline const Vector& normal() const { return n_; }
+
+    virtual bool intersects(const Sphere& s) const;
+    virtual Projection project(const Vector& p) const;
+
+  private:
+    Vector x_;
+    Vector n_;
+
+    friend class Sphere;
+  };
+
+  /** Represents a triangle using three coplanar points. */
+  class Triangle : public Shape {
+  public:
+    Triangle();
+    Triangle(const Vector& x0, const Vector& x1, const Vector& x2);
+
+    /** The first point of the triangle. */
+    inline const Vector& x0() const { return x01_.x0(); }
+
+    /** The second point of the triangle. */
+    inline const Vector& x1() const { return x12_.x0(); }
+
+    /** The third point of the triangle. */
+    inline const Vector& x2() const { return x20_.x0(); }
+
+    /** The normal of the triangle, using counterclockwise orientation. */
+    inline const Vector& normal() const { return p_.normal(); }
+
+    virtual bool intersects(const Sphere& s) const;
+    virtual Projection project(const Vector& p) const;
+
+  private:
+    bool contains(const Vector& p) const;
+
+    LineSegment x01_;
+    LineSegment x12_;
+    LineSegment x20_;
+    Plane p_;
+
+    friend class Wedge;
+  };
+
+  /** Represents a quad using four coplanar points. */
+  class Quad : public Shape {
+  public:
+    Quad();
+    Quad(const Vector& x0, const Vector& x1,
+         const Vector& x2, const Vector& x3);
+
+    /** The first point of the quad. */
+    inline const Vector& x0() const { return x01_.x0(); }
+
+    /** The second point of the quad. */
+    inline const Vector& x1() const { return x12_.x0(); }
+
+    /** The third point of the quad. */
+    inline const Vector& x2() const { return x23_.x0(); }
+
+    /** The fourth point of the quad. */
+    inline const Vector& x3() const { return x30_.x0(); }
+
+    /** The normal of the quad, using counterclockwise orientation. */
+    inline const Vector& normal() const { return p_.normal(); }
+
+    /** The plane of the quad. */
+    inline const Plane& plane() const { return p_; }
+
+    virtual bool intersects(const Sphere& s) const;
+    virtual Projection project(const Vector& p) const;
+
+  private:
+    bool contains(const Vector& p) const;
+
+    LineSegment x01_;
+    LineSegment x12_;
+    LineSegment x23_;
+    LineSegment x30_;
+    Plane p_;
+
+    friend class Wedge;
+  };
+
+  /**
+   * Represents a wedge using the four coplanar points of the top surface.
+   * Points x1 and x2 define the upper edge of the wedge; points x4 and x5 lie
+   * directly beneath them to define the fifth side.
+   */
+  class Wedge : public Shape {
+  public:
+    Wedge();
+    Wedge(const Vector& x0, const Vector& x1,
+          const Vector& x2, const Vector& x3);
+
+    /** Returns the first (lower) point of the wedge. */
+    inline const Vector& x0() const { return top_.x0(); }
+
+    /** Returns the second (upper) point of the wedge. */
+    inline const Vector& x1() const { return top_.x1(); }
+
+    /** Returns the third (upper) point of the wedge. */
+    inline const Vector& x2() const { return top_.x2(); }
+
+    /** Returns the forth (upper) point of the wedge. */
+    inline const Vector& x3() const { return top_.x3(); }
+
+    /** Returns the fifth point of the wedge, below x1. */
+    inline const Vector& x4() const { return right_.x2(); }
+
+    /** Returns the six point of the wedge, below x2. */
+    inline const Vector& x5() const { return right_.x3(); }
+
+    /** Returns the top quad. */
+    inline const Quad& top() const { return top_; }
+
+    /** Returns the right quad. */
+    inline const Quad& right() const { return right_; }
+
+    /** Returns the bottom quad. */
+    inline const Quad& bottom() const { return bottom_; }
+
+    /** Returns the front triangle. */
+    inline const Triangle& front() const { return front_; }
+
+    /** Returns the back triangle. */
+    inline const Triangle& back() const { return back_; }
+
+    virtual bool intersects(const Sphere& s) const;
+    virtual Projection project(const Vector& p) const;
+
+  private:
+    Quad top_;
+    Quad right_;
+    Quad bottom_;
+    Triangle front_;
+    Triangle back_;
+  };
+
+  /**
+   * Represents an axis-aligned bounding box (AABB) using two points. The min
+   * point is the bottom back left point, and the max point is the top front
+   * right point.
+   */
+  class AxisAlignedBox : public Shape {
+  public:
+    AxisAlignedBox();
+    AxisAlignedBox(const Vector& min, const Vector& max);
+
+    /** Returns the min point, x0. */
+    inline const Vector& min() const { return min_; }
+
+    /** Returns the max point, x7. */
+    inline const Vector& max() const { return max_; }
+
+    /** Returns the left bottom back point. */
+    inline const Vector& x0() const { return min_; }
+
+    /** Returns the right bottom back point. */
+    inline Vector x1() const { return Vector(max_.x, min_.y, min_.z); }
+
+    /** Returns the right bottom front point. */
+    inline Vector x2() const { return Vector(max_.x, min_.y, max_.z); }
+
+    /** Returns the left bottom front point. */
+    inline Vector x3() const { return Vector(min_.x, min_.y, max_.z); }
+
+    /** Returns the right top back point. */
+    inline Vector x4() const { return Vector(max_.x, max_.y, min_.z); }
+
+    /** Returns the left top back point. */
+    inline Vector x5() const { return Vector(min_.x, max_.y, min_.z); }
+
+    /** Returns the left top front point. */
+    inline Vector x6() const { return Vector(min_.x, max_.y, max_.z); }
+
+    /** Returns the right top front point. */
+    inline const Vector& x7() const { return max_; }
+
+    /** Returns true if this box contains the specified point. */
+    bool contains(const Vector& p) const;
+
+    virtual bool intersects(const Sphere& s) const;
+    virtual Projection project(const Vector& p) const;
+
+  private:
+    Projection projectOut(const Vector& p) const;
+    Projection projectIn(const Vector& p) const;
+
+    Vector min_;
+    Vector max_;
+  };
+
+  /**
+   * Represents an oriented bounding box (OBB) using eight points.
+   */
+  class Box : public Shape {
+  public:
+    Box();
+    Box(const AxisAlignedBox& aab);
+    Box(const Vector& c, const Vector& x, const Vector& y, const Vector& z);
+    Box(const Vector& x0, const Vector& x1, const Vector& x2, const Vector& x3,
+        const Vector& x4, const Vector& x5, const Vector& x6, const Vector& x7);
+
+    /** Returns the left bottom back point. */
+    inline const Vector& x0() const { return bottom_.x0(); }
+
+    /** Returns the right bottom back point. */
+    inline const Vector& x1() const { return bottom_.x1(); }
+
+    /** Returns the right bottom front point. */
+    inline const Vector& x2() const { return bottom_.x2(); }
+
+    /** Returns the left bottom front point. */
+    inline const Vector& x3() const { return bottom_.x3(); }
+
+    /** Returns the right top back point. */
+    inline const Vector& x4() const { return top_.x0(); }
+
+    /** Returns the left top back point. */
+    inline const Vector& x5() const { return top_.x1(); }
+
+    /** Returns the left top front point. */
+    inline const Vector& x6() const { return top_.x2(); }
+
+    /** Returns the right top front point. */
+    inline const Vector& x7() const { return top_.x3(); }
+
+    /** Returns the left quad. */
+    inline const Quad& left() const { return back_; }
+
+    /** Returns the right quad. */
+    inline const Quad& right() const { return right_; }
+
+    /** Returns the bottom quad. */
+    inline const Quad& bottom() const { return bottom_; }
+
+    /** Returns the top quad. */
+    inline const Quad& top() const { return top_; }
+
+    /** Returns the back quad. */
+    inline const Quad& back() const { return back_; }
+
+    /** Returns the front quad. */
+    inline const Quad& front() const { return front_; }
+
+    virtual bool intersects(const Sphere& s) const;
+    virtual Projection project(const Vector& p) const;
+
+  private:
+    Quad top_;
+    Quad bottom_;
+    Quad left_;
+    Quad right_;
+    Quad front_;
+    Quad back_;
+  };
+
+  /** Represents a cylinder as two points and a radius. */
+  class Cylinder : public Shape {
+  public:
+    Cylinder();
+    Cylinder(const Vector& x0, const Vector& x1, float radius);
+
+    /** Returns the start point of the cylinder's main axis. */
+    inline const Vector& x0() const { return axis_.x0(); }
+
+    /** Returns the end point of the cylinder's main axis. */
+    inline const Vector& x1() const { return axis_.x1(); }
+
+    /** Returns the radius. */
+    inline float radius() const { return r_; }
+
+    /** Returns the length of the cylinder. */
+    inline float length() const { return l_; }
+
+    /** Returns the unit direction vector of the cylinder. */
+    inline const Vector& z() const { return v_; }
+
+    virtual bool intersects(const Sphere& s) const;
+    virtual Projection project(const Vector& p) const;
+
+  private:
+    LineSegment axis_;
+    float l_;
+    float r_;
+    Vector v_;
+  };
+
+};
+
+#endif

src/physics/shape_test.cpp

@@ -0,0 +1,301 @@
+// -*- C++ -*-
+
+#include <iostream>
+#include <math.h>
+#include <stdarg.h>
+#include <stdio.h>
+
+#include "shape.h"
+
+using namespace mbostock;
+
+static int returnCode = 0;
+
+void assertTrue(bool condition, char* message, ...) {
+  if (!condition) {
+    va_list args;
+    va_start(args, message);
+    printf("assertion failed: ");
+    vprintf(message, args);
+    printf("\n");
+    va_end(args);
+  }
+}
+
+static void testLineXYZ() {
+  printf("testLineXYZ...\n");
+  LineSegment x(Vector(0, 0, 0), Vector(2, 0, 0));
+  Projection px = x.project(Vector(1, 1, 0));
+  assertTrue(px.length == 1, "px.length");
+  assertTrue(px.x == Vector(1, 0, 0), "px.x");
+
+  LineSegment y(Vector(1, 1, 1), Vector(1, 4, 1));
+  Projection py = y.project(Vector(2, 1, 1));
+  assertTrue(py.length == 1, "py.length");
+  assertTrue(py.x == Vector(1, 1, 1), "py.x");
+
+  LineSegment z(Vector(-2, -2, -2), Vector(-2, -2, -4));
+  Projection pz = z.project(Vector(-2, -2, -2));
+  assertTrue(pz.length == 0, "pz.length");
+  assertTrue(pz.x == Vector(-2, -2, -2), "pz.x");
+}
+
+static void testLineX() {
+  printf("testLineX...\n");
+  LineSegment x(Vector(0, 0, 0), Vector(2, 0, 0));
+
+  /* Colinear, before x0. */
+  Projection p0 = x.project(Vector(-1, 0, 0));
+  assertTrue(p0.length == 1, "p0.length");
+  assertTrue(p0.x == Vector(0, 0, 0), "p0.x");
+
+  /* Colinear, after x1. */
+  Projection p1 = x.project(Vector(3, 0, 0));
+  assertTrue(p1.length == 1, "p1.length");
+  assertTrue(p1.x == Vector(2, 0, 0), "p1.x");
+
+  /* Colinear, at x0. */
+  Projection p2 = x.project(Vector(0, 0, 0));
+  assertTrue(p2.length == 0, "p2.length");
+  assertTrue(p2.x == Vector(0, 0, 0), "p2.x");
+
+  /* Colinear, at x1. */
+  Projection p3 = x.project(Vector(2, 0, 0));
+  assertTrue(p3.length == 0, "p3.length");
+  assertTrue(p3.x == Vector(2, 0, 0), "p3.x");
+
+  /* Colinear, between x0 and x1. */
+  for (float f = 0; f <= 2; f += .1f) {
+    Projection p4 = x.project(Vector(f, 0, 0));
+    assertTrue(p4.length == 0, "p4.length");
+    assertTrue(p4.x == Vector(f, 0, 0), "p4.x");
+  }
+
+  /* Noncolinear, before x0. */
+  Projection p5 = x.project(Vector(-1, -1, -1));
+  assertTrue(p5.length == sqrtf(3), "p5.length");
+  assertTrue(p5.x == Vector(0, 0, 0), "p5.x");
+
+  /* Noncolinear, after x1. */
+  Projection p6 = x.project(Vector(3, 1, 1));
+  assertTrue(p6.length == sqrtf(3), "p6.length");
+  assertTrue(p6.x == Vector(2, 0, 0), "p6.x");
+
+  /* Noncolinear, between x0 and x1. */
+  Projection p7 = x.project(Vector(1, 1, 1));
+  assertTrue(p7.length == sqrtf(2), "p7.length");
+  assertTrue(p7.x == Vector(1, 0, 0), "p7.x");
+}
+
+static void testPlaneXYZ() {
+  printf("testPlaneXYZ...\n");
+  Plane xz(Vector(0, 0, 0), Vector(0, 1, 0));
+  Projection pxz = xz.project(Vector(1, 2, 3));
+  assertTrue(pxz.length == 2, "pxz.length");
+  assertTrue(pxz.x == Vector(1, 0, 3), "pxz.x");
+
+  Plane xy(Vector(0, 0, 0), Vector(0, 0, 1));
+  Projection pxy = xy.project(Vector(1, 2, 3));
+  assertTrue(pxy.length == 3, "pxy.length");
+  assertTrue(pxy.x == Vector(1, 2, 0), "pxy.x");
+
+  Plane yz(Vector(0, 0, 0), Vector(1, 0, 0));
+  Projection pyz = yz.project(Vector(1, 2, 3));
+  assertTrue(pyz.length == 1, "pyz.length");
+  assertTrue(pyz.x == Vector(0, 2, 3), "pyz.x");
+
+  Projection pyz2 = yz.project(Vector(-1, 2, 3));
+  assertTrue(pyz2.length == 1, "pyz2.length");
+  assertTrue(pyz2.x == Vector(0, 2, 3), "pyz2.x");
+}
+
+static void testTriangle() {
+  printf("testTriangle...\n");
+
+  Triangle xy(Vector(1, 0, 0), Vector(0, 1, 0), Vector(0, 0, 0));
+  assertTrue(xy.normal() == Vector(0, 0, 1), "xy.normal");
+
+  /* At x0. */
+  Projection p0 = xy.project(Vector(0, 0, 0));
+  assertTrue(p0.length == 0, "p0.length");
+  assertTrue(p0.x == Vector(0, 0, 0), "p0.x");
+
+  /* At x1. */
+  Projection p1 = xy.project(Vector(0, 1, 0));
+  assertTrue(p1.length == 0, "p1.length");
+  assertTrue(p1.x == Vector(0, 1, 0), "p1.x");
+
+  /* At x2. */
+  Projection p2 = xy.project(Vector(1, 0, 0));
+  assertTrue(p2.length == 0, "p2.length");
+  assertTrue(p2.x == Vector(1, 0, 0), "p2.x");
+
+  /* In the plane of the triangle, x < x0. */
+  Projection p3 = xy.project(Vector(-1, 0, 0));
+  assertTrue(p3.length == 1, "p3.length");
+  assertTrue(p3.x == Vector(0, 0, 0), "p3.x");
+
+  /* In the plane of the triangle, x > x2. */
+  Projection p4 = xy.project(Vector(2, 0, 0));
+  assertTrue(p4.length == 1, "p4.length");
+  assertTrue(p4.x == Vector(1, 0, 0), "p4.x");
+
+  /* In the plane of the triangle, x = x2, y = x1. */
+  Projection p5 = xy.project(Vector(1, 1, 0));
+  assertTrue(p5.length == sqrtf(2) / 2, "p5.length");
+  assertTrue(p5.x == Vector(.5, .5, 0), "p5.x");
+
+  /* In front of the triangle (+z). */
+  Projection p6 = xy.project(Vector(.3, .3, 1));
+  assertTrue(p6.length == 1, "p6.length");
+  assertTrue(p6.x == Vector(.3, .3, 0), "p6.x");
+}
+
+static void testQuad() {
+  printf("testQuad...\n");
+
+  Quad xy(Vector(0, 0, 0), Vector(1, 0, 0), Vector(1, 1, 0), Vector(0, 1, 0));
+  assertTrue(xy.normal() == Vector(0, 0, 1), "xy.normal");
+
+  /* At x0. */
+  Projection p0 = xy.project(Vector(0, 0, 0));
+  assertTrue(p0.length == 0, "p0.length");
+  assertTrue(p0.x == Vector(0, 0, 0), "p0.x");
+
+  /* At x1. */
+  Projection p1 = xy.project(Vector(0, 1, 0));
+  assertTrue(p1.length == 0, "p1.length");
+  assertTrue(p1.x == Vector(0, 1, 0), "p1.x");
+
+  /* At x2. */
+  Projection p2 = xy.project(Vector(1, 0, 0));
+  assertTrue(p2.length == 0, "p2.length");
+  assertTrue(p2.x == Vector(1, 0, 0), "p2.x");
+
+  /* In the plane of the quad, x < x0. */
+  Projection p3 = xy.project(Vector(-1, 0, 0));
+  assertTrue(p3.length == 1, "p3.length");
+  assertTrue(p3.x == Vector(0, 0, 0), "p3.x");
+
+  /* In the plane of the quad, x > x2. */
+  Projection p4 = xy.project(Vector(2, 0, 0));
+  assertTrue(p4.length == 1, "p4.length");
+  assertTrue(p4.x == Vector(1, 0, 0), "p4.x");
+
+  /* In the plane of the quad, x = x2, y = x1. */
+  Projection p5 = xy.project(Vector(2, 2, 0));
+  assertTrue(p5.length == sqrtf(2), "p5.length");
+  assertTrue(p5.x == Vector(1, 1, 0), "p5.x");
+
+  /* In the plane of the quad. */
+  Projection p6 = xy.project(Vector(.5, .5, 0));
+  assertTrue(p6.length == 0, "p6.length");
+  assertTrue(p6.x == Vector(.5, .5, 0), "p6.x");
+
+  /* In front of the quad (+z). */
+  Projection p7 = xy.project(Vector(.3, .3, 1));
+  assertTrue(p7.length == 1, "p7.length");
+  assertTrue(p7.x == Vector(.3, .3, 0), "p6.x");
+}
+
+static void testRamp() {
+  printf("testRamp...\n");
+
+  Wedge r(Vector(-1, 0, .5),
+          Vector(-.5, .25, .5),
+          Vector(-.5, .25, -.5),
+          Vector(-1, 0, -.5));
+
+  Projection p0 = r.project(Vector(-0.694026, 0.05, 0.35));
+  assertTrue(p0.length > 0, "p0.length");
+}
+
+static void testCylinderIntersects() {
+  std::cout << "testCylinderIntersects...\n";
+  Cylinder c(Vector(0, 0, 0), Vector(0, 2, 0), 1);
+
+  Sphere s1(Vector(0, 3.1, 0), 1);
+  assertTrue(!c.intersects(s1), "c.intersects(s1)");
+
+  Sphere s2(Vector(0, 2.9, 0), 1);
+  assertTrue(c.intersects(s2), "c.intersects(s2)");
+
+  Sphere s3(Vector(0, 3, 0), .9);
+  assertTrue(!c.intersects(s3), "c.intersects(s3)");
+
+  Sphere s4(Vector(0, 3, 0), 1.1);
+  assertTrue(c.intersects(s4), "c.intersects(s4)");
+
+  Sphere s5(Vector(0, -1.1, 0), 1);
+  assertTrue(!c.intersects(s5), "c.intersects(s5)");
+
+  Sphere s6(Vector(0, -.9, 0), 1);
+  assertTrue(c.intersects(s6), "c.intersects(s6)");
+
+  Sphere s7(Vector(0, -1, 0), .9);
+  assertTrue(!c.intersects(s7), "c.intersects(s7)");
+
+  Sphere s8(Vector(0, -1, 0), 1.1);
+  assertTrue(c.intersects(s8), "c.intersects(s8)");
+
+  Sphere s9(Vector(2, 1, 0), .9);
+  assertTrue(!c.intersects(s9), "c.intersects(s9)");
+
+  Sphere s10(Vector(2, 1, 0), 1.1);
+  assertTrue(c.intersects(s10), "c.intersects(s10)");
+
+  Sphere s11(Vector(-2, 1, 0), .9);
+  assertTrue(!c.intersects(s11), "c.intersects(s11)");
+
+  Sphere s12(Vector(-2, 1, 0), 1.1);
+  assertTrue(c.intersects(s12), "c.intersects(s12)");
+
+  Sphere s13(Vector(0, 1, 2), .9);
+  assertTrue(!c.intersects(s13), "c.intersects(s13)");
+
+  Sphere s14(Vector(0, 1, 2), 1.1);
+  assertTrue(c.intersects(s14), "c.intersects(s14)");
+
+  Sphere s15(Vector(0, 1, -2), .9);
+  assertTrue(!c.intersects(s15), "c.intersects(s15)");
+
+  Sphere s16(Vector(0, 1, -2), 1.1);
+  assertTrue(c.intersects(s16), "c.intersects(s16)");
+}
+
+static void testCylinderProject() {
+  std::cout << "testCylinderProject...\n";
+  Cylinder c(Vector(0, 0, 0), Vector(0, 2, 0), 1);
+
+  Projection j1 = c.project(Vector(0, 3.1, 0));
+  assertTrue(Vector(0, 2, 0) == j1.x, "c.project(p1).x");
+  assertTrue(j1.length > 1.f && j1.length < 1.2f, "c.project(p1).length");
+  assertTrue(Vector(0, 1, 0) == j1.normal, "c.project(p1).normal");
+
+  Projection j2 = c.project(Vector(0, -1.1, 0));
+  assertTrue(Vector(0, 0, 0) == j2.x, "c.project(p2).x");
+  assertTrue(j2.length > 1.f && j2.length < 1.2f, "c.project(p2).length");
+  assertTrue(Vector(0, -1, 0) == j2.normal, "c.project(p2).normal");
+
+  Projection j3 = c.project(Vector(2.1, 1, 0));
+  assertTrue(Vector(1, 1, 0) == j3.x, "c.project(p3).x");
+  assertTrue(j3.length > 1.f && j3.length < 1.2f, "c.project(p3).length");
+  assertTrue(Vector(1, 0, 0) == j3.normal, "c.project(p3).normal");
+
+  Projection j4 = c.project(Vector(-2.1, 1, 0));
+  assertTrue(Vector(-1, 1, 0) == j4.x, "c.project(p4).x");
+  assertTrue(j4.length > 1.f && j4.length < 1.2f, "c.project(p4).length");
+  assertTrue(Vector(-1, 0, 0) == j4.normal, "c.project(p4).normal");
+}
+
+int main(int argc, char** argv) {
+  testLineXYZ();
+  testLineX();
+  testPlaneXYZ();
+  testTriangle();
+  testQuad();
+  testRamp();
+  testCylinderIntersects();
+  testCylinderProject();
+  return returnCode;
+}

src/physics/transform.cpp

@@ -0,0 +1,13 @@
+// -*- C++ -*-
+
+#include "transform.h"
+
+using namespace mbostock;
+
+Transform::Transform()
+    : enabled_(true) {
+}
+
+void Transform::enable(bool enabled) {
+  enabled_ = enabled;
+}

src/physics/transform.h

@@ -0,0 +1,25 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_TRANSFORM_H
+#define MBOSTOCK_TRANSFORM_H
+
+namespace mbostock {
+
+  class Transform {
+  public:
+    Transform();
+    virtual ~Transform() {}
+
+    void enable(bool enabled = true);
+    inline bool enabled() const { return enabled_; }
+
+    virtual void reset() = 0;
+    virtual void step() = 0;
+
+  private:
+    bool enabled_;
+  };
+
+}
+
+#endif

src/physics/translation.cpp

@@ -0,0 +1,87 @@
+// -*- C++ -*-
+
+#include "particle.h"
+#include "translation.h"
+#include "vector.h"
+
+using namespace mbostock;
+
+Translation::Translation(const Vector& x0, const Vector& x1,
+                         float speed, float start, float dampen)
+    : x0_(x0), x1_(x1), u_(start), s_(speed * ParticleSimulator::timeStep()),
+      kd_(1.f - dampen), v_((x1 - x0) * s_),
+      x_(x0 * (1.f - u_) + x1 * u_), origin_(x_),
+      mode_(REVERSE), reversed_(false) {
+}
+
+void Translation::reset() {
+  origin_ = x_ = x0_ * (1.f - u_) + x1_ * u_;
+  v_ = (x1_ - x0_) * s_;
+  reversed_ = false;
+}
+
+void Translation::step() {
+  if (!enabled()) {
+    return;
+  }
+  Vector x = origin_ + v_;
+  if (reversed_) {
+    if (v_.dot(x0_ - x) < 0.f) {
+      reversed_ = false;
+      v_ *= -1.f;
+    }
+  } else {
+    if (v_.dot(x1_ - x) < 0.f) {
+      switch (mode_) {
+        case REVERSE: {
+          reversed_ = true;
+          v_ *= -1.f;
+          break;
+        }
+        case RESET: {
+          x_ = x0_;
+          origin_ = x0_;
+          break;
+        }
+        case ONE_WAY: {
+          v_ = Vector::ZERO();
+          break;
+        }
+      }
+    }
+  }
+  x_ += v_;
+  dv_ = (x_ - origin_) * kd_;
+  origin_ += dv_;
+}
+
+void Translation::setMode(Mode m) {
+  mode_ = m;
+}
+
+const Vector& Translation::velocity() const {
+  return enabled() ? dv_ : Vector::ZERO();
+}
+
+Vector Translation::point(const Vector& x) const {
+  return x + origin_;
+}
+
+Vector Translation::pointInverse(const Vector& x) const {
+  return x - origin_;
+}
+
+TranslatingShape::TranslatingShape(const Shape& s, const Translation& t)
+    : shape_(s), translation_(t) {
+}
+
+bool TranslatingShape::intersects(const Sphere& s) const {
+  Sphere ts(translation_.pointInverse(s.x()), s.radius());
+  return shape_.intersects(ts);
+}
+
+Projection TranslatingShape::project(const Vector& x) const {
+  Projection p = shape_.project(translation_.pointInverse(x));
+  p.x = translation_.point(p.x);
+  return p;
+}

src/physics/translation.h

@@ -0,0 +1,70 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_TRANSLATION_H
+#define MBOSTOCK_TRANSLATION_H
+
+#include "shape.h"
+#include "transform.h"
+#include "vector.h"
+
+namespace mbostock {
+
+  class Translation : public Transform {
+  public:
+    Translation(const Vector& x0, const Vector& x1,
+                float s, float u, float kd);
+
+    enum Mode { REVERSE, RESET, ONE_WAY };
+
+    /* Sets the translation mode. The default is REVERSE. */
+    void setMode(Mode m);
+
+    /** Translates the specified point. */
+    Vector point(const Vector& x) const;
+
+    /** Inverse-translates the specified point. */
+    Vector pointInverse(const Vector& x) const;
+
+    /** Advances the translation by one time step. */
+    virtual void step();
+
+    /** Resets the translation. */
+    virtual void reset();
+
+    /** Returns the current velocity. */
+    const Vector& velocity() const;
+
+    /** Returns the current origin. */
+    inline const Vector& origin() const { return origin_; }
+
+  private:
+    void update();
+
+    Vector x0_;
+    Vector x1_;
+    float s_;
+    float u_;
+    float kd_;
+    Vector v_;
+    Vector dv_;
+    Vector x_;
+    Vector origin_;
+    Mode mode_;
+    bool reversed_;
+  };
+
+  class TranslatingShape : public Shape {
+  public:
+    TranslatingShape(const Shape& s, const Translation& t);
+
+    virtual bool intersects(const Sphere& s) const;
+    virtual Projection project(const Vector& x) const;
+
+  private:
+    const Shape& shape_;
+    const Translation& translation_;
+  };
+
+}
+
+#endif

src/physics/vector.cpp

@@ -0,0 +1,52 @@
+// -*- C++ -*-
+
+#include <algorithm>
+#include <stdlib.h>
+
+#include "vector.h"
+
+using namespace mbostock;
+
+const Vector& Vector::ZERO() {
+  static const Vector v;
+  return v;
+}
+
+const Vector& Vector::X() {
+  static const Vector v(1.f, 0.f, 0.f);
+  return v;
+}
+
+const Vector& Vector::Y() {
+  static const Vector v(0.f, 1.f, 0.f);
+  return v;
+}
+
+const Vector& Vector::Z() {
+  static const Vector v(0.f, 0.f, 1.f);
+  return v;
+}
+
+const Vector& Vector::INF() {
+  static const Vector v(INFINITY, INFINITY, INFINITY);
+  return v;
+}
+
+static float randomf() {
+  return random() / (float) RAND_MAX;
+}
+
+Vector Vector::min(const Vector& a, const Vector& b) {
+  return Vector(std::min(a.x, b.x), std::min(a.y, b.y), std::min(a.z, b.z));
+}
+
+Vector Vector::max(const Vector& a, const Vector& b) {
+  return Vector(std::max(a.x, b.x), std::max(a.y, b.y), std::max(a.z, b.z));
+}
+
+Vector Vector::randomVector(float k) {
+  Vector v(randomf() - .5f, randomf() - .5f, randomf() - .5f);
+  return ((v.x == 0.f) && (v.y == 0.f) && (v.z == 0.f))
+    ? Vector(k, 0.f, 0.f) // very unlikely, but still...
+    : v.normalize() * k;
+}

src/physics/vector.h

@@ -0,0 +1,134 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_VECTOR_H
+#define MBOSTOCK_VECTOR_H
+
+#include "math.h"
+
+namespace mbostock {
+
+  class Vector {
+  public:
+    inline Vector() : x(0.f), y(0.f), z(0.f) {}
+    inline Vector(float x, float y) : x(x), y(y), z(0.f) {}
+    inline Vector(float x, float y, float z) : x(x), y(y), z(z) {}
+
+    inline Vector operator*(float k) const {
+      return Vector(x * k, y * k, z * k);
+    }
+
+    inline Vector operator/(float k) const {
+      return Vector(x / k, y / k, z / k);
+    }
+
+    inline Vector operator-(float k) const {
+      return Vector(x - k, y - k, z - k);
+    }
+
+    inline Vector operator+(float k) const {
+      return Vector(x + k, y + k, z + k);
+    }
+
+    inline Vector operator-() const {
+      return Vector(-x, -y, -z);
+    }
+
+    inline Vector& operator*=(float k) {
+      x *= k; y *= k; z *= k;
+      return *this;
+    }
+    inline Vector& operator/=(float k) {
+      x /= k; y /= k; z /= k;
+      return *this;
+    }
+    inline Vector& operator-=(float k) {
+      x -= k; y -= k; z -= k;
+      return *this;
+    }
+    inline Vector& operator+=(float k) {
+      x += k; y += k; z += k;
+      return *this;
+    }
+
+    inline Vector operator-(const Vector& v) const {
+      return Vector(x - v.x, y - v.y, z - v.z);
+    }
+
+    inline Vector operator+(const Vector& v) const {
+      return Vector(x + v.x, y + v.y, z + v.z);
+    }
+
+    inline Vector operator*(const Vector& v) const {
+      return Vector(x * v.x, y * v.y, z * v.z);
+    }
+
+    inline Vector operator/(const Vector& v) const {
+      return Vector(x / v.x, y / v.y, z / v.z);
+    }
+
+    inline Vector& operator-=(const Vector& v) {
+      x -= v.x; y -= v.y; z -= v.z;
+      return *this;
+    }
+
+    inline Vector& operator+=(const Vector& v) {
+      x += v.x; y += v.y; z += v.z;
+      return *this;
+    }
+
+    inline Vector& operator*=(const Vector& v) {
+      x *= v.x; y *= v.y; z *= v.z;
+      return *this;
+    }
+
+    inline Vector& operator/=(const Vector& v) {
+      x /= v.x; y /= v.y; z /= v.z;
+      return *this;
+    }
+
+    inline bool operator==(const Vector& v) const {
+      return (x == v.x) && (y == v.y) && (z == v.z);
+    }
+
+    inline bool operator!=(const Vector& v) const {
+      return (x != v.x) || (y != v.y) || (z != v.z);
+    }
+
+    inline float dot(const Vector& v) const {
+      return x * v.x + y * v.y + z * v.z;
+    }
+
+    inline Vector cross(const Vector& v) const {
+      return Vector(y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x);
+    }
+
+    inline float length() const {
+      return sqrtf(x * x + y * y + z * z);
+    }
+
+    inline float squared() const {
+      return x * x + y * y + z * z;
+    }
+
+    inline Vector normalize() const {
+      return *this / length();
+    }
+
+    static Vector min(const Vector& a, const Vector& b);
+    static Vector max(const Vector& a, const Vector& b);
+    static Vector randomVector(float k);
+
+    static const Vector& ZERO();
+    static const Vector& X();
+    static const Vector& Y();
+    static const Vector& Z();
+    static const Vector& INF();
+
+    float x;
+    float y;
+    float z;
+  };
+
+}
+
+#endif

src/physics/vector_test.cpp

@@ -0,0 +1,240 @@
+// -*- C++ -*-
+
+#include <stdio.h>
+#include <stdarg.h>
+
+#include "vector.h"
+
+using namespace mbostock;
+
+static int returnCode = 0;
+
+void assertTrue(bool condition, char* message, ...) {
+  if (!condition) {
+    va_list args;
+    va_start(args, message);
+    printf("assertion failed: ");
+    vprintf(message, args);
+    printf("\n");
+    va_end(args);
+  }
+}
+
+static void testDefaultConstructor() {
+  printf("testDefaultConstructor...\n");
+  Vector v1;
+  Vector v2;
+  assertTrue(v1 == v2, "v1 != v2");
+  assertTrue(v1.x == 0.f, "v1.x != 0.f");
+  assertTrue(v1.y == 0.f, "v1.y != 0.f");
+  assertTrue(v1.z == 0.f, "v1.z != 0.f");
+}
+
+static void testExplicitConstructor() {
+  printf("testExplicitConstructor...\n");
+  Vector v(1.f, 2.f, 3.f);
+  assertTrue(v.x == 1.f, "v.x != 1.f");
+  assertTrue(v.y == 2.f, "v.y != 2.f");
+  assertTrue(v.z == 3.f, "v.z != 3.f");
+}
+
+static void testAssignment() {
+  printf("testAssignment...\n");
+  Vector v1(1.f, 2.f, 3.f);
+  Vector v2(4.f, 5.f, 6.f);
+  Vector v3(4.f, 5.f, 6.f);
+  v1 = v2;
+  assertTrue(v1 == v2, "v1 != v2");
+  assertTrue(v1 == v3, "v1 != v3");
+  assertTrue(v2 == v3, "v2 != v3");
+}
+
+static void testCopyConstructor() {
+  printf("testCopyConstructor...\n");
+  Vector v1(1.f, 2.f, 3.f);
+  Vector v2(v1);
+  assertTrue(v1 == v2, "v1 != v2");
+  assertTrue(v2 == v1, "v2 != v1");
+}
+
+static void testEqualVector() {
+  printf("testEqualVector...\n");
+  Vector v1a(1.f, 2.f, 3.f);
+  Vector v1b(1.f, 2.f, 3.f);
+  Vector v2(1.f, 2.f, 4.f);
+  Vector v3(1.f, 3.f, 3.f);
+  Vector v4(2.f, 3.f, 3.f);
+  assertTrue(v1a == v1b, "v1a != v1b");
+  assertTrue(v1b == v1a, "v1b != v1a");
+  assertTrue(!(v1a == v2), "v1a == v2");
+  assertTrue(!(v1a == v3), "v1a == v3");
+  assertTrue(!(v1a == v4), "v1a == v4");
+}
+
+static void testNotEqualVector() {
+  printf("testNotEqualVector...\n");
+  Vector v1a(1.f, 2.f, 3.f);
+  Vector v1b(1.f, 2.f, 3.f);
+  Vector v2(1.f, 2.f, 4.f);
+  Vector v3(1.f, 3.f, 3.f);
+  Vector v4(2.f, 3.f, 3.f);
+  assertTrue(!(v1a != v1b), "v1a != v1b");
+  assertTrue(!(v1b != v1a), "v1b != v1a");
+  assertTrue(v1a != v2, "v1a == v2");
+  assertTrue(v1a != v3, "v1a == v3");
+  assertTrue(v1a != v4, "v1a == v4");
+}
+
+static void testMultiplyScalar() {
+  printf("testMultiplyScalar...\n");
+  Vector v1(1.f, 2.f, 3.f);
+  Vector v2(2.f, 4.f, 6.f);
+  assertTrue(v2 == (v1 * 2), "v2 != v1 * 2");
+  assertTrue(v1 == (v2 * .5), "v1 != v2 * .5");
+}
+
+static void testDivideScalar() {
+  printf("testDivideScalar...\n");
+  Vector v1(1.f, 2.f, 3.f);
+  Vector v2(2.f, 4.f, 6.f);
+  assertTrue(v2 == (v1 / .5), "v2 != v1 / .5");
+  assertTrue(v1 == (v2 / 2), "v1 != v2 / 2");
+}
+
+static void testAddScalar() {
+  printf("testAddScalar...\n");
+  Vector v1(1, 2, 3);
+  assertTrue(Vector(2, 3, 4) == (v1 + 1), "v1 + 1");
+  assertTrue(Vector(0, 1, 2) == (v1 + -1), "v1 + -1");
+}
+
+static void testSubtractScalar() {
+  printf("testSubtractScalar...\n");
+  Vector v1(1, 2, 3);
+  assertTrue(Vector(0, 1, 2) == (v1 - 1), "v1 - 1");
+  assertTrue(Vector(2, 3, 4) == (v1 - -1), "v1 - -1");
+}
+
+static void testMultiplyAssignmentScalar() {
+  printf("testMultiplyAssignmentScalar...\n");
+  Vector v1(1, 2, 3);
+  v1 *= 2;
+  assertTrue(Vector(2, 4, 6) == v1, "v1 * 2");
+  v1 *= .5;
+  assertTrue(Vector(1, 2, 3) == v1, "v1 * .5");
+}
+
+static void testDivideAssignmentScalar() {
+  printf("testDivideAssignmentScalar...\n");
+  Vector v1(1, 2, 3);
+  v1 /= .5;
+  assertTrue(Vector(2, 4, 6) == v1, "v1 / .5");
+  v1 /= 2;
+  assertTrue(Vector(1, 2, 3) == v1, "v1 / 2");
+}
+
+static void testAddAssignmentScalar() {
+  printf("testAddAssignmentScalar...\n");
+  Vector v1(1, 2, 3);
+  v1 += 1;
+  assertTrue(Vector(2, 3, 4) == v1, "v1 + 1");
+  v1 += -1;
+  assertTrue(Vector(1, 2, 3) == v1, "v1 + -1");
+}
+
+static void testSubtractAssignmentScalar() {
+  printf("testSubtractAssignmentScalar...\n");
+  Vector v1(1, 2, 3);
+  v1 -= -1;
+  assertTrue(Vector(2, 3, 4) == v1, "v1 - -1");
+  v1 -= 1;
+  assertTrue(Vector(1, 2, 3) == v1, "v1 - 1");
+}
+
+static void testAddVector() {
+  printf("testAddVector...\n");
+  Vector v1(1, 2, 3);
+  Vector v2 = v1 * 2;
+  Vector v3 = v1 + v2;
+  assertTrue(Vector(3, 6, 9) == v3, "v1 + v1 * 2");
+}
+
+static void testSubtractVector() {
+  printf("testSubtractVector...\n");
+  Vector v1(1, 2, 3);
+  Vector v2 = v1 * 2;
+  Vector v3 = v1 - v2;
+  assertTrue(Vector(-1, -2, -3) == v3, "v1 - v1 * 2");
+}
+
+static void testAddAssignmentVector() {
+  printf("testAddAssignmentVector...\n");
+  Vector v1(1, 2, 3);
+  v1 += v1 * 2;
+  assertTrue(Vector(3, 6, 9) == v1, "v1 + v1 * 2");
+}
+
+static void testSubtractAssignmentVector() {
+  printf("testSubtractAssignmentVector...\n");
+  Vector v1(1, 2, 3);
+  v1 -= v1 * 2;
+  assertTrue(Vector(-1, -2, -3) == v1, "v1 - v1 * 2");
+}
+
+static void testNegative() {
+  printf("testNegative...\n");
+  Vector v1(1, 2, 3);
+  Vector v2(-1, -2, -3);
+  assertTrue(v2 == -v1, "v2 != -v1");
+}
+
+static void testDotVector() {
+  printf("testDotVector...\n");
+  Vector v1(1, 2, 3);
+  Vector v2(-2, -3, -4);
+  float d = 1 * -2 + 2 * -3 + 3 * -4;
+  assertTrue(v1.dot(v2) == d, "v1 dot v2");
+}
+
+static void testCrossVector() {
+  printf("testCrossVector...\n");
+  Vector x(1, 0, 0);
+  Vector y(0, 1, 0);
+  Vector z(0, 0, 1);
+  assertTrue(x.cross(y) == z, "x cross y");
+}
+
+static void testLength() {
+  printf("testLength...\n");
+}
+
+static void testNormalize() {
+  printf("testNormalize...\n");
+}
+
+int main(int argc, char** argv) {
+  testDefaultConstructor();
+  testExplicitConstructor();
+  testAssignment();
+  testCopyConstructor();
+  testEqualVector();
+  testNotEqualVector();
+  testMultiplyScalar();
+  testDivideScalar();
+  testAddScalar();
+  testSubtractScalar();
+  testMultiplyAssignmentScalar();
+  testDivideAssignmentScalar();
+  testAddAssignmentScalar();
+  testSubtractAssignmentScalar();
+  testAddVector();
+  testSubtractVector();
+  testAddAssignmentVector();
+  testSubtractAssignmentVector();
+  testNegative();
+  testDotVector();
+  testCrossVector();
+  testLength();
+  testNormalize();
+  return returnCode;
+}

src/player.cpp

@@ -0,0 +1,626 @@
+// -*- C++ -*-
+
+#include <GLUT/glut.h>
+#include <OpenGL/gl.h>
+#include <OpenGL/glu.h>
+#include <iostream>
+#include <math.h>
+#include <stdio.h>
+#include <vector>
+
+#include "material.h"
+#include "model.h"
+#include "physics/constraint.h"
+#include "physics/force.h"
+#include "physics/particle.h"
+#include "physics/shape.h"
+#include "physics/vector.h"
+#include "player.h"
+#include "room.h"
+#include "room_object.h"
+#include "world.h"
+
+using namespace mbostock;
+
+static const float axleRadius = .02f;
+static const float axleLength = .14f;
+static const float bodySize = .15f;
+static const float eyeRadius = .02f;
+static const float wheelRadius = .1f;
+static const float wheelWeight = 1.f;
+static const float counterWeightRadius = .07f;
+static const float counterWeightOffset = .07f;
+static const float counterWeight = .1f;
+
+static const float forwardForce = 12.f;
+static const float backwardForce = 5.f;
+static const float turnForce = 1.f;
+static const float motorFriction = 1.f;
+static const float brakeFriction = 3.f;
+
+static const Vector spokeColor(.5f, .5f, .5f);
+static const Vector eyeColor(.5f, .5f, .5f);
+static const Vector axleMaterialDiffuse(.9f, .8f, .8f);
+static const Vector tireMaterialDiffuse(.1f, .1f, .1f);
+static const Vector bodyMaterialDiffuse(.6f, .2f, .3f);
+
+namespace mbostock {
+
+  /**
+   * A small model which displays just the player's wheel, so that this model
+   * can be compiled into a display list for faster display. This model uses the
+   * parent model's quadric.
+   */
+  class PlayerWheelModel : public Model {
+  public:
+    PlayerWheelModel(const PlayerModel& model);
+
+    virtual void display();
+
+  private:
+    const PlayerModel& model_;
+  };
+
+  /**
+   * A small model which displays just the player's body (and axle, and eyes),
+   * so that this model can be compiled into a display list for faster display.
+   * This model uses the parent model's quadric.
+   */
+  class PlayerBodyModel : public Model {
+  public:
+    PlayerBodyModel(const PlayerModel& model);
+
+    virtual void display();
+
+  private:
+    const PlayerModel& model_;
+  };
+
+}
+
+PlayerWheelModel::PlayerWheelModel(const PlayerModel& model)
+    : model_(model) {
+}
+
+void PlayerWheelModel::display() {
+  /* Spokes. */
+  glDisable(GL_LIGHTING);
+  glColorv(spokeColor);
+  glLineWidth(2.f);
+  gluQuadricDrawStyle(model_.quadric_, GLU_LINE);
+  gluDisk(model_.quadric_, 0, wheelRadius - wheelRadius / 8.f, 8, 1);
+  gluQuadricDrawStyle(model_.quadric_, GLU_FILL);
+  glEnable(GL_LIGHTING);
+
+  /* Tire. */
+  glMaterialv(GL_FRONT_AND_BACK, GL_DIFFUSE, tireMaterialDiffuse);
+  glutSolidTorus(wheelRadius / 4.f, 3 * wheelRadius / 4.f, 16, 32);
+}
+
+PlayerBodyModel::PlayerBodyModel(const PlayerModel& model) 
+  : model_(model) {
+}
+
+void PlayerBodyModel::display() {
+  float l = axleLength + wheelRadius / 4.f + wheelRadius;
+
+  /* Axle endcap. */
+  glMaterialv(GL_FRONT_AND_BACK, GL_DIFFUSE, axleMaterialDiffuse);
+  glPushMatrix();
+  glTranslatef(0.f, 0.f, -l / 2.f);
+  glRotatef(180.f, 0.f, 1.f, 0.f);
+  gluDisk(model_.quadric_, 0.f, axleRadius, 16, 8);
+  glPopMatrix();
+
+  /* Axle. */
+  glPushMatrix();
+  glTranslatef(0.f, 0.f, -l / 2.f);
+  gluCylinder(model_.quadric_, axleRadius, axleRadius, l, 16, 1);
+  glPopMatrix();
+
+  /* Axle endcap. */
+  glPushMatrix();
+  glTranslatef(0.f, 0.f, l / 2.f);
+  gluDisk(model_.quadric_, 0.f, axleRadius, 16, 8);
+  glPopMatrix();
+
+  /* Body. */
+  glMaterialv(GL_FRONT_AND_BACK, GL_DIFFUSE, bodyMaterialDiffuse);
+  glutSolidCube(bodySize);
+
+  /* Eyes. */
+  static const float eyeDepth = .01f;
+  glMaterialv(GL_FRONT_AND_BACK, GL_DIFFUSE, axleMaterialDiffuse);
+  glColorv(eyeColor);
+  glTranslatef(bodySize / 2.f, bodySize / 6.f, bodySize / 4.f);
+  glPushMatrix();
+  glRotatef(90.f, 0.f, 1.f, 0.f);
+  gluCylinder(model_.quadric_, eyeRadius, eyeRadius, eyeDepth, 16, 1);
+  glTranslatef(0.f, 0.f, eyeDepth);
+  gluDisk(model_.quadric_, 0, eyeRadius, 8, 1);
+  glTranslatef(0.f, -.004f, .001f);
+  glMaterialv(GL_FRONT_AND_BACK, GL_DIFFUSE, tireMaterialDiffuse);
+  gluDisk(model_.quadric_, 0, eyeRadius / 2.f, 8, 1);
+  glMaterialv(GL_FRONT_AND_BACK, GL_DIFFUSE, axleMaterialDiffuse);
+  glPopMatrix();
+  glTranslatef(0.f, 0.f, -bodySize / 2.f);
+  glPushMatrix();
+  glRotatef(90.f, 0.f, 1.f, 0.f);
+  gluCylinder(model_.quadric_, eyeRadius, eyeRadius, eyeDepth, 16, 1);
+  glTranslatef(0.f, 0.f, eyeDepth);
+  gluDisk(model_.quadric_, 0, eyeRadius, 8, 1);
+  glTranslatef(0.f, -.004f, .001f);
+  glMaterialv(GL_FRONT_AND_BACK, GL_DIFFUSE, tireMaterialDiffuse);
+  gluDisk(model_.quadric_, 0, eyeRadius / 2.f, 8, 1);
+  glPopMatrix();
+}
+
+PlayerModel::PlayerModel(const Player& player)
+  : player_(player), quadric_(NULL),
+    wheelModel_(Models::compile(new PlayerWheelModel(*this))),
+    bodyModel_(Models::compile(new PlayerBodyModel(*this))) {
+  for (int i = 0; i < 15; i++) {
+    orientation_[i] = 0.f;
+  }
+  orientation_[15] = 1.f;
+}
+
+PlayerModel::~PlayerModel() {
+  delete wheelModel_;
+  delete bodyModel_;
+  if (quadric_ != NULL) {
+    gluDeleteQuadric(quadric_);
+  }
+}
+
+void PlayerModel::initialize() {
+  if (quadric_ != NULL) {
+    gluDeleteQuadric(quadric_);
+  }
+  quadric_ = gluNewQuadric();
+  wheelModel_->initialize();
+  bodyModel_->initialize();
+}
+
+float* PlayerModel::orientation() {
+  const Vector& x = player_.x();
+  const Vector& y = player_.y();
+  const Vector& z = player_.z();
+  orientation_[0] = x.x; orientation_[1] = x.y; orientation_[2] = x.z;
+  orientation_[4] = y.x; orientation_[5] = y.y; orientation_[6] = y.z;
+  orientation_[8] = z.x; orientation_[9] = z.y; orientation_[10] = z.z;
+  return orientation_;
+}
+
+void PlayerModel::display() {
+  Materials::blank().bind();
+
+  glPushMatrix();
+  glTranslatev(player_.origin());
+  glMultMatrixf(orientation());
+
+  /* Left wheel. */
+  if (!World::world()->debug() || player_.leftWheelFriction()) {
+    glPushMatrix();
+    glTranslatef(0.f, 0.f, -axleLength / 2.f - wheelRadius / 2.f);
+    glRotatef(player_.leftWheelAngle(), 0.f, 0.f, 1.f);
+    glRotatef(180.f, 0.f, 1.f, 0.f);
+    wheelModel_->display();
+    glPopMatrix();
+  }
+
+  /* Right wheel. */
+  if (!World::world()->debug() || player_.rightWheelFriction()) {
+    glPushMatrix();
+    glTranslatef(0.f, 0.f, axleLength / 2.f + wheelRadius / 2.f);
+    glRotatef(player_.rightWheelAngle(), 0.f, 0.f, 1.f);
+    wheelModel_->display();
+    glPopMatrix();
+  }
+
+  /* Axes. */
+  if (World::world()->debug()) {
+    displayAxes();
+  }
+
+  /* Body. */
+  bodyModel_->display();
+  glPopMatrix();
+}
+
+void PlayerModel::displayAxes() {
+  glDisable(GL_LIGHTING);
+  glBegin(GL_LINES);
+  glColor3f(1.f, 0.f, 0.f);
+  glVertex3f(0.f, 0.f, 0.f);
+  glVertex3f(2 * bodySize, 0.f, 0.f);
+  glColor3f(0.f, 1.f, 0.f);
+  glVertex3f(0.f, 0.f, 0.f);
+  glVertex3f(0.f, 2 * bodySize, 0.f);
+  glColor3f(0.f, 0.f, 1.f);
+  glVertex3f(0.f, 0.f, 0.f);
+  glVertex3f(0.f, 0.f, 2 * bodySize);
+  glEnd();
+  glEnable(GL_LIGHTING);
+}
+
+Player::Player()
+    : turnState_(NONE), moveState_(NONE),
+      sphere_(Vector::ZERO(), wheelRadius * 2.f), model_(*this) {
+  counterWeight_.inverseMass = 1.f / counterWeight;
+}
+
+float Player::mass() const {
+  return counterWeight + 3.f;
+}
+
+Player::Wheel::Wheel()
+    : angle(0.f), angleStep(0.f) {
+  inverseMass = 1.f / wheelWeight;
+}
+
+void Player::setOrigin(const Vector& origin) {
+  static const float a2 = axleLength / 2.f;
+
+  leftWheel_.position = origin + Vector(0.f, 0.f, -a2);
+  rightWheel_.position = origin + Vector(0.f, 0.f, a2);
+  body_.position = origin;
+  counterWeight_.position = origin + Vector(-counterWeightOffset, 0.f, 0.f);
+
+  leftWheel_.previousPosition = leftWheel_.position;
+  rightWheel_.previousPosition = rightWheel_.position;
+  body_.previousPosition = body_.position;
+  counterWeight_.previousPosition = counterWeight_.position;
+
+  origin_ = origin;
+  x_ = Vector::X();
+  y_ = Vector::Y();
+  z_ = Vector::Z();
+}
+
+void Player::setVelocity(const Vector& velocity) {
+
+  /* TODO: Rotate the player to face the velocity. But this works for now. */
+  if (velocity.x < 0) {
+    static const float a2 = axleLength / 2.f;
+
+    leftWheel_.position = origin_ + Vector(0.f, 0.f, a2);
+    rightWheel_.position = origin_ + Vector(0.f, 0.f, -a2);
+    counterWeight_.position = origin_ + Vector(counterWeightOffset, 0.f, 0.f);
+
+    leftWheel_.previousPosition = leftWheel_.position;
+    rightWheel_.previousPosition = rightWheel_.position;
+    counterWeight_.previousPosition = counterWeight_.position;
+  }
+
+  Vector v = velocity * ParticleSimulator::timeStep();
+  leftWheel_.previousPosition = leftWheel_.position - v;
+  rightWheel_.previousPosition = rightWheel_.position - v;
+  body_.previousPosition = body_.position - v;
+  counterWeight_.previousPosition = counterWeight_.position - v;
+}
+
+void Player::resetForces() {
+  leftWheel_.force = Vector::ZERO();
+  rightWheel_.force = Vector::ZERO();
+  body_.force = Vector::ZERO();
+  counterWeight_.force = Vector::ZERO();
+
+  /* Moving forces. */
+  switch (moveState_) {
+    case FORWARD: {
+      leftWheel_.applyForce(z_, forwardForce);
+      rightWheel_.applyForce(z_, forwardForce);
+      break;
+    }
+    case BACKWARD: {
+      leftWheel_.applyForce(z_, -backwardForce);
+      rightWheel_.applyForce(z_, -backwardForce);
+      break;
+    }
+  }
+
+  /* Turning forces. */
+  switch (turnState_) {
+    case RIGHT: {
+      leftWheel_.applyForce(z_, turnForce);
+      rightWheel_.applyForce(z_, -turnForce);
+      break;
+    }
+    case LEFT: {
+      leftWheel_.applyForce(z_, -turnForce);
+      rightWheel_.applyForce(z_, turnForce);
+      break;
+    }
+  }
+
+  /* Uniform drag forces. */
+  if ((turnState_ != NONE) || (moveState_ == FORWARD)) {
+    leftWheel_.applyQuadraticDrag(motorFriction);
+    rightWheel_.applyQuadraticDrag(motorFriction);
+  } else {
+    leftWheel_.applyLinearDrag(brakeFriction);
+    rightWheel_.applyLinearDrag(brakeFriction);
+  }
+
+  /* Flipping force if upside-down. */
+  if ((leftWheel_.contact && rightWheel_.contact)
+      && (y_.y < (.5 * leftWheel_.contactNormal.y))) {
+    counterWeight_.force -= y_;
+  }
+}
+
+void Player::Wheel::applyContact(const RoomObject& o) {
+  const Vector& normal = Constraints::projection().normal;
+  if ((normal.y > o.slip()) && (!contact || (contactNormal.y < normal.y))) {
+    contactNormal = normal;
+    contactVelocity = o.velocity(position);
+    contact = true;
+  }
+}
+
+void Player::Wheel::applyForce(const Vector& z, float f) {
+  if (friction()) {
+    force += z.cross(contactNormal) * -f;
+  }
+}
+
+void Player::Wheel::applyLinearDrag(float kd) {
+  if (friction()) {
+    Vector v = (position - previousPosition - contactVelocity)
+        / ParticleSimulator::timeStep();
+    float vs = v.squared();
+    if (vs < 1E-3) {
+      kd *= 10.f;
+    } else if (vs < 1E-2) {
+      kd *= 2.f;
+    }
+    force -= v * kd;
+  }
+}
+
+void Player::Wheel::applyQuadraticDrag(float kd) {
+  if (friction()) {
+    Vector v = (position - previousPosition - contactVelocity)
+        / ParticleSimulator::timeStep();
+    force -= v * v.length() * kd;
+  }
+}
+
+void Player::applyForce(UnaryForce& force) {
+  force.apply(leftWheel_);
+  force.apply(rightWheel_);
+  force.apply(body_);
+  force.apply(counterWeight_);
+}
+
+void Player::step(const ParticleSimulator& s) {
+  leftWheel_.contact = false;
+  rightWheel_.contact = false;
+  s.step(leftWheel_);
+  s.step(rightWheel_);
+  s.step(body_);
+  s.step(counterWeight_);
+  sphere_.x() = body_.position;
+}
+
+bool Player::intersects(const Shape& s) {
+  return s.intersects(sphere_);
+}
+
+bool Player::constrainOutside(const RoomObject& o) {
+  const Shape& s = o.shape();
+  bool contact = false;
+  if (s.intersects(sphere_)) {
+    if (Constraints::outside(leftWheel_, s, wheelRadius)) {
+      constrainGlancing(o);
+      leftWheel_.applyContact(o);
+      contact = true;
+    }
+    if (Constraints::outside(rightWheel_, s, wheelRadius)) { 
+      constrainGlancing(o);
+      rightWheel_.applyContact(o);
+      contact = true;
+    }
+    if (Constraints::outside(body_, s, wheelRadius)) {
+      sphere_.x() = body_.position;
+      contact = true;
+    }
+    if (Constraints::outside(counterWeight_, s, counterWeightRadius)) {
+      contact = true;
+    }
+  }
+  return contact;
+}
+
+void Player::constrainGlancing(const RoomObject& o) {
+  const Projection& j = Constraints::projection();
+
+  /*
+   * This constraint only applies if we are moving parallel to the wall, so if
+   * the body particle isn't moving, don't do anything.
+   */
+  if ((body_.position - body_.previousPosition).squared() < 1E-5) {
+    return;
+  }
+
+  /* For glancing collisions, align parallel to the wall. */
+  float n = j.normal.dot(z_);
+  if (fabsf(n) > .8f) {
+    Plane p = Plane(j.x, j.normal);
+    Vector x = p.project((leftWheel_.position + rightWheel_.position) / 2.f).x;
+
+    Vector wp1 = x + j.normal * wheelRadius;
+    Vector wp2 = x + j.normal * (wheelRadius + axleLength);
+    Vector bp = x + j.normal * (wheelRadius + axleLength / 2.f);
+    const Vector& lp = (n > 0.f) ? wp1 : wp2;
+    const Vector& rp = (n > 0.f) ? wp2 : wp1;
+
+    /*
+     * Only allow this constraint to move the player a small amount, as we don't
+     * want this constraint to apply in weird cases like hitting a corner.
+     * Otherwise, this constraint might throw the player large distances.
+     */
+    if (((bp - body_.position).squared()
+           + (lp - leftWheel_.position).squared()
+           + (rp - rightWheel_.position).squared()) > 1E-3) {
+      return;
+    }
+
+    body_.position = bp;
+    leftWheel_.position = lp;
+    rightWheel_.position = rp;
+
+    /*
+     * Include a small fudge factor to push the player away from the wall, so
+     * that the player can turn. Otherwise, the player seemingly gets locked to
+     * the wall surface!
+     */
+    body_.position += j.normal * 1E-4;
+    leftWheel_.position += j.normal * 1E-4;
+    rightWheel_.position += j.normal * 1E-4;
+    body_.previousPosition += j.normal * 1E-4;
+    leftWheel_.previousPosition += j.normal * 1E-4;
+    rightWheel_.previousPosition += j.normal * 1E-4;
+  }
+}
+
+void Player::constrainBody() {
+  Constraints::distance(leftWheel_, rightWheel_, axleLength);
+  Constraints::distance(leftWheel_, body_, axleLength / 2.f);
+  Constraints::distance(rightWheel_, body_, axleLength / 2.f);
+
+  static const float d
+      = sqrtf(axleLength * axleLength / 4.f
+          + counterWeightOffset * counterWeightOffset);
+  Constraints::distance(leftWheel_, counterWeight_, d);
+  Constraints::distance(rightWheel_, counterWeight_, d);
+  Constraints::distance(body_, counterWeight_, counterWeightOffset);
+}
+
+void Player::constrainInternal() {
+  leftWheel_.constrainDirection(z_);
+  rightWheel_.constrainDirection(z_);
+  constrainBody();
+
+  origin_ = body_.position;
+  z_ = (rightWheel_.position - leftWheel_.position) / axleLength;
+  x_ = (body_.position - counterWeight_.position) / counterWeightOffset;
+  y_ = -x_.cross(z_);
+}
+
+void Player::Wheel::constrainDirection(const Vector& z) {
+  if (!friction()) {
+    angle += angleStep;
+    return; // no normal force, no friction
+  }
+
+  /* Don't let upwards-moving platforms launch us into the air. */
+  Vector cv = contactVelocity;
+  if (cv.y > 0) {
+    cv.y = 0.f;
+  }
+
+  /* Cancel the sideways component of the velocity. */
+  Vector v = position - previousPosition - cv;
+  v -= z * z.dot(v);
+  previousPosition = position - v - cv;
+
+  /* Rotate the wheel. */
+  Vector forward = z.cross(contactNormal);
+  Vector vForward = forward * forward.dot(v);
+  float l = vForward.length();
+  if (vForward.cross(contactNormal).dot(z) > 0) {
+    l *= -1;
+  }
+  angleStep = 360.f * l / (M_PI * 2 * wheelRadius);
+  angle += angleStep;
+}
+
+void Player::jiggle() {
+  Vector v = Vector::randomVector(ParticleSimulator::timeStep());
+  v.y = ParticleSimulator::timeStep();
+  leftWheel_.position += v;
+  rightWheel_.position += v;
+  body_.position += v;
+  counterWeight_.position += v;
+}
+
+void Player::move(Direction d) {
+  switch (d) {
+    case LEFT: {
+      if (turnState_ == NONE) {
+        turnState_ = LEFT;
+      } else {
+        turnState_ = NONE;
+      }
+      break;
+    }
+    case RIGHT: {
+      if (turnState_ == NONE) {
+        turnState_ = RIGHT;
+      } else {
+        turnState_ = NONE;
+      }
+      break;
+    }
+    case FORWARD: {
+      if (moveState_ == NONE) {
+        moveState_ = FORWARD;
+      } else {
+        moveState_ = NONE;
+      }
+      break;
+    }
+    case BACKWARD: {
+      if (moveState_ == NONE) {
+        moveState_ = BACKWARD;
+      } else {
+        moveState_ = NONE;
+      }
+      break;
+    }
+  }
+}
+
+void Player::stop() {
+  turnState_ = NONE;
+  moveState_ = NONE;
+}
+
+void Player::stop(Direction d) {
+  switch (d) {
+    case LEFT: {
+      if (turnState_ != NONE) {
+        turnState_ = NONE;
+      } else {
+        turnState_ = RIGHT;
+      }
+      break;
+    }
+    case RIGHT: {
+      if (turnState_ != NONE) {
+        turnState_ = NONE;
+      } else {
+        turnState_ = LEFT;
+      }
+      break;
+    }
+    case FORWARD: {
+      if (moveState_ != NONE) {
+        moveState_ = NONE;
+      } else {
+        moveState_ = BACKWARD;
+      }
+      break;
+    }
+    case BACKWARD: {
+      if (moveState_ != NONE) {
+        moveState_ = NONE;
+      } else {
+        moveState_ = FORWARD;
+      }
+      break;
+    }
+  }
+}

src/player.h

@@ -0,0 +1,117 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_PLAYER_H
+#define MBOSTOCK_PLAYER_H
+
+#include <vector>
+
+#include "model.h"
+#include "physics/particle.h"
+#include "physics/shape.h"
+#include "physics/vector.h"
+
+namespace mbostock {
+
+  class Player;
+  class RoomObject;
+  class UnaryForce;
+
+  class PlayerModel : public Model {
+  public:
+    PlayerModel(const Player& player);
+    virtual ~PlayerModel();
+
+    virtual void initialize();
+    virtual void display();
+
+  private:
+    float* orientation();
+    void displayAxes();
+
+    const Player& player_;
+    GLUquadric* quadric_;
+    float orientation_[16];
+    Model* wheelModel_;
+    Model* bodyModel_;
+
+    friend class PlayerWheelModel;
+    friend class PlayerBodyModel;
+  };
+
+  class Player {
+  public:
+    Player();
+
+    enum Direction { NONE, LEFT, RIGHT, FORWARD, BACKWARD };
+
+    void move(Direction d);
+    void stop(Direction d);
+    void stop();
+    void jiggle();
+    void setOrigin(const Vector& origin);
+    void setVelocity(const Vector& velocity);
+
+    void resetForces();
+    void applyForce(UnaryForce& force);
+    void step(const ParticleSimulator& s);
+    bool intersects(const Shape& s);
+    bool constrainOutside(const RoomObject& o);
+    void constrainInternal();
+
+    float mass() const;
+
+    inline const Vector& x() const { return x_; }
+    inline const Vector& y() const { return y_; }
+    inline const Vector& z() const { return z_; }
+    inline const Vector& origin() const { return origin_; }
+    inline float leftWheelAngle() const { return leftWheel_.angle; }
+    inline float rightWheelAngle() const { return rightWheel_.angle; }
+
+    bool leftWheelFriction() const { return leftWheel_.friction(); }
+    bool rightWheelFriction() const { return rightWheel_.friction(); }
+
+    inline Model& model() { return model_; }
+
+  private:
+    class Wheel : public Particle {
+    public:
+      Wheel();
+
+      inline bool friction() const { return contact; }
+      void applyContact(const RoomObject& o);
+      void applyForce(const Vector& z, float f);
+      void applyLinearDrag(float kd);
+      void applyQuadraticDrag(float kd);
+      void constrainDirection(const Vector& z);
+
+      bool contact;
+      Vector contactNormal;
+      Vector contactVelocity;
+
+      float angle;
+      float angleStep;
+    };
+
+    void constrainGlancing(const RoomObject& o);
+    void constrainBody();
+
+    Wheel leftWheel_;
+    Wheel rightWheel_;
+    Particle body_;
+    Particle counterWeight_;
+
+    Direction turnState_;
+    Direction moveState_;
+
+    Sphere sphere_;
+    Vector origin_;
+    Vector x_;
+    Vector y_;
+    Vector z_;
+
+    PlayerModel model_;
+  };
+
+}
+
+#endif

src/portal.cpp

@@ -0,0 +1,11 @@
+// -*- C++ -*-
+
+#include "physics/vector.h"
+#include "portal.h"
+
+using namespace mbostock;
+
+Portal::Portal(const Vector& min, const Vector& max,
+               int room, int origin, bool reset)
+    : box_(min, max), room_(room), origin_(origin), reset_(reset) {
+}

src/portal.h

@@ -0,0 +1,33 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_PORTAL_H
+#define MBOSTOCK_PORTAL_H
+
+#include "physics/shape.h"
+
+namespace mbostock {
+
+  class Vector;
+
+  class Portal {
+  public:
+    Portal(const Vector& min, const Vector& max,
+           int room, int origin, bool reset);
+
+    inline bool contains(const Vector& p) const { return box_.contains(p); }
+    inline const AxisAlignedBox& bounds() const { return box_; }
+    inline int room() const { return room_; }
+    inline int origin() const { return origin_; }
+    inline bool reset() const { return reset_; }
+
+  private:
+    AxisAlignedBox box_;
+    int room_;
+    int origin_;
+    bool reset_;
+  };
+
+}
+
+
+#endif

src/ramp.cpp

@@ -0,0 +1,35 @@
+// -*- C++ -*-
+
+#include "material.h"
+#include "model.h"
+#include "physics/shape.h"
+#include "physics/vector.h"
+#include "ramp.h"
+#include "room.h"
+
+using namespace mbostock;
+
+Ramp::Ramp(const Vector& x0, const Vector& x1,
+           const Vector& x2, const Vector& x3)
+    : wedge_(x0, x1, x2, x3), model_(wedge_) {
+}
+
+Model& Ramp::model() {
+  return model_;
+}
+
+const Shape& Ramp::shape() const {
+  return wedge_;
+}
+
+void Ramp::setMaterial(const Material& m) {
+  model_.setMaterial(m);
+}
+
+void Ramp::setTopMaterial(const Material& m) {
+  model_.setTopMaterial(m);
+}
+
+float Ramp::slip() const {
+  return model_.material().slip();
+}

src/ramp.h

@@ -0,0 +1,32 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_RAMP_H
+#define MBOSTOCK_RAMP_H
+
+#include "model.h"
+#include "physics/shape.h"
+#include "room_object.h"
+
+namespace mbostock {
+
+  class Ramp : public RoomObject {
+  public:
+    Ramp(const Vector& x0, const Vector& x1,
+         const Vector& x2, const Vector& x3);
+
+    virtual Model& model();
+    virtual const Shape& shape() const;
+    virtual float slip() const;
+
+    void setMaterial(const Material& m);
+    void setTopMaterial(const Material& m);
+
+  private:
+    const Wedge wedge_;
+    WedgeModel model_;
+  };
+
+}
+
+
+#endif

src/resource.cpp

@@ -0,0 +1,27 @@
+// -*- C++ -*-
+
+#include <fstream>
+#include <ios>
+#include <iostream>
+
+#include "resource.h"
+
+using namespace mbostock;
+
+const char* Resources::path() {
+  return "Contents/Resources/";
+}
+
+const char* Resources::readFile(const char* p) {
+  std::string fullPath(path());
+  fullPath.append(p);
+  std::ifstream file(fullPath.c_str());
+  file.seekg(0, std::ios::end);
+  std::ifstream::pos_type size = file.tellg();
+  file.seekg(0, std::ios::beg);
+  char* buffer = new char[1 + size];
+  file.read(buffer, size);
+  buffer[size] = '\0';
+  file.close();
+  return buffer;
+}

src/resource.h

@@ -0,0 +1,19 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_RESOURCE_H
+#define MBOSTOCK_RESOURCE_H
+
+namespace mbostock {
+
+  class Resources {
+  public:
+    static const char* path();
+    static const char* readFile(const char* path);
+
+  private:
+    Resources();
+  };
+
+}
+
+#endif

src/room.cpp

@@ -0,0 +1,175 @@
+// -*- C++ -*-
+
+#include <algorithm>
+#include <iostream>
+
+#include "lighting.h"
+#include "physics/transform.h"
+#include "physics/vector.h"
+#include "portal.h"
+#include "room.h"
+#include "room_force.h"
+#include "room_object.h"
+#include "sound.h"
+#include "trail.h"
+#include "world.h"
+
+using namespace mbostock;
+
+class RoomStaticModel : public Model {
+public:
+  RoomStaticModel(const Room& room) : room_(room) {}
+
+  virtual void initialize() {
+    staticObjects_.clear();
+    std::vector<RoomObject*>::const_iterator i;
+    for (i = room_.objects().begin(); i != room_.objects().end(); i++) {
+      RoomObject* o = *i;
+      if (!o->dynamic()) {
+        o->model().initialize();
+        staticObjects_.push_back(o);
+      }
+    }
+  }
+
+  virtual void display() {
+    std::vector<RoomObject*>::const_iterator i;
+    for (i = staticObjects_.begin(); i != staticObjects_.end(); i++) {
+      (*i)->model().display();
+    }
+  }
+
+private:
+  const Room& room_;
+  std::vector<RoomObject*> staticObjects_;
+};
+
+RoomModel::RoomModel(const Room& room)
+    : room_(room), staticModel_(Models::compile(new RoomStaticModel(room))) {
+}
+
+RoomModel::~RoomModel() {
+  delete staticModel_;
+}
+
+void RoomModel::initialize() {
+  room_.lighting().initialize();
+  staticModel_->initialize();
+  dynamicObjects_.clear();
+  std::vector<RoomObject*>::const_iterator i;
+  for (i = room_.objects().begin(); i != room_.objects().end(); i++) {
+    RoomObject* o = *i;
+    if (o->dynamic()) {
+      o->model().initialize();
+      dynamicObjects_.push_back(o);
+    }
+  }
+}
+
+void RoomModel::display() {
+  if (World::world()->paused()) {
+    World::world()->pauseLighting().display();
+  } else {
+    room_.lighting().display();
+  }
+  staticModel_->display();
+  std::vector<RoomObject*>::const_iterator i;
+  for (i = dynamicObjects_.begin(); i != dynamicObjects_.end(); i++) {
+    (*i)->model().display();
+  }
+}
+
+Room::Room()
+    : lighting_(&(Lightings::standard())), music_(NULL),
+      cameraBounds_(-Vector::INF(), Vector::INF()),
+      model_(*this), trail_(NULL) {
+}
+
+Room::~Room() {
+  std::vector<RoomObject*>::const_iterator io;
+  for (io = objects_.begin(); io != objects_.end(); io++) {
+    delete (*io);
+  }
+  std::vector<RoomOrigin*>::const_iterator ir;
+  for (ir = origins_.begin(); ir != origins_.end(); ir++) {
+    delete (*ir);
+  }
+  std::vector<Portal*>::const_iterator ip;
+  for (ip = portals_.begin(); ip != portals_.end(); ip++) {
+    delete (*ip);
+  }
+  std::vector<RoomForce*>::const_iterator ii;
+  for (ii = forces_.begin(); ii != forces_.end(); ii++) {
+    delete (*ii);
+  }
+  std::vector<Transform*>::const_iterator iz;
+  for (iz = transforms_.begin(); iz != transforms_.end(); iz++) {
+    delete (*iz);
+  }
+  std::vector<Trail*>::const_iterator it;
+  for (it = trails_.begin(); it != trails_.end(); it++) {
+    delete (*it);
+  }
+  if (trail_ != NULL) {
+    delete trail_;
+  }
+}
+
+void Room::nextTrail(const Vector& origin) {
+  if (trail_ != NULL) {
+    trails_.push_back(trail_);
+  }
+  trail_ = new Trail(origin);
+}
+
+void Room::setMusic(const Sound& music) {
+  music_ = &music;
+}
+
+void Room::setLighting(const Lighting& lighting) {
+  lighting_ = &lighting;
+}
+
+void Room::setCameraBounds(const Vector& min, const Vector& max) {
+  cameraBounds_ = AxisAlignedBox(min, max);
+}
+
+void Room::resetForces() {
+  std::vector<RoomObject*>::const_iterator i;
+  for (i = objects_.begin(); i != objects_.end(); i++) {
+    (*i)->resetForces();
+  }
+}
+
+void Room::applyForce(UnaryForce& force) {
+  std::vector<RoomObject*>::const_iterator i;
+  for (i = objects_.begin(); i != objects_.end(); i++) {
+    (*i)->applyForce(force);
+  }
+}
+
+void Room::step(const ParticleSimulator& s) {
+  std::vector<Transform*>::const_iterator ir;
+  for (ir = transforms_.begin(); ir != transforms_.end(); ir++) {
+    (*ir)->step();
+  }
+  std::vector<RoomObject*>::const_iterator i;
+  for (i = objects_.begin(); i != objects_.end(); i++) {
+    (*i)->step(s);
+  }
+}
+
+void Room::reset() {
+  std::vector<RoomObject*>::const_iterator i;
+  for (i = objects_.begin(); i != objects_.end(); i++) {
+    (*i)->reset();
+  }
+  std::vector<Transform*>::const_iterator iz;
+  for (iz = transforms_.begin(); iz != transforms_.end(); iz++) {
+    (*iz)->reset();
+  }
+}
+
+RoomOrigin::RoomOrigin(const Vector& position, const Vector& velocity)
+    : position_(position), velocity_(velocity) {
+}

src/room.h

@@ -0,0 +1,99 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_ROOM_H
+#define MBOSTOCK_ROOM_H
+
+#include <vector>
+
+#include "model.h"
+
+namespace mbostock {
+
+  class Lighting;
+  class ParticleSimulator;
+  class Portal;
+  class Room;
+  class RoomForce;
+  class RoomObject;
+  class RoomOrigin;
+  class Shape;
+  class Sound;
+  class Trail;
+  class Transform;
+  class UnaryForce;
+
+  class RoomModel : public Model {
+  public:
+    RoomModel(const Room& room);
+    virtual ~RoomModel();
+
+    virtual void initialize();
+    virtual void display();
+
+  private:
+    const Room& room_;
+    std::vector<RoomObject*> dynamicObjects_;
+    Model* staticModel_;
+  };
+
+  class Room {
+  public:
+    Room();
+    ~Room();
+
+    inline void addForce(RoomForce* o) { forces_.push_back(o); }
+    inline void addOrigin(RoomOrigin* o) { origins_.push_back(o); }
+    inline void addObject(RoomObject* o) { objects_.push_back(o); }
+    inline void addPortal(Portal* p) { portals_.push_back(p); }
+    inline void addTransform(Transform* r) { transforms_.push_back(r); }
+
+    inline const std::vector<RoomForce*>& forces() const { return forces_; }
+    inline const std::vector<RoomOrigin*>& origins() const { return origins_; }
+    inline const std::vector<RoomObject*>& objects() const { return objects_; }
+    inline const std::vector<Portal*>& portals() const { return portals_; }
+    inline const std::vector<Trail*>& trails() const { return trails_; }
+
+    inline Model& model() { return model_; }
+    inline Trail& trail() { return *trail_; }
+    inline const Lighting& lighting() const { return *lighting_; }
+    inline const AxisAlignedBox& cameraBounds() const { return cameraBounds_; }
+    inline const Sound* music() const { return music_; }
+
+    void setMusic(const Sound& music);
+    void setLighting(const Lighting& lighting);
+    void setCameraBounds(const Vector& min, const Vector& max);
+    void resetForces();
+    void applyForce(UnaryForce& force);
+    void step(const ParticleSimulator& s);
+    void reset();
+    void nextTrail(const Vector& origin);
+
+  private:
+    std::vector<RoomForce*> forces_;
+    std::vector<RoomOrigin*> origins_;
+    std::vector<RoomObject*> objects_;
+    std::vector<Portal*> portals_;
+    std::vector<Trail*> trails_;
+    std::vector<Transform*> transforms_;
+    const Lighting* lighting_;
+    const Sound* music_;
+    AxisAlignedBox cameraBounds_;
+    RoomModel model_;
+    Trail* trail_;
+  };
+
+  class RoomOrigin {
+  public:
+    RoomOrigin(const Vector& position, const Vector& velocity);
+
+    inline const Vector& position() const { return position_; }
+    inline const Vector& velocity() const { return velocity_; }
+
+  private:
+    Vector position_;
+    Vector velocity_;
+  };
+
+}
+
+#endif

src/room_force.cpp

@@ -0,0 +1,20 @@
+// -*- C++ -*-
+
+#include "room_force.h"
+#include "physics/particle.h"
+#include "physics/vector.h"
+
+using namespace mbostock;
+
+ConstantRoomForce::ConstantRoomForce(const Vector& min, const Vector& max,
+                                     const Vector& f)
+    : box_(min, max), force_(f) {
+}
+
+const Shape& ConstantRoomForce::shape() const {
+  return box_;
+}
+
+Vector ConstantRoomForce::force(const Particle& p) {
+  return box_.contains(p.position) ? force_ : Vector::ZERO();
+}

src/room_force.h

@@ -0,0 +1,45 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_ROOM_FORCE_H
+#define MBOSTOCK_ROOM_FORCE_H
+
+#include "physics/force.h"
+#include "physics/shape.h"
+#include "physics/vector.h"
+
+namespace mbostock {
+
+  class Particle;
+  class Shape;
+
+  /**
+   * Defines a force that applies to the player when the player intersects a
+   * given shape. For example, this can be used to represent a fan blowing on
+   * the player, which the force of the fan stronger near the fan blades,
+   * falling off to zero outside of the fan's stream.
+   */
+  class RoomForce : public UnaryForce {
+  public:
+    virtual const Shape& shape() const = 0;
+  };
+
+  /**
+   * A simple room force that applies a constant force within an axis-aligned
+   * box. More elaborate forces can be approximated by composing multiple (even
+   * overlapping) constant room forces.
+   */
+  class ConstantRoomForce : public RoomForce {
+  public:
+    ConstantRoomForce(const Vector& min, const Vector& max, const Vector& f);
+
+    virtual const Shape& shape() const;
+    virtual Vector force(const Particle& p);
+
+  private:
+    AxisAlignedBox box_;
+    Vector force_;
+  };
+
+}
+
+#endif

src/room_object.cpp

@@ -0,0 +1,40 @@
+// -*- C++ -*-
+
+#include "physics/vector.h"
+#include "room_object.h"
+
+using namespace mbostock;
+
+bool RoomObject::dynamic() const {
+  return false;
+}
+
+Vector RoomObject::velocity(const Vector& x) const {
+  return Vector::ZERO();
+}
+
+float RoomObject::slip() const {
+  return 0.f;
+}
+
+void RoomObject::resetForces() {
+}
+
+void RoomObject::applyForce(UnaryForce& force) {
+}
+
+void RoomObject::applyWeight(float w, const Vector& x) {
+}
+
+void RoomObject::step(const ParticleSimulator& s) {
+}
+
+void RoomObject::constrainInternal() {
+}
+
+void RoomObject::reset() {
+}
+
+bool DynamicRoomObject::dynamic() const {
+  return true;
+}

src/room_object.h

@@ -0,0 +1,36 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_ROOM_OBJECT_H
+#define MBOSTOCK_ROOM_OBJECT_H
+
+namespace mbostock {
+
+  class Model;
+  class ParticleSimulator;
+  class Shape;
+  class UnaryForce;
+  class Vector;
+
+  class RoomObject {
+  public:
+    virtual Model& model() = 0;
+    virtual const Shape& shape() const = 0;
+    virtual bool dynamic() const;
+    virtual Vector velocity(const Vector& x) const;
+    virtual float slip() const;
+    virtual void resetForces();
+    virtual void applyForce(UnaryForce& force);
+    virtual void applyWeight(float w, const Vector& x);
+    virtual void step(const ParticleSimulator& s);
+    virtual void constrainInternal();
+    virtual void reset();
+  };
+
+  class DynamicRoomObject : public RoomObject {
+  public:
+    virtual bool dynamic() const;
+  };
+
+}
+
+#endif

src/rotating.cpp

@@ -0,0 +1,46 @@
+// -*- C++ -*-
+
+#include <math.h>
+
+#include "physics/particle.h"
+#include "rotating.h"
+
+using namespace mbostock;
+
+RotatingModel::RotatingModel(Model& m, const Rotation& r)
+    : model_(m), rotation_(r) {
+}
+
+void RotatingModel::initialize() {
+  model_.initialize();
+}
+
+void RotatingModel::display() {
+  glPushMatrix();
+  glTranslatev(rotation_.origin());
+  glRotatev(rotation_.angle(), rotation_.axis());
+  glTranslatev(-rotation_.origin());
+  model_.display();
+  glPopMatrix();
+}
+
+RotatingRoomObject::RotatingRoomObject(RoomObject* o, const Rotation& r)
+    : TransformingRoomObject(o), rotation_(r),
+      shape_(o->shape(), r), model_(o->model(), r) {
+}
+
+Model& RotatingRoomObject::model() {
+  return model_;
+}
+
+const Shape& RotatingRoomObject::shape() const {
+  return shape_;
+}
+
+Vector RotatingRoomObject::velocity(const Vector& x) const {
+  Vector v = rotation_.velocity(x);
+  if (object_->dynamic()) {
+    v += rotation_.vector(object_->velocity(rotation_.pointInverse(x)));
+  }
+  return v;
+}

src/rotating.h

@@ -0,0 +1,57 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_ROTATING_H
+#define MBOSTOCK_ROTATING_H
+
+#include "model.h"
+#include "physics/rotation.h"
+#include "physics/vector.h"
+#include "room_force.h"
+#include "room_object.h"
+#include "transforming.h"
+
+namespace mbostock {
+
+  class RotatingModel : public Model {
+  public:
+    RotatingModel(Model& m, const Rotation& r);
+
+    virtual void initialize();
+    virtual void display();
+
+  private:
+    Model& model_;
+    const Rotation& rotation_;
+  };
+
+  class RotatingRoomObject : public TransformingRoomObject {
+  public:
+    RotatingRoomObject(RoomObject* o, const Rotation& r);
+
+    virtual Model& model();
+    virtual const Shape& shape() const;
+    virtual Vector velocity(const Vector& x) const;
+
+  private:
+    const Rotation& rotation_;
+    RotatingShape shape_;
+    RotatingModel model_;
+  };
+
+  class RotatingRoomForce : public RoomForce {
+  public:
+    RotatingRoomForce(RoomForce* f, const Rotation& r);
+    virtual ~RotatingRoomForce();
+
+    virtual const Shape& shape() const;
+    virtual Vector force(const Particle& p);
+
+  private:
+    RoomForce* force_;
+    const Rotation& rotation_;
+    RotatingShape shape_;
+  };
+
+}
+
+#endif

src/seesaw.cpp

@@ -0,0 +1,117 @@
+// -*- C++ -*-
+
+#include <OpenGL/gl.h>
+
+#include "material.h"
+#include "physics/constraint.h"
+#include "physics/force.h"
+#include "physics/particle.h"
+#include "physics/shape.h"
+#include "physics/vector.h"
+#include "seesaw.h"
+
+using namespace mbostock;
+
+static const float centerOffset = .1f;
+
+Seesaw::Seesaw(const Vector& min, const Vector& max, float mass)
+    : origin_((min + max) / 2.f), size_(max - min), drag_(1.f),
+      model_(box_) {
+  left_.inverseMass = 1.f / (mass * .1f);
+  right_.inverseMass = 1.f / (mass * .1f);
+  center_.inverseMass = 1.f / (mass * .8f);
+  reset();
+}
+
+const Shape& Seesaw::shape() const {
+  return box_;
+}
+
+Model& Seesaw::model() {
+  return model_;
+}
+
+void Seesaw::resetForces() {
+  left_.force = Vector::ZERO();
+  right_.force = Vector::ZERO();
+  center_.force = Vector::ZERO();
+  drag_.apply(left_);
+  drag_.apply(right_);
+  drag_.apply(center_);
+}
+
+void Seesaw::applyForce(UnaryForce& force) {
+  force.apply(left_);
+  force.apply(right_);
+  force.apply(center_);
+}
+
+void Seesaw::applyWeight(float w, const Vector& x) {
+  if (x.x < origin_.x) {
+    left_.force.y -= w * (2 * (origin_.x - x.x) / size_.x);
+  } else {
+    right_.force.y += w * (2 * (origin_.x - x.x) / size_.x);
+  }
+}
+
+void Seesaw::step(const ParticleSimulator& s) {
+  s.step(left_);
+  s.step(right_);
+  s.step(center_);
+  updateBox();
+}
+
+void Seesaw::updateBox() {
+  Vector x = (right_.position - left_.position) / size_.x;
+  Vector y = -x.cross(Vector::Z()) * size_.y / 2.f;
+  Vector z = Vector::Z() * size_.z / 2.f;
+  box_ = Box(
+      left_.position - z - y,
+      right_.position - z - y,
+      right_.position + z - y,
+      left_.position + z - y,
+      right_.position - z + y,
+      left_.position - z + y,
+      left_.position + z + y,
+      right_.position + z + y);
+}
+
+void Seesaw::constrainInternal() {
+  Vector v = origin_ - (left_.position + right_.position) / 2.f;
+  right_.position += v;
+  left_.position += v;
+
+  Constraints::distance(right_, origin_, size_.x / 2.f);
+  Constraints::distance(left_, origin_, size_.x / 2.f);
+  Constraints::distance(right_, left_, size_.x);
+  Constraints::distance(center_, origin_, centerOffset);
+
+  float d = sqrtf(size_.x * size_.x / 4.f + centerOffset * centerOffset);
+  Constraints::distance(left_, center_, d);
+  Constraints::distance(right_, center_, d);
+}
+
+void Seesaw::reset() {
+  left_.position = origin_;
+  left_.position.x -= size_.x / 2.f;
+  right_.position = origin_;
+  right_.position.x += size_.x / 2.f;
+  center_.position = origin_;
+  center_.position.y -= centerOffset;
+  left_.previousPosition = left_.position;
+  right_.previousPosition = right_.position;
+  center_.previousPosition = center_.position;
+  updateBox();
+}
+
+void Seesaw::setMaterial(const Material& m) {
+  model_.setMaterial(m);
+}
+
+void Seesaw::setTopMaterial(const Material& m) {
+  model_.setTopMaterial(m);
+}
+
+float Seesaw::slip() const {
+  return model_.material().slip();
+}

src/seesaw.h

@@ -0,0 +1,51 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_SEESAW_H
+#define MBOSTOCK_SEESAW_H
+
+#include "model.h"
+#include "physics/force.h"
+#include "physics/particle.h"
+#include "physics/shape.h"
+#include "physics/vector.h"
+#include "room_object.h"
+
+namespace mbostock {
+
+  class Material;
+
+  class Seesaw : public DynamicRoomObject {
+  public:
+    Seesaw(const Vector& min, const Vector& max, float mass);
+
+    virtual Model& model();
+    virtual const Shape& shape() const;
+    virtual void resetForces();
+    virtual void applyForce(UnaryForce& force);
+    virtual void applyWeight(float w, const Vector& x);
+    virtual void step(const ParticleSimulator& s);
+    virtual void constrainInternal();
+    virtual void reset();
+    virtual float slip() const;
+
+    void setMaterial(const Material& m);
+    void setTopMaterial(const Material& m);
+
+  private:
+    void updateBox();
+
+    const Vector origin_;
+    const Vector size_;
+    LinearDragForce drag_;
+
+    Box box_;
+    Particle left_;
+    Particle right_;
+    Particle center_;
+
+    BoxModel model_;
+  };
+
+}
+
+#endif

src/shader.cpp

@@ -0,0 +1,92 @@
+// -*- C++ -*-
+
+#include <stdlib.h>
+
+#include "model.h"
+#include "resource.h"
+#include "shader.h"
+
+using namespace mbostock;
+
+static Shader* defaultShader_ = NULL;
+static Shader* normalShader_ = NULL;
+static Shader* wireframeShader_ = NULL;
+
+class DefaultShader : public Shader {
+public:
+  virtual void display(Model& model) {
+    model.display();
+  }
+};
+
+class WireframeShader : public Shader {
+public:
+  virtual void display(Model& model) {
+    glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
+    model.display();
+    glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
+  }
+};
+
+Shader* Shaders::defaultShader() {
+  if (defaultShader_ == NULL) {
+    defaultShader_ = new DefaultShader();
+  }
+  return defaultShader_;
+}
+
+Shader* Shaders::wireframeShader() {
+  if (wireframeShader_ == NULL) {
+    wireframeShader_ = new WireframeShader();
+  }
+  return wireframeShader_;
+}
+
+Shader* Shaders::normalShader() {
+  if (normalShader_ == NULL) {
+    normalShader_ = new GlslShader("normal.vert", "normal.frag");
+  }
+  return normalShader_;
+}
+
+GlslShader::GlslShader(const char* vertexPath, const char* fragmentPath)
+  : vertexPath_(vertexPath), fragmentPath_(fragmentPath),
+    program_(GL_NONE) {
+}
+
+GlslShader::~GlslShader() {
+  if (program_ != GL_NONE) {
+    glDeleteProgram(program_);
+  }
+}
+
+void GlslShader::initialize() {
+  if (program_ != GL_NONE) {
+    glDeleteProgram(program_);
+  }
+  program_ = glCreateProgram();
+  attach(vertexPath_, GL_VERTEX_SHADER);
+  attach(fragmentPath_, GL_FRAGMENT_SHADER);
+  link();
+}
+
+void GlslShader::attach(const char* path, GLenum shaderType) {
+  GLuint shader = glCreateShader(shaderType);
+  const char* source = Resources::readFile(path);
+  glShaderSource(shader, 1, &source, NULL);
+  delete[] source;
+  glCompileShader(shader);
+  glAttachShader(program_, shader);
+}
+
+void GlslShader::link() {
+  glLinkProgram(program_);
+}
+
+void GlslShader::display(Model& model) {
+  glDisable(GL_LIGHTING);
+  glUseProgram(program_);
+  model.display();
+  glUseProgram(GL_NONE);
+  glEnable(GL_LIGHTING);
+}

src/shader.h

@@ -0,0 +1,49 @@
+// -*- C++ -*-
+
+#ifndef _SHADER_H
+#define _SHADER_H
+
+#include <OpenGL/gl.h>
+
+namespace mbostock {
+
+  class Model;
+
+  class Shader {
+  public:
+    virtual ~Shader() {};
+
+    virtual void initialize() {}
+    virtual void display(Model& model) = 0;
+  };
+
+  class GlslShader : public Shader {
+  public:
+    GlslShader(const char* vertexPath, const char* fragmentPath);
+    virtual ~GlslShader();
+
+    virtual void initialize();
+    virtual void display(Model& model);
+
+  private:
+    void attach(const char* path, GLenum shaderType);
+    void link();
+
+    const char* vertexPath_;
+    const char* fragmentPath_;
+    GLuint program_;
+  };
+
+  class Shaders {
+  public:
+    static Shader* defaultShader();
+    static Shader* normalShader();
+    static Shader* wireframeShader();
+
+  private:
+    Shaders();
+  };
+
+}
+
+#endif

src/simulation.cpp

@@ -0,0 +1,45 @@
+// -*- C++ -*-
+
+#include <SDL/sdl.h>
+
+#include "simulation.h"
+
+using namespace mbostock;
+
+/*
+ * If for some reason we haven't been able to run the simulation in a long time
+ * (e.g., the computer went to sleep for two hours!), then we don't want to run
+ * the simulation a zillion times to catch up. Instead we just drop some frames
+ * and hope it gets better.
+ */
+static const uint32_t maxSkippedMs = 500;
+
+Simulation::Simulation(uint32_t timeStepMs)
+  : timeStepMs_(timeStepMs), skippedMs_(0), lastTimeMs_(0), paused_(false) {
+}
+
+uint32_t Simulation::elapsedMillis() {
+  uint32_t currentTimeMs = SDL_GetTicks();
+  uint32_t differenceMs = currentTimeMs - lastTimeMs_;
+  lastTimeMs_ = currentTimeMs;
+  return differenceMs;
+}
+
+void Simulation::togglePaused() {
+  paused_ = !paused_;
+}
+
+void Simulation::simulate() {
+  if (paused_ || (lastTimeMs_ == 0)) {
+    lastTimeMs_ = SDL_GetTicks();
+    return;
+  }
+  skippedMs_ += elapsedMillis();
+  if (skippedMs_ > maxSkippedMs) {
+    skippedMs_ = timeStepMs_;
+  }
+  while (skippedMs_ >= timeStepMs_) {
+    step();
+    skippedMs_ -= timeStepMs_;
+  }
+}

src/simulation.h

@@ -0,0 +1,34 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_SIMULATION_H
+#define MBOSTOCK_SIMULATION_H
+
+#include <stdint.h>
+
+namespace mbostock {
+
+  class Simulation {
+  public:
+    Simulation(uint32_t timeStepMs);
+    virtual ~Simulation() {}
+
+    virtual void togglePaused();
+    inline bool paused() const { return paused_; }
+
+    void simulate();
+
+  protected:
+    virtual void step() = 0;
+
+  private:
+    uint32_t elapsedMillis();
+
+    uint32_t timeStepMs_;
+    uint32_t skippedMs_;
+    uint32_t lastTimeMs_;
+    bool paused_;
+  };
+
+}
+
+#endif

src/sound.cpp

@@ -0,0 +1,165 @@
+// -*- C++ -*-
+
+#include <SDL/SDL_error.h>
+#include <SDL_mixer/SDL_mixer.h>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "resource.h"
+#include "sound.h"
+
+using namespace mbostock;
+
+namespace mbostock {
+
+  class SoundImpl : public Sound {
+  public:
+    SoundImpl(const char* path);
+
+    inline const char* path() const { return path_.c_str(); }
+
+  private:
+    std::string path_;
+  };
+
+  class MusicImpl : public SoundImpl {
+  public:
+    MusicImpl(const char* path);
+    virtual ~MusicImpl();
+
+    virtual void play(int loops) const;
+    virtual void stop() const;
+
+  private:
+    Mix_Music* music_;
+  };
+
+  class ChunkImpl : public SoundImpl {
+  public:
+    ChunkImpl(const char* path);
+    virtual ~ChunkImpl();
+
+    virtual void play(int loops) const;
+    virtual void stop() const;
+
+  private:
+    Mix_Chunk* chunk_;
+    mutable int channel_;
+  };
+
+}
+
+SoundImpl::SoundImpl(const char* path)
+    : path_(path) {
+}
+
+MusicImpl::MusicImpl(const char* path)
+    : SoundImpl(path) {
+  std::string fullPath(Resources::path());
+  fullPath.append(path);
+  music_ = Mix_LoadMUS(fullPath.c_str());
+  if (music_ == NULL) {
+    std::cerr << SDL_GetError() << "\n";
+  }
+}
+
+MusicImpl::~MusicImpl() {
+  if (music_ != NULL) {
+    Mix_FreeMusic(music_);
+  }
+}
+
+void MusicImpl::play(int loops) const {
+  stop();
+  if (music_ != NULL) {
+    Mix_PlayMusic(music_, loops);
+  }
+}
+
+void MusicImpl::stop() const {
+  Mix_HaltMusic();
+}
+
+ChunkImpl::ChunkImpl(const char* path)
+    : SoundImpl(path), channel_(-1) {
+  std::string fullPath(Resources::path());
+  fullPath.append(path);
+  chunk_ = Mix_LoadWAV(fullPath.c_str());
+  if (chunk_ == NULL) {
+    std::cerr << SDL_GetError() << "\n";
+  }
+}
+
+ChunkImpl::~ChunkImpl() {
+  if (chunk_ != NULL) {
+    Mix_FreeChunk(chunk_);
+  }
+}
+
+void ChunkImpl::play(int loops) const {
+  stop();
+  if (chunk_ != NULL) {
+    channel_ = Mix_PlayChannel(-1, chunk_, loops);
+  }
+}
+
+void ChunkImpl::stop() const {
+  if (channel_ != -1) {
+    Mix_HaltChannel(channel_);
+    channel_ = -1;
+  }
+}
+
+static std::vector<SoundImpl*>& sounds() {
+  static std::vector<SoundImpl*> v;
+  return v;
+}
+
+void Sounds::initialize() {
+  Mix_OpenAudio(44100, AUDIO_S16SYS, 2, 1024);
+}
+
+void Sounds::dispose() {
+  Mix_HaltMusic();
+  Mix_HaltChannel(-1);
+  std::vector<SoundImpl*>::const_iterator i;
+  for (i = sounds().begin(); i != sounds().end(); i++) {
+    delete *i;
+  }
+  Mix_CloseAudio();
+}
+
+Sound& Sounds::fromFile(const char* path) {
+  /* First check to see if we've loaded this sound already. */
+  std::vector<SoundImpl*>::const_iterator i;
+  for (i = sounds().begin(); i != sounds().end(); i++) {
+    SoundImpl& s = **i;
+    if (strcmp(s.path(), path) == 0) {
+      return s;
+    }
+  }
+
+  /* If not, load the new sound. */
+  std::string spath(path);
+  SoundImpl* sound =
+      ((spath.rfind(".mid", std::string::npos, 4) != -1)
+       || (spath.rfind(".ogg", std::string::npos, 4) != -1)
+       || (spath.rfind(".mp3", std::string::npos, 4) != -1))
+          ? (SoundImpl*) new MusicImpl(path)
+          : (SoundImpl*) new ChunkImpl(path);
+  sounds().push_back(sound);
+  return *sound;
+}
+
+/* It appears that Mix_PauseMusic is unsupported for MIDI. :\ */
+
+void Sounds::pause() {
+  Mix_PauseMusic();
+  Mix_Pause(-1);
+}
+
+void Sounds::resume() {
+  Mix_ResumeMusic();
+  Mix_Resume(-1);
+}

src/sound.h

@@ -0,0 +1,29 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_SOUND_H
+#define MBOSTOCK_SOUND_H
+
+namespace mbostock {
+
+  class Sound {
+  public:
+    virtual ~Sound() {}
+    virtual void play(int loops = 0) const = 0;
+    virtual void stop() const = 0;
+  };
+
+  class Sounds {
+  public:
+    static void initialize();
+    static void pause();
+    static void resume();
+    static void dispose();
+    static Sound& fromFile(const char* path);
+
+  private:
+    Sounds();
+  };
+
+}
+
+#endif

src/switch.cpp

@@ -0,0 +1,49 @@
+// -*- C++ -*-
+
+#include "material.h"
+#include "physics/transform.h"
+#include "physics/vector.h"
+#include "switch.h"
+
+using namespace mbostock;
+
+Switch::Switch(const Vector& min, const Vector& max)
+    : AxisAlignedBlock(min, max), activeMaterial_(NULL) {
+}
+
+bool Switch::dynamic() const {
+  return (activeMaterial_ != NULL);
+}
+
+void Switch::addTarget(Transform& t) {
+  t.enable(false);
+  targets_.push_back(&t);
+}
+
+void Switch::setActiveMaterial(const Material& m) {
+  activeMaterial_ = &m;
+  inactiveMaterial_ = &model_.material();
+  inactiveTopMaterial_ = &model_.topMaterial();
+}
+
+void Switch::reset() {
+  if (activeMaterial_ != NULL) {
+    setMaterial(*inactiveMaterial_);
+    setTopMaterial(*inactiveTopMaterial_);
+  }
+  std::vector<Transform*>::const_iterator i;
+  for (i = targets_.begin(); i != targets_.end(); i++) {
+    (*i)->enable(false);
+  }
+}
+
+void Switch::applyWeight(float w, const Vector& x) {
+  if (activeMaterial_ != NULL) {
+    setMaterial(*activeMaterial_);
+    setTopMaterial(*activeMaterial_);
+  }
+  std::vector<Transform*>::const_iterator i;
+  for (i = targets_.begin(); i != targets_.end(); i++) {
+    (*i)->enable();
+  }
+}

src/switch.h

@@ -0,0 +1,34 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_SWITCH_H
+#define MBOSTOCK_SWITCH_H
+
+#include <vector>
+
+#include "block.h"
+
+namespace mbostock {
+
+  class Transform;
+
+  class Switch : public AxisAlignedBlock {
+  public:
+    Switch(const Vector& min, const Vector& max);
+
+    virtual bool dynamic() const;
+    virtual void applyWeight(float w, const Vector& x);
+    virtual void reset();
+
+    void addTarget(Transform& t);
+    void setActiveMaterial(const Material& m);
+
+  private:
+    std::vector<Transform*> targets_;
+    const Material* inactiveMaterial_;
+    const Material* inactiveTopMaterial_;
+    const Material* activeMaterial_;
+  };
+
+}
+
+#endif

src/texture.cpp

@@ -0,0 +1,99 @@
+// -*- C++ -*-
+
+#include <OpenGL/gl.h>
+#include <OpenGL/glu.h>
+#include <SDL/sdl.h>
+#include <SDL_image/SDL_image.h>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "resource.h"
+#include "texture.h"
+
+using namespace mbostock;
+
+class TextureImpl : public Texture {
+public:
+  TextureImpl(const char* path);
+  virtual ~TextureImpl();
+
+  void initialize();
+  virtual void bind() const;
+
+private:
+  GLuint id_;
+  bool alpha_;
+  const std::string path_;
+};
+
+TextureImpl::TextureImpl(const char* path)
+    : id_(GL_NONE), alpha_(false), path_(path) {
+}
+
+TextureImpl::~TextureImpl() {
+  if (id_ != GL_NONE) {
+    glDeleteTextures(1, &id_);
+  }
+}
+
+void TextureImpl::initialize() {
+  if (id_ != GL_NONE) {
+    glDeleteTextures(1, &id_);
+  }
+
+  std::string path(Resources::path());
+  path.append(path_);
+  SDL_Surface* image = IMG_Load(path.c_str());
+  if (image == NULL) {
+    std::cerr << "Couldn't load " << path << ": " << SDL_GetError() << "\n";
+    id_ = GL_NONE;
+    return;
+  }
+
+	glGenTextures(1, &id_);
+	glBindTexture(GL_TEXTURE_2D, id_);
+  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_NEAREST);
+  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
+
+  alpha_ = (image->format->BytesPerPixel == 4);
+  GLenum byteOrder = alpha_
+          ? ((image->format->Rmask == 0x000000FF) ? GL_RGBA : GL_BGRA)
+          : ((image->format->Rmask == 0x000000FF) ? GL_RGB : GL_BGR);
+
+	glTexImage2D(GL_TEXTURE_2D, 0,
+      image->format->BytesPerPixel, image->w, image->h, 0,
+      byteOrder, GL_UNSIGNED_BYTE, image->pixels);
+
+  gluBuild2DMipmaps(GL_TEXTURE_2D,
+      image->format->BytesPerPixel, image->w, image->h,
+      byteOrder, GL_UNSIGNED_BYTE, image->pixels);
+
+  SDL_FreeSurface(image);
+}
+
+void TextureImpl::bind() const {
+  if (alpha_) {
+    glEnable(GL_BLEND);
+    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
+  }
+  glBindTexture(GL_TEXTURE_2D, id_);
+}
+
+static std::vector<TextureImpl*>& textures() {
+  static std::vector<TextureImpl*> v;
+  return v;
+}
+
+const Texture& Textures::fromFile(const char* path) {
+  TextureImpl* texture = new TextureImpl(path);
+  textures().push_back(texture);
+  return *texture;
+}
+
+void Textures::initialize() {
+  std::vector<TextureImpl*>::const_iterator i;
+  for (i = textures().begin(); i != textures().end(); i++) {
+    (*i)->initialize();
+  }
+}

src/texture.h

@@ -0,0 +1,28 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_TEXTURE_H
+#define MBOSTOCK_TEXTURE_H
+
+#include <OpenGL/gl.h>
+
+namespace mbostock {
+
+  class Texture {
+  public:
+    virtual ~Texture() {}
+
+    virtual void bind() const = 0;
+  };
+
+  class Textures {
+  public:
+    static const Texture& fromFile(const char* path);
+    static void initialize();
+
+  private:
+    Textures();
+  };
+
+}
+
+#endif

src/trail.cpp

@@ -0,0 +1,38 @@
+// -*- C++ -*-
+
+#include "trail.h"
+
+using namespace mbostock;
+
+static const float minsq = .1f * .1f;
+
+Trail::Trail(const Vector& origin) {
+  points_.push_back(origin);
+}
+
+bool Trail::add(const Vector& p) {
+  if ((points_.back() - p).squared() >= minsq) {
+    points_.push_back(p);
+    return true;
+  }
+  return false;
+}
+
+TrailModel::TrailModel(const Trail& trail)
+    : trail_(trail) {
+}
+
+void TrailModel::display() {
+  glDisable(GL_LIGHTING);
+  glEnable(GL_BLEND);
+  glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
+  glColor4f(.6f, .2f, .3f, .5f);
+  glBegin(GL_LINE_STRIP);
+  std::vector<Vector>::const_iterator i;
+  for (i = trail_.points().begin(); i != trail_.points().end(); i++) {
+    glVertexv(*i);
+  }
+  glEnd();
+  glDisable(GL_BLEND);
+  glEnable(GL_LIGHTING);
+}

src/trail.h

@@ -0,0 +1,36 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_TRAIL_H
+#define MBOSTOCK_TRAIL_H
+
+#include <vector>
+
+#include "model.h"
+#include "physics/vector.h"
+
+namespace mbostock {
+
+  class Trail {
+  public:
+    Trail(const Vector& origin);
+
+    bool add(const Vector& p);
+    inline const std::vector<Vector>& points() const { return points_; }
+
+  private:
+    std::vector<Vector> points_;
+  };
+
+  class TrailModel : public Model {
+  public:
+    TrailModel(const Trail& trail);
+
+    virtual void display();
+
+  private:
+    const Trail& trail_;
+  };
+
+}
+
+#endif

src/transforming.cpp

@@ -0,0 +1,41 @@
+// -*- C++ -*-
+
+#include "transforming.h"
+
+using namespace mbostock;
+
+TransformingRoomObject::TransformingRoomObject(RoomObject* o)
+    : object_(o) {
+}
+
+TransformingRoomObject::~TransformingRoomObject() {
+  delete object_;
+}
+
+void TransformingRoomObject::step(const ParticleSimulator& s) {
+  object_->step(s);
+}
+
+float TransformingRoomObject::slip() const {
+  return object_->slip();
+}
+
+void TransformingRoomObject::resetForces() {
+  object_->resetForces();
+}
+
+void TransformingRoomObject::applyForce(UnaryForce& force) {
+  object_->applyForce(force); // TODO transform the force
+}
+
+void TransformingRoomObject::applyWeight(float w, const Vector& x) {
+  object_->applyWeight(w, x); // TODO transform the point
+}
+
+void TransformingRoomObject::constrainInternal() {
+  object_->constrainInternal();
+}
+
+void TransformingRoomObject::reset() {
+  object_->reset();
+}

src/transforming.h

@@ -0,0 +1,29 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_TRANSFORMING_H
+#define MBOSTOCK_TRANSFORMING_H
+
+#include "room_object.h"
+
+namespace mbostock {
+
+  class TransformingRoomObject : public DynamicRoomObject {
+  public:
+    TransformingRoomObject(RoomObject* o);
+    virtual ~TransformingRoomObject();
+
+    virtual float slip() const;
+    virtual void resetForces();
+    virtual void applyForce(UnaryForce& force);
+    virtual void applyWeight(float w, const Vector& x);
+    virtual void step(const ParticleSimulator& s);
+    virtual void constrainInternal();
+    virtual void reset();
+
+  protected:
+    RoomObject* object_;
+  };
+
+}
+
+#endif

src/translating.cpp

@@ -0,0 +1,43 @@
+// -*- C++ -*-
+
+#include "physics/particle.h"
+#include "physics/vector.h"
+#include "translating.h"
+
+using namespace mbostock;
+
+TranslatingModel::TranslatingModel(Model& m, const Translation& t)
+    : model_(m), translation_(t) {
+}
+
+void TranslatingModel::initialize() {
+  model_.initialize();
+}
+
+void TranslatingModel::display() {
+  glPushMatrix();
+  glTranslatev(translation_.origin());
+  model_.display();
+  glPopMatrix();
+}
+
+TranslatingRoomObject::TranslatingRoomObject(RoomObject* o, const Translation& t)
+    : TransformingRoomObject(o), translation_(t),
+      shape_(o->shape(), t), model_(o->model(), t) {
+}
+
+Model& TranslatingRoomObject::model() {
+  return model_;
+}
+
+const Shape& TranslatingRoomObject::shape() const {
+  return shape_;
+}
+
+Vector TranslatingRoomObject::velocity(const Vector& x) const {
+  Vector v = translation_.velocity();
+  if (object_->dynamic()) {
+    v += object_->velocity(translation_.pointInverse(x));
+  }
+  return v;
+}

src/translating.h

@@ -0,0 +1,43 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_TRANSLATING_H
+#define MBOSTOCK_TRANSLATING_H
+
+#include "model.h"
+#include "physics/shape.h"
+#include "physics/translation.h"
+#include "physics/vector.h"
+#include "room_object.h"
+#include "transforming.h"
+
+namespace mbostock {
+
+  class TranslatingModel : public Model {
+  public:
+    TranslatingModel(Model& m, const Translation& t);
+
+    virtual void initialize();
+    virtual void display();
+
+  private:
+    Model& model_;
+    const Translation& translation_;
+  };
+
+  class TranslatingRoomObject : public TransformingRoomObject {
+  public:
+    TranslatingRoomObject(RoomObject* o, const Translation& t);
+
+    virtual Model& model();
+    virtual const Shape& shape() const;
+    virtual Vector velocity(const Vector& x) const;
+
+  private:
+    const Translation& translation_;
+    TranslatingShape shape_;
+    TranslatingModel model_;
+  };
+
+}
+
+#endif

src/tube.cpp

@@ -0,0 +1,35 @@
+// -*- C++ -*-
+
+#include "material.h"
+#include "model.h"
+#include "physics/particle.h"
+#include "physics/shape.h"
+#include "physics/vector.h"
+#include "room.h"
+#include "tube.h"
+
+using namespace mbostock;
+
+Tube::Tube(const Vector& x0, const Vector& x1, const Vector& y, float radius)
+    : cylinder_(x0, x1, radius), y_(y), model_(cylinder_, y_) {
+}
+
+Model& Tube::model() {
+  return model_;
+}
+
+const Shape& Tube::shape() const {
+  return cylinder_;
+}
+
+void Tube::setMaterial(const Material& m) {
+  model_.setMaterial(m);
+}
+
+void Tube::setCapMaterial(const Material& m) {
+  model_.setCapMaterial(m);
+}
+
+float Tube::slip() const {
+  return model_.material().slip();
+}

src/tube.h

@@ -0,0 +1,31 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_TUBE_H
+#define MBOSTOCK_TUBE_H
+
+#include "model.h"
+#include "physics/shape.h"
+#include "room_object.h"
+
+namespace mbostock {
+
+  class Tube : public RoomObject {
+  public:
+    Tube(const Vector& x0, const Vector& x1, const Vector& y, float radius);
+
+    virtual Model& model();
+    virtual const Shape& shape() const;
+    virtual float slip() const;
+
+    void setMaterial(const Material& m);
+    void setCapMaterial(const Material& m);
+
+  private:
+    const Cylinder cylinder_;
+    Vector y_;
+    CylinderModel model_;
+  };
+
+}
+
+#endif

src/wall.cpp

@@ -0,0 +1,60 @@
+// -*- C++ -*-
+
+#include "material.h"
+#include "model.h"
+#include "physics/shape.h"
+#include "room.h"
+#include "wall.h"
+
+using namespace mbostock;
+
+Wall::Wall(const Vector& x0, const Vector& x1,
+           const Vector& x2, const Vector& x3)
+    : quad_(x0, x1, x2, x3), model_(quad_) {
+}
+
+Model& Wall::model() {
+  return model_;
+}
+
+const Shape& Wall::shape() const {
+  return quad_;
+}
+
+void Wall::setMaterial(const Material& m) {
+  model_.setMaterial(m);
+}
+
+float Wall::slip() const {
+  return model_.material().slip();
+}
+
+void Wall::setTexCoords(const Vector& t0, const Vector& t1,
+                        const Vector& t2, const Vector& t3) {
+  model_.setTexCoords(t0, t1, t2, t3);
+}
+
+TriWall::TriWall(const Vector& x0, const Vector& x1, const Vector& x2)
+    : triangle_(x0, x1, x2), model_(triangle_) {
+}
+
+Model& TriWall::model() {
+  return model_;
+}
+
+const Shape& TriWall::shape() const {
+  return triangle_;
+}
+
+void TriWall::setTexCoords(const Vector& t0, const Vector& t1,
+                           const Vector& t2) {
+  model_.setTexCoords(t0, t1, t2);
+}
+
+void TriWall::setMaterial(const Material& m) {
+  model_.setMaterial(m);
+}
+
+float TriWall::slip() const {
+  return model_.material().slip();
+}

src/wall.h

@@ -0,0 +1,51 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_WALL_H
+#define MBOSTOCK_WALL_H
+
+#include "model.h"
+#include "physics/shape.h"
+#include "room_object.h"
+
+namespace mbostock {
+
+  class Material;
+
+  class Wall : public RoomObject {
+  public:
+    Wall(const Vector& x0, const Vector& x1,
+         const Vector& x2, const Vector& x3);
+
+    virtual Model& model();
+    virtual const Shape& shape() const;
+    virtual float slip() const;
+
+    void setTexCoords(const Vector& t0, const Vector& t1,
+                      const Vector& t2, const Vector& t3);
+    void setMaterial(const Material& m);
+
+  private:
+    const Quad quad_;
+    QuadModel model_;
+  };
+
+  class TriWall : public RoomObject {
+  public:
+    TriWall(const Vector& x0, const Vector& x1, const Vector& x2);
+
+    virtual Model& model();
+    virtual const Shape& shape() const;
+    virtual float slip() const;
+
+    void setTexCoords(const Vector& t0, const Vector& t1, const Vector& t2);
+    void setMaterial(const Material& m);
+
+  private:
+    const Triangle triangle_;
+    TriangleModel model_;
+  };
+
+}
+
+
+#endif

src/world.cpp

@@ -0,0 +1,236 @@
+// -*- C++ -*-
+
+#include <OpenGL/gl.h>
+
+#include "material.h"
+#include "portal.h"
+#include "room.h"
+#include "room_force.h"
+#include "room_object.h"
+#include "sound.h"
+#include "trail.h"
+#include "world.h"
+
+using namespace mbostock;
+
+/* Fog. */
+static const float fogDensity = 0.02f;
+static const float fogColor[] = { 0.f, 0.f, 0.f, 1.f };
+
+static const float gravity = 10.f;
+static const float minY = -50.f;
+
+static World* world_ = NULL;
+
+WorldModel::WorldModel(World& world)
+    : world_(world) {
+}
+
+void WorldModel::initialize() {
+  glFogi(GL_FOG_MODE, GL_EXP2);
+  glFogfv(GL_FOG_COLOR, fogColor);
+  glFogf(GL_FOG_DENSITY, fogDensity);
+
+  glEnable(GL_CULL_FACE);
+  glEnable(GL_DEPTH_TEST);
+  glEnable(GL_FOG);
+  glEnable(GL_LIGHTING);
+  glEnable(GL_LINE_SMOOTH);
+  glEnable(GL_TEXTURE_2D);
+
+  std::vector<Room*>::const_iterator i;
+  for (i = world_.rooms().begin(); i != world_.rooms().end(); i++) {
+    (*i)->model().initialize();
+  }
+  world_.player().model().initialize();
+  world_.pauseLighting().initialize();
+}
+
+void WorldModel::display() {
+  world_.room().model().display();
+
+  if (world_.debug()) {
+    std::vector<Trail*>::const_iterator i;
+    for (i = world_.room().trails().begin(); i != world_.room().trails().end();
+         i++) {
+      TrailModel m(**i);
+      m.initialize();
+      m.display();
+    }
+    TrailModel m(world_.room().trail());
+    m.initialize();
+    m.display();
+  }
+
+  world_.player().model().display();
+}
+
+World::World()
+    : Simulation(roundf(ParticleSimulator::timeStep() * 1000.f)),
+      simulator_(1.f), gravity_(gravity), room_(NULL), debug_(false),
+      model_(*this) {
+  world_ = this;
+  pauseLighting_.light(0).setDiffuse(.1f, .1f, .1f, 1.f);
+  pauseLighting_.light(0).setSpecular(.1f, .1f, .1f, 1.f);
+}
+
+World::~World() {
+  std::vector<Room*>::const_iterator ir;
+  for (ir = rooms_.begin(); ir != rooms_.end(); ir++) {
+    delete (*ir);
+  }
+  std::vector<Material*>::const_iterator im;
+  for (im = materials_.begin(); im != materials_.end(); im++) {
+    delete (*im);
+  }
+}
+
+World* World::world() {
+  return world_;
+}
+
+void World::addRoom(Room* r) {
+  if (room_ == NULL) {
+    setRoom(r, r->origins()[0]);
+  }
+  rooms_.push_back(r);
+}
+
+void World::addMaterial(Material* m) {
+  materials_.push_back(m);
+}
+
+void World::addLighting(Lighting* l) {
+  lightings_.push_back(l);
+}
+
+void World::togglePaused() {
+  Simulation::togglePaused();
+  if (paused()) {
+    Sounds::pause();
+  } else {
+    Sounds::resume();
+  }
+}
+
+void World::toggleDebug() {
+  debug_ = !debug_;
+}
+
+void World::setRoom(Room* r, RoomOrigin* origin) {
+  bool sameMusic = false;
+  if ((room_ != NULL) && (room_->music() != NULL)) {
+    if (room_->music() == r->music()) {
+      sameMusic = true;
+    } else {
+      room_->music()->stop();
+    }
+  }
+
+  room_ = r;
+  room_->nextTrail(origin->position());
+  player_.setOrigin(origin->position());
+  player_.setVelocity(origin->velocity());
+
+  if (!sameMusic && (room_->music() != NULL)) {
+    room_->music()->play(-1);
+  }
+}
+
+void World::step() {
+  std::vector<RoomObject*>::const_iterator i;
+
+  /* Reset forces. */
+  player_.resetForces();
+  room_->resetForces();
+
+  /* Apply gravity and contact forces. */
+  player_.applyForce(gravity_);
+  room_->applyForce(gravity_);
+  for (i = contactObjects_.begin(); i != contactObjects_.end(); i++) {
+    (*i)->applyWeight(gravity * player_.mass(), player_.origin());
+  }
+
+  /* Apply localized forces. */
+  std::vector<RoomForce*>::const_iterator f;
+  for (f = room_->forces().begin(); f != room_->forces().end(); f++) {
+    RoomForce& force = **f;
+    if (player_.intersects(force.shape())) {
+      player_.applyForce(force);
+    }
+  }
+
+  /* Run the simulation. */
+  player_.step(simulator_);
+  room_->step(simulator_);
+
+  /* If the player landed on a portal, move to the associated room. */
+  std::vector<Portal*>::const_iterator p;
+  for (p = room_->portals().begin(); p != room_->portals().end(); p++) {
+    Portal& portal = **p;
+    if (portal.contains(player_.origin())) {
+      if (portal.reset()) {
+        room_->reset();
+      }
+      Room* room = rooms_[portal.room()];
+      setRoom(room, room->origins()[portal.origin()]);
+      contactObjects_.clear();
+      return;
+    }
+  }
+
+  /* Apply constraints, detect contacts. */
+  contactObjects_.clear();
+  for (i = room_->objects().begin(); i != room_->objects().end(); i++) {
+    RoomObject& object = **i;
+    object.constrainInternal();
+    if (player_.constrainOutside(object)) {
+      contactObjects_.push_back(&object);
+    }
+  }
+  player_.constrainInternal();
+  room_->trail().add(player_.origin());
+
+  /* Reset if the player falls into a chasm. */
+  if (player_.origin().y < minY) {
+    resetPlayer();
+  }
+}
+
+void World::resetPlayer() {
+  RoomOrigin* origin = room_->origins()[0];
+  player_.setOrigin(origin->position());
+  player_.setVelocity(origin->velocity());
+  room_->reset();
+  room_->nextTrail(origin->position());
+}
+
+void World::nextRoom() {
+  std::vector<Room*>::const_iterator i;
+  Room* nextRoom = rooms_[0];
+  for (i = rooms_.begin(); i != rooms_.end(); i++) {
+    if ((*i) == room_) {
+      i++;
+      if (i != rooms_.end()) {
+        nextRoom = *i;
+      }
+      break;
+    }
+  }
+  setRoom(nextRoom, nextRoom->origins()[0]);
+}
+
+void World::previousRoom() {
+  std::vector<Room*>::reverse_iterator i;
+  Room* nextRoom = rooms_[rooms_.size() - 1];
+  for (i = rooms_.rbegin(); i != rooms_.rend(); i++) {
+    if ((*i) == room_) {
+      i++;
+      if (i != rooms_.rend()) {
+        nextRoom = *i;
+      }
+      break;
+    }
+  }
+  setRoom(nextRoom, nextRoom->origins()[0]);
+}

src/world.h

@@ -0,0 +1,80 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_WORLD_H
+#define MBOSTOCK_WORLD_H
+
+#include <vector>
+
+#include "lighting.h"
+#include "model.h"
+#include "physics/force.h"
+#include "player.h"
+#include "simulation.h"
+
+namespace mbostock {
+
+  class Material;
+  class Room;
+  class RoomObject;
+  class RoomOrigin;
+  class World;
+
+  class WorldModel : public Model {
+  public:
+    WorldModel(World& world);
+
+    virtual void initialize();
+    virtual void display();
+
+  private:
+    World& world_;
+  };
+
+  class World : public Simulation {
+  public:
+    World();
+    virtual ~World();
+
+    static World* world();
+
+    void addRoom(Room* r);
+    void addMaterial(Material* m);
+    void addLighting(Lighting* l);
+
+    inline Player& player() { return player_; }
+    inline const std::vector<Room*>& rooms() const { return rooms_; }
+    inline Room& room() const { return *room_; }
+    inline Model& model() { return model_; }
+    inline const Lighting& pauseLighting() const { return pauseLighting_; }
+
+    void resetPlayer();
+    void nextRoom();
+    void previousRoom();
+
+    virtual void togglePaused();
+
+    void toggleDebug();
+    inline bool debug() const { return debug_; }
+
+  protected:
+    virtual void step();
+
+  private:
+    void setRoom(Room* room, RoomOrigin* origin);
+
+    ParticleSimulator simulator_;
+    GravitationalForce gravity_;
+    Player player_;
+    Lighting pauseLighting_;
+    std::vector<Lighting*> lightings_;
+    std::vector<Material*> materials_;
+    std::vector<Room*> rooms_;
+    std::vector<RoomObject*> contactObjects_;
+    Room* room_;
+    bool debug_;
+    WorldModel model_;
+  };
+
+}
+
+#endif

src/worlds.cpp

@@ -0,0 +1,798 @@
+#include <TinyXML/tinyxml.h>
+#include <iostream>
+#include <list>
+#include <map>
+#include <math.h>
+#include <string>
+#include <vector>
+
+#include "ball.h"
+#include "block.h"
+#include "escalator.h"
+#include "fan.h"
+#include "lighting.h"
+#include "material.h"
+#include "portal.h"
+#include "ramp.h"
+#include "resource.h"
+#include "room.h"
+#include "room_force.h"
+#include "rotating.h"
+#include "seesaw.h"
+#include "sound.h"
+#include "switch.h"
+#include "translating.h"
+#include "tube.h"
+#include "wall.h"
+#include "world.h"
+#include "worlds.h"
+
+namespace mbostock {
+
+  class Portal;
+  class RoomForce;
+  class RoomObject;
+  class RoomOrigin;
+  class World;
+
+  class SwitchTargets {
+  public:
+    SwitchTargets(Switch* s);
+
+    inline Switch* svitch() const { return switch_; }
+
+    void addTarget(const char* targetName);
+    inline const std::vector<std::string>& targets() const { return targets_; }
+
+  private:
+    Switch* switch_;
+    std::vector<std::string> targets_;
+  };
+
+  class XmlWorldBuilder {
+  public:
+    XmlWorldBuilder();
+
+    World* parseWorld(const char* path);
+
+  private:
+    static Vector parseVector(TiXmlElement* e, const Vector& d);
+    static Vector parseVector(TiXmlElement* e);
+    static bool parseBool(TiXmlElement* e, const char* name, bool d);
+    static bool parseBool(TiXmlElement* e, const char* name);
+
+    void parseLightings(TiXmlElement* e);
+    void parseLighting(TiXmlElement* e);
+    void parseLightGlobalAmbient(Lighting* l, TiXmlElement* e);
+    void parseLights(Lighting* l, TiXmlElement* e);
+    void parseLight(Light& l, TiXmlElement* e);
+    void parseLightAmbient(Light& l, TiXmlElement* e);
+    void parseLightDiffuse(Light& l, TiXmlElement* e);
+    void parseLightSpecular(Light& l, TiXmlElement* e);
+    void parseLightPosition(Light& l, TiXmlElement* e);
+    void parseLightSpotDirection(Light& l, TiXmlElement* e);
+    void parseLightAttributes(Light& l, TiXmlElement* e);
+
+    void parseMaterials(TiXmlElement* e);
+    void parseMaterial(TiXmlElement* e);
+    void parseMaterialParameter(Material* m, TiXmlElement* e);
+
+    void parseRooms(TiXmlElement* e);
+    void parseRoom(TiXmlElement* e);
+    void parseRoomLighting(Room* r, TiXmlElement* e);
+    void parseRoomMusic(Room* r, TiXmlElement* e);
+    void parseRoomCameraBounds(Room* r, TiXmlElement* e);
+    void parseRoomTopLevelObjects(Room* r, TiXmlElement* e);
+    void parseRoomObjects(Room* r, TiXmlElement* e);
+
+    RoomOrigin* parseRoomOrigin(TiXmlElement* e);
+    Portal* parseRoomPortal(TiXmlElement* e);
+
+    RoomObject* parseRoomObject(TiXmlElement* e);
+    RoomObject* parseRoomBlock(TiXmlElement* e);
+    RoomObject* parseRoomAxisAlignedBlock(TiXmlElement* e);
+    RoomObject* parseRoomOrientedBlock(TiXmlElement* e);
+    RoomObject* parseRoomWall(TiXmlElement* e);
+    RoomObject* parseRoomQuadWall(TiXmlElement* e);
+    RoomObject* parseRoomTriWall(TiXmlElement* e);
+    RoomObject* parseRoomEscalator(TiXmlElement* e);
+    RoomObject* parseRoomSeesaw(TiXmlElement* e);
+    RoomObject* parseRoomRamp(TiXmlElement* e);
+    RoomObject* parseRoomTube(TiXmlElement* e);
+    RoomObject* parseRoomBall(TiXmlElement* e);
+    RoomObject* parseRoomFan(TiXmlElement* e);
+    RoomObject* parseRoomSwitch(TiXmlElement* e);
+
+    RoomForce* parseRoomConstantForce(TiXmlElement* e);
+
+    void parseRotation(Room* r, TiXmlElement* e);
+    void parseTranslation(Room* r, TiXmlElement* e);
+    RoomObject* applyTransforms(RoomObject* o);
+
+    void parseRoomSwitchTargets(Switch* s, TiXmlElement* e);
+    void resolveSwitchTargets();
+
+    int findRoom(const char* name);
+    int findRoomOrigin(const char* name);
+
+    TiXmlDocument document_;
+    World* world_;
+    std::map<std::string, Lighting*> lightings_;
+    std::map<std::string, Material*> materials_;
+    std::map<std::string, Transform*> transforms_;
+    std::list<Transform*> activeTransforms_;
+    std::vector<SwitchTargets*> switchTargets_;
+  };
+
+}
+
+using namespace mbostock;
+
+SwitchTargets::SwitchTargets(Switch* s)
+    : switch_(s) {
+}
+
+void SwitchTargets::addTarget(const char* targetName) {
+  targets_.push_back(targetName);
+}
+
+XmlWorldBuilder::XmlWorldBuilder()
+    : world_(NULL) {
+}
+
+World* XmlWorldBuilder::parseWorld(const char* path) {
+  std::string fullPath(Resources::path());
+  fullPath.append(path);
+  if (!document_.LoadFile(fullPath.c_str())) {
+    std::cerr << "Error loading world \"" << path << "\": ";
+    std::cerr << document_.ErrorDesc() << "\n";
+    return NULL;
+  }
+
+  world_ = new World();
+  TiXmlElement* e = document_.FirstChildElement("world");
+  parseLightings(e);
+  parseMaterials(e);
+  parseRooms(e);
+  resolveSwitchTargets();
+  return world_;
+}
+
+void XmlWorldBuilder::parseLightings(TiXmlElement* e) {
+  for (TiXmlElement* l = e->FirstChildElement("lighting"); l != NULL;
+       l = l->NextSiblingElement("lighting")) {
+    parseLighting(l);
+  }
+}
+
+void XmlWorldBuilder::parseLighting(TiXmlElement* e) {
+  Lighting* l = new Lighting();
+  lightings_[e->Attribute("name")] = l;
+  parseLightGlobalAmbient(l, e->FirstChildElement("ambient"));
+  parseLights(l, e);
+  world_->addLighting(l);
+}
+
+void XmlWorldBuilder::parseLightGlobalAmbient(Lighting* l, TiXmlElement* e) {
+  if (e != NULL) {
+    float r = 0.f, g = 0.f, b = 0.f, a = 1.f;
+    e->QueryFloatAttribute("r", &r);
+    e->QueryFloatAttribute("g", &g);
+    e->QueryFloatAttribute("b", &b);
+    e->QueryFloatAttribute("a", &a);
+    l->setGlobalAmbient(r, g, b, a);
+  }
+}
+
+void XmlWorldBuilder::parseLights(Lighting* l, TiXmlElement* e) {
+  int i = 0;
+  for (TiXmlElement* le = e->FirstChildElement("light");
+       le != NULL; le = le->NextSiblingElement("light")) {
+    parseLight(l->light(i++), le);
+  }
+}
+
+void XmlWorldBuilder::parseLight(Light& l, TiXmlElement* e) {
+  parseLightAmbient(l, e->FirstChildElement("ambient"));
+  parseLightSpecular(l, e->FirstChildElement("specular"));
+  parseLightDiffuse(l, e->FirstChildElement("diffuse"));
+  parseLightPosition(l, e->FirstChildElement("position"));
+  parseLightSpotDirection(l, e->FirstChildElement("spot-direction"));
+  parseLightAttributes(l, e);
+  l.enable();
+}
+
+void XmlWorldBuilder::parseLightAmbient(Light& l, TiXmlElement* e) {
+  if (e != NULL) {
+    float r = 0.f, g = 0.f, b = 0.f, a = 1.f;
+    e->QueryFloatAttribute("r", &r);
+    e->QueryFloatAttribute("g", &g);
+    e->QueryFloatAttribute("b", &b);
+    e->QueryFloatAttribute("a", &a);
+    l.setAmbient(r, g, b, a);
+  }
+}
+
+void XmlWorldBuilder::parseLightSpecular(Light& l, TiXmlElement* e) {
+  if (e != NULL) {
+    float r = 0.f, g = 0.f, b = 0.f, a = 1.f;
+    e->QueryFloatAttribute("r", &r);
+    e->QueryFloatAttribute("g", &g);
+    e->QueryFloatAttribute("b", &b);
+    e->QueryFloatAttribute("a", &a);
+    l.setSpecular(r, g, b, a);
+  }
+}
+
+void XmlWorldBuilder::parseLightDiffuse(Light& l, TiXmlElement* e) {
+  if (e != NULL) {
+    float r = 0.f, g = 0.f, b = 0.f, a = 1.f;
+    e->QueryFloatAttribute("r", &r);
+    e->QueryFloatAttribute("g", &g);
+    e->QueryFloatAttribute("b", &b);
+    e->QueryFloatAttribute("a", &a);
+    l.setDiffuse(r, g, b, a);
+  }
+}
+
+void XmlWorldBuilder::parseLightPosition(Light& l, TiXmlElement* e) {
+  if (e != NULL) {
+    float x, y, z, w;
+    e->QueryFloatAttribute("x", &x);
+    e->QueryFloatAttribute("y", &y);
+    e->QueryFloatAttribute("z", &z);
+    e->QueryFloatAttribute("w", &w);
+    l.setPosition(x, y, z, w);
+  }
+}
+
+void XmlWorldBuilder::parseLightSpotDirection(Light& l, TiXmlElement* e) {
+  if (e != NULL) {
+    float x, y, z;
+    e->QueryFloatAttribute("x", &x);
+    e->QueryFloatAttribute("y", &y);
+    e->QueryFloatAttribute("z", &z);
+    l.setSpotDirection(x, y, z);
+  }
+}
+
+void XmlWorldBuilder::parseLightAttributes(Light& l, TiXmlElement* e) {
+  if (e->Attribute("spot-exponent")) {
+    float f;
+    e->QueryFloatAttribute("spot-exponent", &f);
+    l.setSpotExponent(f);
+  }
+  if (e->Attribute("constant-attenuation")) {
+    float f;
+    e->QueryFloatAttribute("constant-attenuation", &f);
+    l.setConstantAttenuation(f);
+  }
+  if (e->Attribute("linear-attenuation")) {
+    float f;
+    e->QueryFloatAttribute("linear-attenuation", &f);
+    l.setLinearAttenuation(f);
+  }
+  if (e->Attribute("quadratic-attenuation")) {
+    float f;
+    e->QueryFloatAttribute("quadratic-attenuation", &f);
+    l.setQuadraticAttenuation(f);
+  }
+}
+
+void XmlWorldBuilder::parseMaterials(TiXmlElement* e) {
+  for (TiXmlElement* m = e->FirstChildElement("material"); m != NULL;
+       m = m->NextSiblingElement("material")) {
+    parseMaterial(m);
+  }
+}
+
+void XmlWorldBuilder::parseMaterial(TiXmlElement* e) {
+  Material* m = new Material();
+  materials_[e->Attribute("name")] = m;
+  float s = 90.f;
+  e->QueryFloatAttribute("slip-angle", &s);
+  m->setSlipAngle(s);
+  for (TiXmlElement* p = e->FirstChildElement(); p != NULL;
+       p = p->NextSiblingElement()) {
+    parseMaterialParameter(m, p);
+  }
+  world_->addMaterial(m);
+}
+
+void XmlWorldBuilder::parseMaterialParameter(Material* m, TiXmlElement* e) {
+  const std::string& name = e->ValueStr();
+  if (name == "texture") {
+    m->setTexture(e->Attribute("path"));
+    return;
+  }
+
+  float r = 0.f, g = 0.f, b = 0.f, a = 1.f;
+  e->QueryFloatAttribute("r", &r);
+  e->QueryFloatAttribute("g", &g);
+  e->QueryFloatAttribute("b", &b);
+  e->QueryFloatAttribute("a", &a);
+
+  if (name == "ambient") {
+    m->setAmbient(r, g, b, a);
+  } else if (name == "diffuse") {
+    m->setDiffuse(r, g, b, a);
+  } else if (name == "emission") {
+    m->setEmission(r, g, b, a);
+  } else if (name == "specular") {
+    m->setSpecular(r, g, b, a);
+  }
+}
+
+void XmlWorldBuilder::parseRooms(TiXmlElement* e) {
+  for (TiXmlElement* r = e->FirstChildElement("room"); r != NULL;
+       r = r->NextSiblingElement("room")) {
+    parseRoom(r);
+  }
+}
+
+void XmlWorldBuilder::parseRoom(TiXmlElement* e) {
+  Room* r = new Room();
+  parseRoomLighting(r, e);
+  parseRoomMusic(r, e);
+  parseRoomCameraBounds(r, e);
+  parseRoomTopLevelObjects(r, e);
+  world_->addRoom(r);
+}
+
+void XmlWorldBuilder::parseRoomLighting(Room* r, TiXmlElement* e) {
+  const char* lighting = e->Attribute("lighting");
+  if (lighting != NULL) {
+    r->setLighting(*lightings_[lighting]);
+  }
+}
+
+void XmlWorldBuilder::parseRoomMusic(Room* r, TiXmlElement* e) {
+  const char* music = e->Attribute("music");
+  if (music != NULL) {
+    r->setMusic(Sounds::fromFile(music));
+  }
+}
+
+void XmlWorldBuilder::parseRoomCameraBounds(Room* r, TiXmlElement* e) {
+  Vector min = parseVector(e->FirstChildElement("camera-min"), -Vector::INF());
+  Vector max = parseVector(e->FirstChildElement("camera-max"), Vector::INF());
+  r->setCameraBounds(min, max);
+}
+
+void XmlWorldBuilder::parseRoomTopLevelObjects(Room* r, TiXmlElement* e) {
+  for (TiXmlElement* o = e->FirstChildElement(); o != NULL;
+       o = o->NextSiblingElement()) {
+    const std::string& name = o->ValueStr();
+    if (name == "origin") {
+      r->addOrigin(parseRoomOrigin(o));
+    } else if (name == "portal") {
+      r->addPortal(parseRoomPortal(o));
+    }
+  }
+  parseRoomObjects(r, e);
+}
+
+void XmlWorldBuilder::parseRoomObjects(Room* r, TiXmlElement* e) {
+  for (TiXmlElement* o = e->FirstChildElement(); o != NULL;
+       o = o->NextSiblingElement()) {
+    const std::string& name = o->ValueStr();
+    if (name == "constant-force") {
+      r->addForce(parseRoomConstantForce(o));
+    } else if (name == "rotation") {
+      parseRotation(r, o);
+    } else if (name == "translation") {
+      parseTranslation(r, o);
+    } else {
+      RoomObject* ro = parseRoomObject(o);
+      if (ro != NULL) { // ignore unknown objects
+        r->addObject(applyTransforms(ro));
+      }
+    }
+  }
+}
+
+void XmlWorldBuilder::parseRotation(Room* r, TiXmlElement* e) {
+  float s = 10.f, a = 0.f;
+  e->QueryFloatAttribute("speed", &s);
+  e->QueryFloatAttribute("angle", &a);
+  Rotation* z = new Rotation(
+      parseVector(e->FirstChildElement("origin")),
+      parseVector(e->FirstChildElement("axis")),
+      s, a);
+  const char* name = e->Attribute("name");
+  if (name != NULL) {
+    transforms_[name] = z;
+  }
+  activeTransforms_.push_back(z);
+  parseRoomObjects(r, e);
+  activeTransforms_.pop_back();
+  r->addTransform(z);
+}
+
+void XmlWorldBuilder::parseTranslation(Room* r, TiXmlElement* e) {
+  float s = 10.f, u = 0.f, kd = 0.f;
+  e->QueryFloatAttribute("speed", &s);
+  e->QueryFloatAttribute("start", &u);
+  e->QueryFloatAttribute("dampen", &kd);
+  Translation* z = new Translation(
+      parseVector(e->FirstChildElement("x0")),
+      parseVector(e->FirstChildElement("x1")),
+      s, u, kd);
+  const char* mode = e->Attribute("mode");
+  if (mode != NULL) {
+    std::string modeString = mode;
+    if (modeString == "one-way") {
+      z->setMode(Translation::ONE_WAY);
+    } else if (modeString == "reset") {
+      z->setMode(Translation::RESET);
+    }
+  }
+  const char* name = e->Attribute("name");
+  if (name != NULL) {
+    transforms_[name] = z;
+  }
+  activeTransforms_.push_back(z);
+  parseRoomObjects(r, e);
+  activeTransforms_.pop_back();
+  r->addTransform(z);
+}
+
+RoomObject* XmlWorldBuilder::applyTransforms(RoomObject* o) {
+  std::list<Transform*>::reverse_iterator i;
+  for (i = activeTransforms_.rbegin(); i != activeTransforms_.rend(); i++) {
+    Transform* z = *i;
+    Translation* t = dynamic_cast<Translation*>(z);
+    if (t != NULL) {
+      o = new TranslatingRoomObject(o, *t);
+      continue;
+    }
+    Rotation* r = dynamic_cast<Rotation*>(z);
+    if (r != NULL) {
+      o = new RotatingRoomObject(o, *r);
+      continue;
+    }
+  }
+  return o;
+}
+
+RoomOrigin* XmlWorldBuilder::parseRoomOrigin(TiXmlElement* e) {
+  return new RoomOrigin(
+      parseVector(e->FirstChildElement("position")),
+      parseVector(e->FirstChildElement("velocity")));
+}
+
+Portal* XmlWorldBuilder::parseRoomPortal(TiXmlElement* e) {
+  return new Portal(
+      parseVector(e->FirstChildElement("min")),
+      parseVector(e->FirstChildElement("max")),
+      findRoom(e->Attribute("origin")),
+      findRoomOrigin(e->Attribute("origin")),
+      parseBool(e, "reset", false));
+}
+
+RoomObject* XmlWorldBuilder::parseRoomObject(TiXmlElement* e) {
+  const std::string& name = e->ValueStr();
+  if (name == "block") {
+    return parseRoomBlock(e);
+  } else if (name == "wall") {
+    return parseRoomWall(e);
+  } else if (name == "escalator") {
+    return parseRoomEscalator(e);
+  } else if (name == "seesaw") {
+    return parseRoomSeesaw(e);
+  } else if (name == "ramp") {
+    return parseRoomRamp(e);
+  } else if (name == "tube") {
+    return parseRoomTube(e);
+  } else if (name == "ball") {
+    return parseRoomBall(e);
+  } else if (name == "fan") {
+    return parseRoomFan(e);
+  } else if (name == "switch") {
+    return parseRoomSwitch(e);
+  }
+  return NULL;
+}
+
+RoomObject* XmlWorldBuilder::parseRoomBlock(TiXmlElement* e) {
+  return (e->FirstChildElement("min") != NULL)
+      ? parseRoomAxisAlignedBlock(e)
+      : parseRoomOrientedBlock(e);
+}
+
+RoomObject* XmlWorldBuilder::parseRoomAxisAlignedBlock(TiXmlElement* e) {
+  AxisAlignedBlock* b = new AxisAlignedBlock(
+      parseVector(e->FirstChildElement("min")),
+      parseVector(e->FirstChildElement("max")));
+  const char* material = e->Attribute("material");
+  if (material != NULL) {
+    b->setMaterial(*materials_[material]);
+  }
+  const char* topMaterial = e->Attribute("top-material");
+  if (topMaterial != NULL) {
+    b->setTopMaterial(*materials_[topMaterial]);
+  }
+  return b;
+}
+
+RoomObject* XmlWorldBuilder::parseRoomOrientedBlock(TiXmlElement* e) {
+  Block* b = new Block(
+      parseVector(e->FirstChildElement("c")),
+      parseVector(e->FirstChildElement("x")),
+      parseVector(e->FirstChildElement("y")),
+      parseVector(e->FirstChildElement("z")));
+  const char* material = e->Attribute("material");
+  if (material != NULL) {
+    b->setMaterial(*materials_[material]);
+  }
+  const char* topMaterial = e->Attribute("top-material");
+  if (topMaterial != NULL) {
+    b->setTopMaterial(*materials_[topMaterial]);
+  }
+  return b;
+}
+
+RoomObject* XmlWorldBuilder::parseRoomWall(TiXmlElement* e) {
+  return (e->FirstChildElement("x3") == NULL)
+      ? parseRoomTriWall(e)
+      : parseRoomQuadWall(e);
+}
+
+RoomObject* XmlWorldBuilder::parseRoomQuadWall(TiXmlElement* e) {
+  Wall* w = new Wall(
+      parseVector(e->FirstChildElement("x0")),
+      parseVector(e->FirstChildElement("x1")),
+      parseVector(e->FirstChildElement("x2")),
+      parseVector(e->FirstChildElement("x3")));
+  TiXmlElement* t = e->FirstChildElement("tex-coords");
+  if (t != NULL) {
+    w->setTexCoords(
+        parseVector(t->FirstChildElement("t0")),
+        parseVector(t->FirstChildElement("t1")),
+        parseVector(t->FirstChildElement("t2")),
+        parseVector(t->FirstChildElement("t3")));
+  }
+  const char* material = e->Attribute("material");
+  if (material != NULL) {
+    w->setMaterial(*materials_[material]);
+  }
+  return w;
+}
+
+RoomObject* XmlWorldBuilder::parseRoomTriWall(TiXmlElement* e) {
+  TriWall* w = new TriWall(
+      parseVector(e->FirstChildElement("x0")),
+      parseVector(e->FirstChildElement("x1")),
+      parseVector(e->FirstChildElement("x2")));
+  TiXmlElement* t = e->FirstChildElement("tex-coords");
+  if (t != NULL) {
+    w->setTexCoords(
+        parseVector(t->FirstChildElement("t0")),
+        parseVector(t->FirstChildElement("t1")),
+        parseVector(t->FirstChildElement("t2")));
+  }
+  const char* material = e->Attribute("material");
+  if (material != NULL) {
+    w->setMaterial(*materials_[material]);
+  }
+  return w;
+}
+
+RoomObject* XmlWorldBuilder::parseRoomEscalator(TiXmlElement* e) {
+  Escalator* o = new Escalator(
+      parseVector(e->FirstChildElement("min")),
+      parseVector(e->FirstChildElement("max")),
+      parseVector(e->FirstChildElement("v")));
+  const char* material = e->Attribute("material");
+  if (material != NULL) {
+    o->setMaterial(*materials_[material]);
+  }
+  const char* topMaterial = e->Attribute("top-material");
+  if (topMaterial != NULL) {
+    o->setTopMaterial(*materials_[topMaterial]);
+  }
+  return o;
+}
+
+RoomObject* XmlWorldBuilder::parseRoomSeesaw(TiXmlElement* e) {
+  float mass = 1.f;
+  e->QueryFloatAttribute("mass", &mass);
+  Seesaw* s = new Seesaw(
+      parseVector(e->FirstChildElement("min")),
+      parseVector(e->FirstChildElement("max")),
+      mass);
+  const char* material = e->Attribute("material");
+  if (material != NULL) {
+    s->setMaterial(*materials_[material]);
+  }
+  const char* topMaterial = e->Attribute("top-material");
+  if (topMaterial != NULL) {
+    s->setTopMaterial(*materials_[topMaterial]);
+  }
+  return s;
+}
+
+RoomObject* XmlWorldBuilder::parseRoomRamp(TiXmlElement* e) {
+  Ramp* r = new Ramp(
+      parseVector(e->FirstChildElement("x0")),
+      parseVector(e->FirstChildElement("x1")),
+      parseVector(e->FirstChildElement("x2")),
+      parseVector(e->FirstChildElement("x3")));
+  const char* material = e->Attribute("material");
+  if (material != NULL) {
+    r->setMaterial(*materials_[material]);
+  }
+  const char* topMaterial = e->Attribute("top-material");
+  if (topMaterial != NULL) {
+    r->setTopMaterial(*materials_[topMaterial]);
+  }
+  return r;
+}
+
+RoomObject* XmlWorldBuilder::parseRoomTube(TiXmlElement* e) {
+  float r = 1.f;
+  e->QueryFloatAttribute("radius", &r);
+  Tube* t = new Tube(
+      parseVector(e->FirstChildElement("x0")),
+      parseVector(e->FirstChildElement("x1")),
+      parseVector(e->FirstChildElement("y"), Vector::Y()),
+      r);
+  const char* material = e->Attribute("material");
+  if (material != NULL) {
+    t->setMaterial(*materials_[material]);
+  }
+  const char* capMaterial = e->Attribute("cap-material");
+  if (capMaterial != NULL) {
+    t->setCapMaterial(*materials_[capMaterial]);
+  }
+  return t;
+}
+
+RoomObject* XmlWorldBuilder::parseRoomBall(TiXmlElement* e) {
+  float r = 1.f;
+  e->QueryFloatAttribute("radius", &r);
+  Ball* b = new Ball(parseVector(e->FirstChildElement("x")), r);
+  const char* material = e->Attribute("material");
+  if (material != NULL) {
+    b->setMaterial(*materials_[material]);
+  }
+  return b;
+}
+
+RoomObject* XmlWorldBuilder::parseRoomFan(TiXmlElement* e) {
+  float r, s;
+  e->QueryFloatAttribute("radius", &r);
+  e->QueryFloatAttribute("speed", &s);
+  Fan* f = new Fan(
+      parseVector(e->FirstChildElement("x")),
+      parseVector(e->FirstChildElement("v")),
+      r, s);
+  const char* material = e->Attribute("material");
+  if (material != NULL) {
+    f->setMaterial(*materials_[material]);
+  }
+  return f;
+}
+
+RoomObject* XmlWorldBuilder::parseRoomSwitch(TiXmlElement* e) {
+  Switch* b = new Switch(
+      parseVector(e->FirstChildElement("min")),
+      parseVector(e->FirstChildElement("max")));
+  const char* material = e->Attribute("material");
+  if (material != NULL) {
+    b->setMaterial(*materials_[material]);
+  }
+  const char* topMaterial = e->Attribute("top-material");
+  if (topMaterial != NULL) {
+    b->setTopMaterial(*materials_[topMaterial]);
+  }
+  const char* activeMaterial = e->Attribute("active-material");
+  if (activeMaterial != NULL) {
+    b->setActiveMaterial(*materials_[activeMaterial]);
+  }
+  parseRoomSwitchTargets(b, e);
+  return b;
+}
+
+void XmlWorldBuilder::parseRoomSwitchTargets(Switch* s, TiXmlElement* e) {
+  SwitchTargets* q = new SwitchTargets(s);
+  for (TiXmlElement* t = e->FirstChildElement("target"); t != NULL;
+       t = t->NextSiblingElement("target")) {
+    q->addTarget(t->Attribute("name"));
+  }
+  switchTargets_.push_back(q);
+}
+
+void XmlWorldBuilder::resolveSwitchTargets() {
+  std::vector<SwitchTargets*>::const_iterator i;
+  for (i = switchTargets_.begin(); i != switchTargets_.end(); i++) {
+    SwitchTargets* t = *i;
+    Switch* s = t->svitch();
+    std::vector<std::string>::const_iterator is;
+    for (is = t->targets().begin(); is != t->targets().end(); is++) {
+      s->addTarget(*transforms_[*is]);
+    }
+  }
+  switchTargets_.clear();
+}
+
+RoomForce* XmlWorldBuilder::parseRoomConstantForce(TiXmlElement* e) {
+  return new ConstantRoomForce(
+      parseVector(e->FirstChildElement("min")),
+      parseVector(e->FirstChildElement("max")),
+      parseVector(e->FirstChildElement("force")));
+}
+
+int XmlWorldBuilder::findRoom(const char* name) {
+  int i = 0;
+  std::string s = name;
+  std::string roomName = s.substr(0, s.find('.'));
+  TiXmlElement* e = document_.FirstChildElement("world");
+  for (TiXmlElement* r = e->FirstChildElement("room"); r != NULL;
+       r = r->NextSiblingElement("room"), i++) {
+    const char* t = r->Attribute("name");
+    if ((t != NULL) && (roomName == t)) {
+      return i;
+    }
+  }
+  std::cerr << "Error: could not find room " << name << "\n";
+  return -1;
+}
+
+int XmlWorldBuilder::findRoomOrigin(const char* name) {
+  std::string s = name;
+  int i = s.find('.');
+  std::string roomName = s.substr(0, i);
+  std::string originName = s.substr(i + 1);
+  TiXmlElement* e = document_.FirstChildElement("world");
+  for (TiXmlElement* r = e->FirstChildElement("room"); r != NULL;
+       r = r->NextSiblingElement("room")) {
+    const char* t = r->Attribute("name");
+    if ((t != NULL) && (roomName == t)) {
+      int j = 0;
+      for (TiXmlElement* o = r->FirstChildElement("origin"); o != NULL;
+           o = o->NextSiblingElement("origin"), j++) {
+        const char* u = o->Attribute("name");
+        if ((u != NULL) && (originName == u)) {
+          return j;
+        }
+      }
+      std::cerr << "Error: could not find origin " << roomName << "."
+                << originName << "\n";
+      return -1;
+    }
+  }
+  std::cerr << "Error: could not find room " << roomName << "\n";
+  return -1;
+}
+
+Vector XmlWorldBuilder::parseVector(TiXmlElement* e, const Vector& d) {
+  return (e == NULL) ? d : parseVector(e);
+}
+
+Vector XmlWorldBuilder::parseVector(TiXmlElement* e) {
+  Vector v;
+  e->QueryFloatAttribute("x", &v.x);
+  e->QueryFloatAttribute("y", &v.y);
+  e->QueryFloatAttribute("z", &v.z);
+  return v;
+}
+
+bool XmlWorldBuilder::parseBool(TiXmlElement* e, const char* name) {
+  static const std::string TRUE = "true";
+  return TRUE == e->Attribute(name);
+}
+
+bool XmlWorldBuilder::parseBool(TiXmlElement* e, const char* name, bool d) {
+  const char* value = e->Attribute(name);
+  if (value == NULL) {
+    return d;
+  }
+  static const std::string TRUE = "true";
+  return TRUE == value;
+}
+
+World* Worlds::fromFile(const char* path) {
+  XmlWorldBuilder builder;
+  return builder.parseWorld(path);
+}

src/worlds.h

@@ -0,0 +1,20 @@
+// -*- C++ -*-
+
+#ifndef MBOSTOCK_WORLDS_H
+#define MBOSTOCK_WORLDS_H
+
+namespace mbostock {
+
+  class World;
+
+  class Worlds {
+  public:
+    static World* fromFile(const char* path);
+
+  private:
+    Worlds();
+  };
+
+}
+
+#endif