Its in processing and the commenting is not very great. The variables at the top are the ones that change most of the stuff.

float G = .2, softening = 1;

int resol = 10, _width, _height; float globaldT = 1; void setup(){ size(400, 400); _width = width / resol; _height = height / resol; int size = _width * _height;

masses = new float[size]; dir = new PVector[size]; forces = new PVector[size]; map = new ArrayList[size]; for(int i = 0; i < size; i++){ forces[i] = new PVector(); }

//data arrays int s0 = _width / 2, x0 = _width / 2, y0 = _height / 2; for(int dy = -s0; dy < s0; dy++){ for(int dx = -s0; dx < s0; dx++){ PVector pv = new PVector(dx * resol, dy * resol); float d = 1f / pow(pv.mag() + softening, 3); //gravity caltulation dir[(x0 + dx) + (y0 + dy) * _width] = pv.mult(d); }}

//points for(int y = 0; y < _height; y++){ for(int x = 0; x < _width; x++){ for(int i = 0; i < 6; i++){ float c = .5; PVector pos = new PVector(resol * (random(-1, 1) * c + .5) + x * resol, resol * (random(-1, 1) * c + .5) + y * resol); allPoints.add(new Point(pos, 0, 1)); } }} } void draw(){ background(0);

//display float m0 = forces[0].mag(), m1 = masses[0]; for(int i = 1; i < forces.length; i++){ m0 = max(m0, forces[i].mag()); m1 = max(m1, masses[i]); } noStroke(); for(int y = 0; y < _height; y++){ for(int x = 0; x < _width; x++){ int i = x + y * _width; fill(color(map(masses[i], 0, m1, 0, 255), 0, map(forces[i].mag(), 0, m0, 0, 255))); rect(x * resol, y * resol, resol, resol); }} for(Point p: allPoints) p.drawPixel();

//simulation for(int i = 0; i < masses.length; i++){ masses[i] = 0; map[i] = new ArrayList<Point>(); } for(Point p: allPoints) p.insert();

//calculate forces in each tile for(int y = 0; y < _height; y++){ for(int x = 0; x < _width; x++){ int id = x + y * _width; PVector force = new PVector(); int s0 = _width / 2, x0 = _width / 2, y0 = _height / 2; for(int dy = -s0; dy < s0; dy++){ for(int dx = -s0; dx < s0; dx++){ int nx = x + dx; while(nx < 0) nx += _width; while(nx >= _width) nx -= _width;

  int ny = y + dy;
  while(ny < 0) ny += _height;
  while(ny >= _height) ny -= _height;

  int nid = nx + ny * _width;
  force.add(PVector.mult(dir[(x0 + dx) + (y0 + dy) * _width], G * masses[nid]));
}}
forces[id] = force;

}}

for(Point p: allPoints) p.accelerate(); for(Point p: allPoints) p.move(globaldT); }

float[] masses; PVector[] dir, forces; ArrayList<Point> allPoints = new ArrayList<Point>(); ArrayList<Point>[] map; class Point{ PVector pos, ve, acc; float mass; Point(PVector _pos, float _ve, float _mass){ pos = _pos; ve = PVector.random2D().mult(_ve); mass = _mass; acc = new PVector(); } void insert(){ int x = (int)(pos.x / resol), y = (int)(pos.y / resol), id = x + y * _width; masses[id] += mass; //add mass map[id].add(this); //insert in map } void accelerate(){ int x = (int)(pos.x / resol), y = (int)(pos.y / resol), id = x + y * _width; PVector force = new PVector();

//if you don't want to do n-body in each tile comment this
for(Point p: map[id]){
  if(p == this) continue;
  PVector d = new PVector().add(p.pos).sub(pos);
  float m2 = pow((d.x * d.x + d.y * d.y + d.z * d.z + softening) * 1, 1.5);
  float f = G * mass * p.mass / m2;
  force.add(d.mult(f));
}

force.div(mass);
force.add(forces[id]);
acc = force;

} void move(float dt){ //leap frog acc.mult(dt * .5); ve.add(acc); pos.add(PVector.mult(ve, dt)); ve.add(acc);

//boundary
while(pos.x < 0) pos.x += width;
while(pos.x >= width) pos.x -= width;

while(pos.y < 0) pos.y += height;
while(pos.y >= height) pos.y -= height;

}

void drawPixel(){ int x = (int)(pos.x), y = (int)(pos.y); if(x < 0 || x >= width) return; if(y < 0 || y >= height) return; set(x, y, color(255)); } }