+static void calc_contrib_bounds(const Emitter *em, Rect *rect)
+{
+ 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));
+
+ 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 simstep()
+{
+ // 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
+ 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);
+
+ 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 - em->pos;
+ float dsq = dot(dir, dir);
+ float radsq = emradius * emradius;
+ if(dsq < radsq) {
+ dsq = radsq;
+ }
+
+ gptr[x] -= em->mass / dsq;
+ }
+
+ startpos.y += GRID_DELTA;
+ gptr += GRID_SIZE;
+ }
+ }
+
+ // update texture
+ assert(glGetError() == GL_NO_ERROR);
+ grid_tex->bind();
+ glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, GRID_SIZE, GRID_SIZE, GL_LUMINANCE,
+ GL_FLOAT, grid);
+ assert(glGetError() == GL_NO_ERROR);
+}
+
+static void update()
+{
+ if(emit_place_pending) {
+ emit_place_pending = false;
+ Emitter *em = new Emitter;
+ em->pos = emit_place_pos;
+ em->mass = 100;
+ em->rate = 1;
+ em->chunk = 0.05;
+ em->angle = -1;
+ em->spread = 0;
+ emitters.push_back(em);
+
+ Rect cbox;
+ calc_contrib_bounds(em, &cbox);
+ printf("bounds: %d,%d %dx%d\n", cbox.x, cbox.y, cbox.width, cbox.height);
+ }
+
+ // update simulation
+ static float interval;
+ interval += frame_dt;
+ if(interval >= SIM_DT) {
+ interval -= SIM_DT;
+ simstep();
+ assert(glGetError() == GL_NO_ERROR);
+ }
+
+ // update projection matrix
+ proj_matrix.perspective(deg_to_rad(60.0f), win_aspect, 0.5, 5000.0);
+
+ // update view matrix
+ Vec3 targ;
+ if(targ_pos) {
+ targ.x = targ_pos->x;
+ targ.z = targ_pos->y;
+ }
+
+ float theta = -deg_to_rad(cam_theta);
+ Vec3 camdir = Vec3(sin(theta) * cam_dist, pow(cam_dist * 0.1, 2.0) + 0.5, cos(theta) * cam_dist);
+ Vec3 campos = targ + camdir;
+
+ view_matrix.inv_lookat(campos, targ, Vec3(0, 1, 0));
+}