panning
[ld42_outofspace] / src / gamescr.cc
index d51c16a..3eec738 100644 (file)
@@ -6,10 +6,11 @@
 #include "opengl.h"
 #include "texture.h"
 #include "sdr.h"
+#include "mesh.h"
+#include "meshgen.h"
 
 /* NOTES:
  * - whistle hhgg music
- * - colliding particles merge
  * - select objects and center camera on them
  */
 
@@ -18,6 +19,10 @@ struct Particle {
        float mass;
        Vec2 pos;
        Vec2 vel;
+
+       float vis_height;
+
+       Particle *next, *prev;
 };
 
 struct Emitter {
@@ -25,6 +30,8 @@ struct Emitter {
        float mass;
        float rate, chunk;
        float angle, spread;
+
+       float spawn_pending;
 };
 
 struct Rect {
@@ -32,10 +39,20 @@ struct Rect {
        int width, height;
 };
 
+struct QuadMesh {
+       Vec3 *v;
+       Vec2 *uv;
+       uint16_t *idx;
+
+       int num_verts, num_idx, num_quads;
+
+       unsigned int vbo_v, vbo_uv, ibo;
+};
+
 #define SIM_DT         0.016
 
 #define GRID_SIZE      2048
-#define GRID_BITS      12
+#define GRID_BITS      11
 #define GRID_X(idx)    (((idx) >> GRID_BITS) & (GRID_SIZE - 1))
 #define GRID_Y(idx)    ((idx) & (GRID_SIZE - 1))
 #define GRID_DELTA     ((float)FIELD_SIZE / (float)GRID_SIZE)
@@ -44,28 +61,62 @@ struct Rect {
 #define MIN_CAM_DIST   1.0f
 #define MAX_CAM_DIST   350.0f
 
-#define MASS_TO_RAD(m) log((m) + 1.0)
+#define MASS_TO_RADIUS(m)      log((m) + 1.0)
 
 #define CONTRIB_THRES  0.005
 #define CONTRIB_RANGE(m) sqrt((m) / CONTRIB_THRES)
 
+/* gravitational strength */
+#define GRAV_STR       16.0f
+
+static int pos_to_grid_x_noclamp(float x);
+static int pos_to_grid_y_noclamp(float y);
 static int pos_to_grid(float x, float y);
 static Vec2 grid_to_pos(int gx, int gy);
 
+static void calc_contrib_bounds(const Vec2 &pos, float mass, Rect *rect);
+static void add_influence(const Vec2 &pos, float mass, float radius, const Rect &cbox);
+
+static Vec2 calc_field_grad(int gidx);
+
+static void destroy_quadmesh(QuadMesh *m);
+static void draw_quadmesh(const QuadMesh *m, unsigned int prim = GL_QUADS);
+static void gen_quad_plane(QuadMesh *mesh, float width, float height, int usub, int vsub);
+
+static void spawn_particle(const Vec2 &pos, const Vec2 &vel, float mass);
+static void add_particle(Particle *p);
+static void remove_particle(Particle *p);
+static Particle *alloc_particle();
+void free_particle(Particle *p);
+
+static bool pause;
+
 static float grid[GRID_SIZE * GRID_SIZE];
-static Particle grid_part[GRID_SIZE * GRID_SIZE];
+static Particle *grid_part[GRID_SIZE * GRID_SIZE];
 static Texture *grid_tex;
 
+static Particle *plist;
+
 static std::vector<Emitter*> emitters;
 
 static Texture *gvis_tex;      // texture tile for visualizing a grid
 static unsigned int field_sdr;
+static int tess_level = 64;
+static float field_scale = 16.0f;
+
+static QuadMesh field_mesh;
 
+static Mesh *pmesh;
+static unsigned int particle_sdr;
+
+static Vec2 cam_pos;
 static float cam_theta;
 static float cam_dist = 100.0f;
-static Vec2 *targ_pos;
+static Vec2 *targ_pos = &cam_pos;
 static Mat4 view_matrix, proj_matrix;
 
+static bool wireframe;
+
 // emitter placement data (filled by event handlers, completed in update)
 static bool emit_place_pending;
 static Vec2 emit_place_pos;
@@ -97,78 +148,155 @@ bool GameScreen::init()
        }
        set_uniform_int(field_sdr, "gvis_tex", 0);
        set_uniform_int(field_sdr, "field_tex", 1);
-       set_uniform_int(field_sdr, "tess_level", 2);
+       set_uniform_float(field_sdr, "gvis_scale", FIELD_SIZE / 32.0f);
+       set_uniform_int(field_sdr, "tess_level", tess_level);
+       set_uniform_float(field_sdr, "field_scale", field_scale);
+
+       gen_quad_plane(&field_mesh, FIELD_SIZE, FIELD_SIZE, 32, 32);
+
+       pmesh = new Mesh;
+       gen_geosphere(pmesh, 1.0, 2);
+
+       if(!(particle_sdr = create_program_load("sdr/sph.v.glsl", "sdr/sph.p.glsl"))) {
+               return false;
+       }
 
        // XXX DBG
        emit_place_pos = Vec2(0, 0);
        emit_place_pending = true;
 
+       assert(glGetError() == GL_NO_ERROR);
        return true;
 }
 
 void GameScreen::destroy()
 {
+       free_program(field_sdr);
        delete gvis_tex;
+       delete grid_tex;
+       destroy_quadmesh(&field_mesh);
+       delete pmesh;
 }
 
-static void calc_contrib_bounds(const Emitter *em, Rect *rect)
+static void simstep()
 {
-       int gidx = pos_to_grid(em->pos.x, em->pos.y);
-       int gx = GRID_X(gidx);
-       int gy = GRID_Y(gidx);
-       int maxrange = (int)ceil(CONTRIB_RANGE(em->mass));
+       if(pause) return;
+
+       // move existing particles
+       Particle *p = plist;
+       while(p) {
+               // calculate the field gradient at the particle position
+               int gidx = pos_to_grid(p->pos.x, p->pos.y);
+               Vec2 grad = calc_field_grad(gidx) * GRAV_STR;
+
+               p->vel += grad * SIM_DT;
+               p->pos += p->vel * SIM_DT;
+
+               // if it moved outside of the simulation field, remove it
+               int gx = pos_to_grid_x_noclamp(p->pos.x);
+               int gy = pos_to_grid_y_noclamp(p->pos.y);
+               if(gx < 0 || gx >= GRID_SIZE || gy < 0 || gy >= GRID_SIZE) {
+                       Particle *next = p->next;
+                       grid_part[gidx] = 0;
+                       remove_particle(p);
+                       free_particle(p);
+                       p = next;
+                       continue;
+               }
 
-       int sx = gx - maxrange;
-       int sy = gy - maxrange;
-       int ex = gx + maxrange;
-       int ey = gy + maxrange;
+               // find the grid cell it's moving to
+               int gidx_next = pos_to_grid(p->pos.x, p->pos.y);
+               p->vis_height = 0.0f;//-grid[gidx_next] * field_scale;
 
-       if(ex > GRID_SIZE) ex = GRID_SIZE;
-       if(ey > GRID_SIZE) ey = GRID_SIZE;
+               if(gidx_next == gidx) {
+                       p = p->next;
+                       continue;
+               }
 
-       rect->x = sx < 0 ? 0 : sx;
-       rect->y = sy < 0 ? 0 : sy;
-       rect->width = ex - sx;
-       rect->height = ey - sy;
-}
+               Particle *destp = grid_part[gidx_next];
+               if(destp) {
+                       assert(destp != p);
+                       // another particle at the destination, merge them
+                       destp->vel += p->vel;
+                       destp->mass += p->mass;
+                       destp->radius = MASS_TO_RADIUS(destp->mass);
+
+                       // ... and remove it
+                       grid_part[gidx] = 0;
+                       Particle *next = p->next;
+                       remove_particle(p);
+                       free_particle(p);
+                       p = next;
+               } else {
+                       // destination is empty, go there
+                       if(gidx != gidx_next) {
+                               grid_part[gidx] = 0;
+                               grid_part[gidx_next] = p;
+                       }
 
-static void simstep()
-{
-       // calculate gravitational field - assume field within radius constant: m / r^2
+                       p = p->next;
+               }
+       }
 
-       // first clear the field, and then add contributions
-       memset(grid, 0, sizeof grid);
+       // TODO destroy particles which left the simulation field
 
-       // contribution of emitters
+       // spawn particles
        int num_emitters = emitters.size();
        for(int i=0; i<num_emitters; i++) {
                Emitter *em = emitters[i];
-               Rect cbox;
-               calc_contrib_bounds(em, &cbox);
-               float emradius = MASS_TO_RAD(em->mass);
-
-               float *gptr = grid + cbox.y * GRID_SIZE + cbox.x;
-               Vec2 startpos = grid_to_pos(cbox.x, cbox.y);
+               em->spawn_pending += em->rate * SIM_DT;
+               while(em->spawn_pending >= 1.0f && em->mass > 0.0f) {
+                       Vec2 pvel;      // XXX calc eject velocity
 
-               for(int y=0; y<cbox.height; y++) {
-                       for(int x=0; x<cbox.width; x++) {
-                               Vec2 cellpos = Vec2(startpos.x + (float)x * GRID_DELTA, startpos.y);
+                       float angle = (float)rand() / (float)RAND_MAX * (M_PI * 2.0);
+                       float emradius = MASS_TO_RADIUS(em->mass);
+                       Vec2 ppos = em->pos + Vec2(cos(angle), sin(angle)) * emradius * 1.00001;
 
-                               Vec2 dir = cellpos - em->pos;
-                               float dsq = dot(dir, dir);
-                               float radsq = emradius * emradius;
-                               if(dsq < radsq) {
-                                       dsq = radsq;
-                               }
+                       float pmass = em->chunk > em->mass ? em->mass : em->chunk;
+                       spawn_particle(ppos, pvel, pmass);
+                       em->mass -= pmass;
 
-                               gptr[x] -= em->mass / dsq;
-                       }
+                       em->spawn_pending -= 1.0f;
+               }
+       }
 
-                       startpos.y += GRID_DELTA;
-                       gptr += GRID_SIZE;
+       // remove dead emitters
+       std::vector<Emitter*>::iterator it = emitters.begin();
+       while(it != emitters.end()) {
+               Emitter *em = *it;
+
+               if(em->mass <= 0.0f) {
+                       printf("emitter depleted\n");
+                       it = emitters.erase(it);
+                       delete em;
+               } else {
+                       it++;
                }
        }
 
+       // calculate gravitational field - assume field within radius constant: m / r^2
+       // first clear the field, and then add contributions
+       memset(grid, 0, sizeof grid);
+
+       // contribution of emitters
+       for(int i=0; i<num_emitters; i++) {
+               Emitter *em = emitters[i];
+               Rect cbox;
+               calc_contrib_bounds(em->pos, em->mass, &cbox);
+               float emradius = MASS_TO_RADIUS(em->mass);
+
+               add_influence(em->pos, -em->mass, emradius, cbox);
+       }
+
+       // contribution of particles
+       p = plist;
+       while(p) {
+               Rect cbox;
+               calc_contrib_bounds(p->pos, p->mass, &cbox);
+               add_influence(p->pos, p->mass, p->radius, cbox);
+               p = p->next;
+       }
+
        // update texture
        assert(glGetError() == GL_NO_ERROR);
        grid_tex->bind();
@@ -184,14 +312,15 @@ static void update()
                Emitter *em = new Emitter;
                em->pos = emit_place_pos;
                em->mass = 100;
-               em->rate = 1;
-               em->chunk = 0.05;
+               em->rate = 10;
+               em->chunk = 0.001 * em->mass;
                em->angle = -1;
                em->spread = 0;
+               em->spawn_pending = 0;
                emitters.push_back(em);
 
                Rect cbox;
-               calc_contrib_bounds(em, &cbox);
+               calc_contrib_bounds(em->pos, em->mass, &cbox);
                printf("bounds: %d,%d %dx%d\n", cbox.x, cbox.y, cbox.width, cbox.height);
        }
 
@@ -233,43 +362,48 @@ void GameScreen::draw()
        glMatrixMode(GL_MODELVIEW);
        glLoadMatrixf(view_matrix[0]);
 
-       glPushAttrib(GL_ENABLE_BIT);
-       glDisable(GL_LIGHTING);
-       glDisable(GL_CULL_FACE);
+       // draw gravitational field
+       if(wireframe) {
+               glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
+
+               float amb[] = {0.5, 0.5, 0.5, 1.0};
+               glLightModelfv(GL_LIGHT_MODEL_AMBIENT, amb);
+       } else {
+               float amb[] = {0.01, 0.01, 0.01, 1.0};
+               glLightModelfv(GL_LIGHT_MODEL_AMBIENT, amb);
+       }
 
-       glEnable(GL_TEXTURE_2D);
        bind_texture(gvis_tex, 0);
        bind_texture(grid_tex, 1);
 
        glUseProgram(field_sdr);
+       glPatchParameteri(GL_PATCH_VERTICES, 4);
+       draw_quadmesh(&field_mesh, GL_PATCHES);
 
-       float maxu = FIELD_SIZE / 32.0f;
-       float maxv = FIELD_SIZE / 32.0f;
-       float hsz = FIELD_SIZE * 0.5f;
+       if(wireframe) {
+               glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
+       }
 
-       glBegin(GL_QUADS);
-       glColor3f(1, 1, 1);
+       // draw particles
+       glUseProgram(particle_sdr);
 
-       glMultiTexCoord2f(0, 0, 0);
-       glMultiTexCoord2f(1, 0, 0);
-       glVertex3f(-hsz, 0, -hsz);
+       Particle *p = plist;
+       while(p) {
+               int gidx = pos_to_grid(p->pos.x, p->pos.y);
 
-       glMultiTexCoord2f(0, maxu, 0);
-       glMultiTexCoord2f(1, 1, 0);
-       glVertex3f(hsz, 0, -hsz);
+               glPushMatrix();
+               glTranslatef(p->pos.x, p->vis_height, p->pos.y);
+               glScalef(p->radius, p->radius, p->radius);
 
-       glMultiTexCoord2f(0, maxu, maxv);
-       glMultiTexCoord2f(1, 1, 1);
-       glVertex3f(hsz, 0, hsz);
+               pmesh->draw();
 
-       glMultiTexCoord2f(0, 0, maxv);
-       glMultiTexCoord2f(1, 0, 1);
-       glVertex3f(-hsz, 0, hsz);
-       glEnd();
+               glPopMatrix();
+               p = p->next;
+       }
 
        glUseProgram(0);
 
-       glPopAttrib();
+       assert(glGetError() == GL_NO_ERROR);
 }
 
 void GameScreen::reshape(int x, int y)
@@ -286,6 +420,47 @@ void GameScreen::keyboard(int key, bool pressed)
                        pop_screen();
                        break;
 
+               case ']':
+                       if(tess_level < glcaps.max_tess_level) {
+                               tess_level++;
+                               printf("tessellation level: %d\n", tess_level);
+                               set_uniform_int(field_sdr, "tess_level", tess_level);
+                               glUseProgram(0);
+                       }
+                       break;
+
+               case '[':
+                       if(tess_level > 1) {
+                               tess_level--;
+                               printf("tessellation level: %d\n", tess_level);
+                               set_uniform_int(field_sdr, "tess_level", tess_level);
+                               glUseProgram(0);
+                       }
+                       break;
+
+               case '=':
+                       field_scale += 0.5;
+                       printf("field scale: %f\n", field_scale);
+                       set_uniform_float(field_sdr, "field_scale", field_scale);
+                       break;
+
+               case '-':
+                       field_scale -= 0.5;
+                       if(field_scale < 0.0f) {
+                               field_scale = 0.0f;
+                       }
+                       printf("field scale: %f\n", field_scale);
+                       set_uniform_float(field_sdr, "field_scale", field_scale);
+                       break;
+
+               case 'w':
+                       wireframe = !wireframe;
+                       break;
+
+               case ' ':
+                       pause = !pause;
+                       break;
+
                default:
                        break;
                }
