パーティクル同士の衝突と反射（3D）

/**
* collision with many particles in 3D
*
* @author aa_debdeb
* @date 2017/01/05
*/

float MIN_RADIUS = 10;
float MAX_RADIUS = 50;

float e = 1.0;
float k = 2.0;

float rotX = 0;
float rotY = 0;

ArrayList<Particle> particles;

void setup(){
  size(500, 500, P3D);
  noStroke();
  frameRate(30);
  particles = new ArrayList<Particle>();
  while(particles.size() < 30){
    float radius = random(MIN_RADIUS, MAX_RADIUS);
    float locSize = 200.0;
    float locAng1 = random(PI);
    float locAng2 = random(TWO_PI);
    PVector loc = new PVector(locSize * sin(locAng1) * cos(locAng2), locSize * cos(locAng1), locSize * sin(locAng1) * sin(locAng2));
    boolean isOverlapping = false;
    for(Particle p: particles){
      if(PVector.dist(loc, p.loc) <= radius + p.radius){
        isOverlapping = true;
        break;
      }
    }
    if(!isOverlapping){
      float velSize = 5.0;
      float velAng1 = random(PI);
      float velAng2 = random(TWO_PI);
      PVector vel = new PVector(velSize * sin(velAng1) * cos(velAng2), velSize * cos(velAng1), velSize * sin(velAng1) * sin(velAng2));
      particles.add(new Particle(loc, vel, radius));
    }
  }  
}

void draw(){
  background(230);
  translate(width / 2, height / 2);
  lights();
  for(Particle p: particles){
    p.render();
    p.move();
  }
  
  for(Particle p1: particles){
    for(Particle p2: particles){
      if(p1 == p2){continue;}
      float d = PVector.dist(p1.loc, p2.loc);
      if(d <= p1.radius + p2.radius){
        PVector p12 = PVector.sub(p2.loc, p1.loc);
        PVector n = PVector.div(p12, p12.mag());
        PVector v12 = PVector.sub(p2.vel, p1.vel);
        PVector vn1 = PVector.mult(n, PVector.dot(p1.vel, n));
        PVector vt1 = PVector.sub(p1.vel, vn1);
        PVector t = PVector.div(vt1, vt1.mag());
        float spring = -k * (p1.radius + p2.radius - d);
        float j = (1 + e) * (p1.mass * p2.mass / (p1.mass + p2.mass)) * PVector.dot(v12, n);
        PVector impulse = PVector.mult(n, j + spring); 
        p1.nvel.add(impulse);
      }
    }
  }
  
  for(Particle p: particles){
    p.updateVel();
  }
  
}

class Particle{
  
  PVector loc, vel, nvel;
  float radius, mass;
  color c;
  
  Particle(PVector loc, PVector vel, float radius){
    this.loc = loc;
    this.vel = vel;
    this.nvel = new PVector(vel.x, vel.y, vel.z);
    this.radius = radius;
    this.mass = 1;
    c = color(lerpColor(color(255, 140, 0), color(255, 0, 170), pow((radius - MIN_RADIUS) / (MAX_RADIUS - MIN_RADIUS), 3)), 200);
  }
  
  void move(){
    PVector center = new PVector(0, 0, 0);
    PVector acc = PVector.sub(center, loc);
    acc.limit(1.0);
    vel.add(acc);
    vel.limit(7.0);
    nvel = new PVector(vel.x, vel.y, vel.z);
    loc.add(vel);
  }
  
  void render(){
    fill(c);
    pushMatrix();
    translate(loc.x, loc.y, loc.z);
    sphere(radius);
    popMatrix();
  }
  
  void updateVel(){
    vel = nvel;
    nvel = new PVector(vel.x, vel.y, vel.z);
  }
}