EDIT: code now hosted on sourceforge: https://sourceforge.net/projects/octree ... particles/
============
i'm going to post in a new thread here as i should have done before. not sure i can move my other posts into here, but i wouldn't complain if an admin did
Note this is an O(N) time n-body simulator. That means that the simulation time only grows linearly with the number of particles. That means you can scale it to HUGEMONGOUS particle counts. and it fits well onto non-uniform memory architectures (NUMA), so hardware-wise, with a few modifications you can put it on massive computer clusters without being bottlenecked by i/o - the i/o pattern stays grid-like, and the i/o demand per channel grows very slowly to not-at-all with the simulation size.
Code for the Particle class:
Code: Select all
public class Particle {
static final float ATOMIC_CHARGE = 1;
static final float ELECTRON_MASS = 1;
static final float PROTON_MASS = 1;
static final float RC2 = 1; //reciprocal speed of light squared?
float x,y,z;
float dx,dy,dz;
public float ddx,ddy,ddz;
public float oddx,oddy,oddz;
public int element_number = 0;
public float mass = 0;
public float charge = 0;
public static float[] charge_density = new float[]{-ATOMIC_CHARGE/ELECTRON_MASS,ATOMIC_CHARGE/PROTON_MASS};
boolean rigid = false;
public Particle(int element_number2, float x2, float y2, float z2) {
element_number = element_number2;
x = x2;
y = y2;
z = z2;
}
public Particle(int element_number2, float x2, float y2, float z2, float mass, boolean rigid) {
element_number = element_number2;
x = x2;
y = y2;
z = z2;
this.mass = mass;
this.rigid = rigid;
this.charge = mass*charge_density[element_number];
}
public float[] getNetForceOnParticle(Particle p) {
float deltax = x-p.x;
float deltay = y-p.y;
float deltaz = z-p.z;
float rds = (float)1.0/(deltax*deltax + deltay*deltay + deltaz*deltaz); //reciprocal distance squared
float reciprocal_norm = (float)Math.sqrt(rds);
float electrostatic_force = charge*p.charge*rds;
//electric force vector
float ex = electrostatic_force*deltax*reciprocal_norm;
float ey = electrostatic_force*deltay*reciprocal_norm;
float ez = electrostatic_force*deltaz*reciprocal_norm;
//magnetic force vector
float deltadx = dx-p.dx;
float deltady = dy-p.dy;
float deltadz = dz-p.dz;
float bx = RC2*(deltady*ez + deltadz*ey);
float by = RC2*(deltadz*ex + deltadx*ez);
float bz = RC2*(deltadx*ey + deltady*ex);
float mx = -(deltady*bz + deltadz*by);
float my = -(deltadz*bx + deltadx*bz);
float mz = -(deltadx*by + deltady*bx);
return new float[]{
ex+mx,
ey+my,
ez+mz,
};
}
public void combine(Particle p) {
float rsm = (float)1.0/(mass+p.mass); //reciprocal sum of masses;
x = rsm*(x*mass+p.x*p.mass);
y = rsm*(y*mass+p.y*p.mass);
z = rsm*(z*mass+p.z*p.mass);
dx = rsm*(dx*mass+p.dx*p.mass);
dy = rsm*(dy*mass+p.dy*p.mass);
dz = rsm*(dz*mass+p.dz*p.mass);
mass += p.mass;
charge += p.charge;
}
public void separate(Particle p) {
float rsm = (float)1.0/(mass-p.mass); //reciprocal sum of masses;
x = rsm*(x*mass-p.x*p.mass);
y = rsm*(y*mass-p.y*p.mass);
z = rsm*(z*mass-p.z*p.mass);
dx = rsm*(dx*mass-p.dx*p.mass);
dy = rsm*(dy*mass-p.dy*p.mass);
dz = rsm*(dz*mass-p.dz*p.mass);
mass -= p.mass;
charge -= p.charge;
}
public void applyForcesToParticle(Particle p) {
float[] forces = getNetForceOnParticle(p);
p.ddx += forces[0];
p.ddy += forces[1];
p.ddz += forces[2];
}
public void integrateTimeStep(float dt) {
if( rigid)
return;
//do initial setup calculations
float dtt = dt*dt/2;
float dttt = dt*dtt/3;
float rdt = (float)1.0/dt;
ddx /= mass;
ddy /= mass;
ddz /= mass;
float dddx = (ddx-oddx)*rdt;
float dddy = (ddy-oddy)*rdt;
float dddz = (ddz-oddz)*rdt;
//adjustLastTimeStep
x += dddx*dttt;
y += dddy*dttt;
z += dddz*dttt;
dx += dddx*dtt;
dy += dddy*dtt;
dz += dddz*dtt;
//doThisTimeStep
x += dx*dt + ddx*dtt;
y += dy*dt + ddy*dtt;
z += dz*dt + ddz*dtt;
dx += ddx;
dy += ddy;
dz += ddz;
//store force for next time step
oddx = ddx;
oddy = ddy;
oddz = ddz;
ddx = 0;
ddy = 0;
ddz = 0;
}
}
Code: Select all
import java.util.*;
public class Octet {
static float threshold = (float)0.1;
float x,y,z; //the pseudo-particles store the inertias (dx,dy,dz)
float delta; //radius
float dv; //volume
Octet parent = null;
Octet[][][] children;
Vector<Particle> particles;
Vector<Particle> pseudo_particles;
float num_particles = 0;
public static void main() {
Octet rootOctet = new Octet(0,0,0,128,(float)0.125,null);
Vector<Particle> particles = new Vector<Particle>();
//TODO: generate or load a particle list here
iterate( particles, rootOctet, (float)0.01);
}
public Octet(float x, float y, float z, float delta, float min_delta, Octet parent) {
this.x = x;
this.y = y;
this.z = z;
this.delta = delta;
this.dv = 8*delta*delta*delta; //since delta is only half the length of a side, we have to multiply by 2*2*2.
this.parent = parent;
float half_delta = delta/2;
if( delta > min_delta) {
children = new Octet[2][2][2];
for( int dz = 0; dz < 2; dz++) {
for( int dy = 0; dy < 2; dy++) {
for( int dx = 0; dx < 2; dx++) {
children[dz][dy][dx] = new Octet(
x-half_delta+delta*dx,
y-half_delta+delta*dy,
z-half_delta+delta*dz,
half_delta,
min_delta,
this
);
}
}
}
} else {
particles = new Vector<Particle>();
}
}
public void addParticleDown(Particle p) {
num_particles++;
if( particles == null) {
children[p.z > z ? 1 : 0][p.y > y ? 1 : 0][p.x > x ? 1 : 0].addParticleDown(p);
} else {
particles.add(p);
}
//now add it at this level
for( Particle pp : pseudo_particles) {
if( pp.element_number == p.element_number) {
pp.combine(p);
return;
}
}
Particle pp = new Particle(p.element_number,x,y,z);
pp.combine(p);
pseudo_particles.add(p);
}
public void removeParticleDown(Particle p) {
num_particles--;
for( Particle pp : pseudo_particles) {
if( pp.element_number == p.element_number) {
pp.separate(p);
if( pp.mass == 0)
pseudo_particles.remove(pp);
break;
}
}
if( particles == null) {
children[p.z > z ? 1 : 0][p.y > y ? 1 : 0][p.x > x ? 1 : 0].removeParticleDown(p);
} else {
particles.remove(p);
}
}
public static void iterate(Vector<Particle> particles, Octet rootOctet, float deltat) {
for( Particle p : particles) {
rootOctet.addParticleDown(p);
p.ddx = 0;
p.ddy = 0;
p.ddz = 0;
}
for( Particle p : particles) {
rootOctet.applyForcesToParticle(p);
}
for( Particle p : particles) {
rootOctet.removeParticleDown(p);
}
for( Particle p : particles) {
p.integrateTimeStep(deltat);
}
}
public void applyForcesToParticle(Particle p) {
//if at bottom, just apply forces
if( particles != null) {
for( Particle ps : particles) {
ps.applyForcesToParticle(p);
}
return;
}
if( getNetPsuedoParticleForceOnParticle(p) > threshold) {
for( int z = 0; z < 2; z++) {
for( int y = 0; y < 2; y++) {
for( int x = 0; x < 2; x++) {
children[z][y][x].applyForcesToParticle(p);
}
}
}
return;
} else {
for( Particle ps : pseudo_particles) {
ps.applyForcesToParticle(p);
}
return;
}
}
float getDistanceSquared(Particle p) {
float deltax = p.x-x;
float deltay = p.y-y;
float deltaz = p.z-z;
return deltax*deltax + deltay*deltay + deltaz*deltaz;
}
//TODO: give this function a better implementation.
public float getNetPsuedoParticleForceOnParticle(Particle p) {
return num_particles/getDistanceSquared(p);
}
}