@@ -307,6 +482,13 @@ void GameScreen::mmotion(int x, int y)
        prev_x = x;
        prev_y = y;
 
+       if(game_bnstate(0)) {
+               float pan_speed = pow(cam_dist, 1.5) * 0.00035; // magic
+               Vec2 dir = rotate(Vec2(dx, dy) * pan_speed, deg_to_rad(cam_theta));
+               cam_pos.x -= dir.x;
+               cam_pos.y -= dir.y;
+       }
+
        if(game_bnstate(2)) {
                cam_theta += dx * 0.5;
        }
@@ -319,11 +501,20 @@ void GameScreen::mwheel(int dir, int x, int y)
        if(cam_dist > MAX_CAM_DIST) cam_dist = MAX_CAM_DIST;
 }
 
+static int pos_to_grid_x_noclamp(float x)
+{
+       return ((x / (float)FIELD_SIZE) + 0.5f) * (float)GRID_SIZE;
+}
+
+static int pos_to_grid_y_noclamp(float y)
+{
+       return ((y / (float)FIELD_SIZE) + 0.5f) * (float)GRID_SIZE;
+}
 
 static int pos_to_grid(float x, float y)
 {
-       int gx = ((x / (float)FIELD_SIZE) + 0.5f) * (float)GRID_SIZE;
-       int gy = ((y / (float)FIELD_SIZE) + 0.5f) * (float)GRID_SIZE;
+       int gx = pos_to_grid_x_noclamp(x);
+       int gy = pos_to_grid_y_noclamp(y);
 
        if(gx < 0) gx = 0;
        if(gx >= GRID_SIZE) gx = GRID_SIZE - 1;
@@ -340,3 +531,239 @@ static Vec2 grid_to_pos(int gx, int gy)
 
        return Vec2(x, y);
 }
