斥力を働かせ合う円の中のパーティクル

/**
* particles with reputation in circle
*
* @author aa_debdeb
* @date 2016/06/16
*/

float G = 1000.0;
float MAX_VEL = 5.0;
float MAX_ACC = 1.0;
float L_RADIOUS = 200;
float S_RADIOUS = 10;
ArrayList<Particle> particles;

void setup(){
  size(500, 500);
  particles = new ArrayList<Particle>();
}

void draw(){
  background(200);
  translate(width / 2, height / 2);
  for(Particle p1: particles){
    for(Particle p2: particles){
      if(p1 != p2){
        p1.addForce(p2);
      }
    }
  }
  noStroke();
  fill(60);
  ellipse(0, 0, (L_RADIOUS + S_RADIOUS) * 2, (L_RADIOUS + S_RADIOUS) * 2);
  for(Particle p: particles){
    p.update();
    p.display();
  }
}

void mousePressed(){
  if(sqrt(sq(mouseX - width / 2) + sq(mouseY - height / 2)) < L_RADIOUS){
    particles.add(new Particle(mouseX - width / 2, mouseY - height / 2));  
  }
}

class Particle{

  PVector pos, vel, acc;
  
  Particle(float x, float y){
    pos = new PVector(x, y);
    vel = new PVector(0, 0);
    acc = new PVector(0, 0);
  }
  
  void addForce(Particle p){
    float d = sq(pos.dist(p.pos));
    PVector force = PVector.sub(pos, p.pos);
    force.normalize();
    force.mult(G / d);
    acc.add(force);  
  }
  
  void display(){
    noStroke();
    fill(128, 255, 0, 200);
    ellipse(pos.x, pos.y, S_RADIOUS * 2, S_RADIOUS * 2);
    fill(255, 200);
    ellipse(pos.x, pos.y, 5, 5);    
  }
  
  void update(){
    acc.limit(MAX_ACC);
    vel.add(acc);
    vel.limit(MAX_VEL);
    acc = new PVector(0, 0);
    pos.add(vel);
    if(pos.mag() > L_RADIOUS){
      float posAng = atan2(pos.y, pos.x);
      float velAng = atan2(vel.y, vel.x);
      float velSize = vel.mag();
      vel.x = velSize * cos(2 * posAng - velAng - PI);
      vel.y = velSize * sin(2 * posAng - velAng - PI);
      pos.add(vel);
    }
  }
}