+ if(pause) return;
+
+ // distribute forces
+ Particle *p = plist;
+ while(p) {
+ int num_emitters = emitters.size();
+ for(int i=0; i<num_emitters; i++) {
+ Emitter *em = emitters[i];
+ Vec2 dir = p->pos - em->pos;
+ p->vel += dir * em->mass * GRAV_STR * 0.01 / dot(dir, dir) * SIM_DT;
+ }
+
+ Particle *q = plist;
+ while(q) {
+ if(p != q) {
+ Vec2 dir = q->pos - p->pos;
+ float accel = GRAV_STR * q->mass / dot(dir, dir);
+ p->vel += dir * accel * SIM_DT;
+ }
+ q = q->next;
+ }
+ p = p->next;
+ }
+
+ // move existing particles
+ 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;
+ }
+
+ // 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(gidx_next == gidx) {
+ p = p->next;
+ continue;
+ }
+
+ 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;
+ }
+
+ p = p->next;
+ }
+ }
+
+ // TODO destroy particles which left the simulation field
+