+
+static void calc_contrib_bounds(const Vec2 &pos, float mass, Rect *rect)
+{
+       int gidx = pos_to_grid(pos.x, pos.y);
+       int gx = GRID_X(gidx);
+       int gy = GRID_Y(gidx);
+       int maxrange = (int)ceil(CONTRIB_RANGE(mass));
+
+       int sx = gx - maxrange;
+       int sy = gy - maxrange;
+       int ex = gx + maxrange;
+       int ey = gy + maxrange;
+
+       if(ex > GRID_SIZE) ex = GRID_SIZE;
+       if(ey > GRID_SIZE) ey = GRID_SIZE;
+
+       rect->x = sx < 0 ? 0 : sx;
+       rect->y = sy < 0 ? 0 : sy;
+       rect->width = ex - sx;
+       rect->height = ey - sy;
+}
+
+static void add_influence(const Vec2 &pos, float mass, float radius, const Rect &cbox)
+{
+       float *gptr = grid + cbox.y * GRID_SIZE + cbox.x;
+       Vec2 startpos = grid_to_pos(cbox.x, cbox.y);
+
+       for(int y=0; y<cbox.height; y++) {
+               for(int x=0; x<cbox.width; x++) {
+                       Vec2 cellpos = Vec2(startpos.x + (float)x * GRID_DELTA, startpos.y);
+
+                       Vec2 dir = cellpos - pos;
+                       float dsq = dot(dir, dir);
+                       float radsq = radius * radius;
+                       if(dsq < radsq) {
+                               dsq = radsq;
+                       }
+
+                       gptr[x] += mass / dsq;
+               }
+
+               startpos.y += GRID_DELTA;
+               gptr += GRID_SIZE;
+       }
+}
+
+static Vec2 calc_field_grad(int gidx)
+{
+       int gx = GRID_X(gidx);
+       int gy = GRID_Y(gidx);
+
+       int nidx = ((gx + 1 >= GRID_SIZE ? gx : gx + 1) << GRID_BITS) | gy;
+       int pidx = ((gx > 0 ? gx - 1 : gx) << GRID_BITS) | gy;
+       float dfdx = grid[nidx] - grid[pidx];
+
+       nidx = (gx << GRID_BITS) | (gy + 1 >= GRID_SIZE ? gy : gy + 1);
+       pidx = (gx << GRID_BITS) | (gy > 0 ? gy - 1 : gy);
+       float dfdy = grid[nidx] - grid[pidx];
+
+       return Vec2(dfdx, dfdy);
+}
+
+
+// ---- quad mesh operations ----
+
+static void destroy_quadmesh(QuadMesh *m)
+{
+       glDeleteBuffers(1, &m->vbo_v);
+       glDeleteBuffers(1, &m->vbo_uv);
+       glDeleteBuffers(1, &m->ibo);
+
+       delete [] m->v;
+       delete [] m->uv;
+       delete [] m->idx;
+}
+
+static void draw_quadmesh(const QuadMesh *m, unsigned int prim)
+{
+       glEnableClientState(GL_VERTEX_ARRAY);
+       glEnableClientState(GL_TEXTURE_COORD_ARRAY);
+
+       glBindBuffer(GL_ARRAY_BUFFER, m->vbo_v);
+       glVertexPointer(3, GL_FLOAT, 0, 0);
+
+       glBindBuffer(GL_ARRAY_BUFFER, m->vbo_uv);
+       glTexCoordPointer(2, GL_FLOAT, 0, 0);
+
+       glBindBuffer(GL_ARRAY_BUFFER, 0);
+
+       glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, m->ibo);
+       glDrawElements(prim, m->num_idx, GL_UNSIGNED_SHORT, 0);
+       glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
+
+       glDisableClientState(GL_VERTEX_ARRAY);
+       glDisableClientState(GL_TEXTURE_COORD_ARRAY);
+}
+
+static void gen_quad_plane(QuadMesh *m, float width, float height, int usub, int vsub)
+{
+       Vec3 *vptr;
+       Vec2 *uvptr;
+       uint16_t *iptr;
+
+       if(usub < 1) usub = 1;
+       if(vsub < 1) vsub = 1;
+
+       int uverts = usub + 1;
+       int vverts = vsub + 1;
+       m->num_verts = uverts * vverts;
+       m->num_quads = usub * vsub;
+       m->num_idx = m->num_quads * 4;
+
+       vptr = m->v = new Vec3[m->num_verts];
+       uvptr = m->uv = new Vec2[m->num_verts];
+       iptr = m->idx = new uint16_t[m->num_idx];
+
+       float du = 1.0f / (float)usub;
+       float dv = 1.0f / (float)vsub;
+
+       float u = 0.0f;
+       for(int i=0; i<uverts; i++) {
+               float x = (u - 0.5f) * width;
+               float v = 0.0f;
+               for(int j=0; j<vverts; j++) {
+                       float y = (v - 0.5f) * height;
+
+                       *vptr++ = Vec3(x, 0, y);
+                       *uvptr++ = Vec2(u, v);
+
+                       if(i < usub && j < vsub) {
+                               int idx = i * vverts + j;
+
+                               *iptr++ = idx;
+                               *iptr++ = idx + 1;
+                               *iptr++ = idx + vverts + 1;
+                               *iptr++ = idx + vverts;
+                       }
+
+                       v += dv;
+               }
+               u += du;
+       }
+
+       glGenBuffers(1, &m->vbo_v);
+       glBindBuffer(GL_ARRAY_BUFFER, m->vbo_v);
+       glBufferData(GL_ARRAY_BUFFER, m->num_verts * 3 * sizeof(float), m->v, GL_STATIC_DRAW);
+
+       glGenBuffers(1, &m->vbo_uv);
+       glBindBuffer(GL_ARRAY_BUFFER, m->vbo_uv);
+       glBufferData(GL_ARRAY_BUFFER, m->num_verts * 2 * sizeof(float), m->uv, GL_STATIC_DRAW);
+
+       glGenBuffers(1, &m->ibo);
+       glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, m->ibo);
+       glBufferData(GL_ELEMENT_ARRAY_BUFFER, m->num_idx * sizeof(uint16_t), m->idx, GL_STATIC_DRAW);
+
+       glBindBuffer(GL_ARRAY_BUFFER, 0);
+       glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
+}
+
+static void spawn_particle(const Vec2 &pos, const Vec2 &vel, float mass)
+{
+       int gidx = pos_to_grid(pos.x, pos.y);
+
+       if(grid_part[gidx]) {
+               // merge with existing
+               Particle *p = grid_part[gidx];
+               p->vel += vel;
+               p->mass += mass;
+               p->radius = MASS_TO_RADIUS(p->mass);
+
+       } else {
+               Particle *p = alloc_particle();
+               p->pos = pos;
+               p->vel = vel;
+               p->mass = mass;
+               p->radius = MASS_TO_RADIUS(mass);
+               grid_part[gidx] = p;
+
+               add_particle(p);
+       }
+}
+
+static void add_particle(Particle *p)
+{
+       if(plist) plist->prev = p;
+
+       p->next = plist;
+       p->prev = 0;
+       plist = p;
+}
+
+static void remove_particle(Particle *p)
+{
+       assert(plist->prev == 0);
+
+       if(plist == p) {
+               assert(p->prev == 0);
+               plist = p->next;
+       }
+       if(p->prev) {
+               p->prev->next = p->next;
+       }
+       if(p->next) {
+               p->next->prev = p->prev;
+       }
+       p->prev = p->next = 0;
+}
+
+// particle allocator
+#define MAX_PFREE_SIZE 256
+static Particle *pfree_list;
+static int pfree_size;
+
+static Particle *alloc_particle()
+{
+       if(pfree_list) {
+               Particle *p = pfree_list;
+               pfree_list = pfree_list->next;
+               pfree_size--;
+               return p;
+       }
+
+       return new Particle;
+}
+
+void free_particle(Particle *p)
+{
+       if(pfree_size < MAX_PFREE_SIZE) {
+               p->next = pfree_list;
+               p->prev = 0;
+               pfree_list = p;
+               pfree_size++;
+       } else {
+               delete p;
+       }
+}