I'm attempting to simulate magnetic dipoles (bar magnets) by using two magnetically charged particles and constraining them to a center point. It isn't quite working. Two bar magnets work fine until, for some reason, one of them "collapses" into one point.

Code:

ArrayList dipoles;

void setup () {
  size(600, 600);
  background(255, 255, 255);
  noFill();
  
  dipoles = new ArrayList();
}

void draw () {
  background(255, 255, 255);
  
  for (int i=0; i<dipoles.size(); i++) {
    Dipole d = (Dipole) dipoles.get(i);
    d.update(dipoles);
  }
}

void mouseClicked () {
  Dipole d = new Dipole(new PVector(mouseX, mouseY));
  dipoles.add(d);
} 

Code:

class Dipole {

  Node[] ends;
  public PVector position;

  Dipole (PVector p) {
    position = p;

    ends = new Node[2];
    ends[0] = new Node( new PVector(position.x, position.y + 20), 1.0 );
    ends[1] = new Node( new PVector(position.x, position.y - 20), -1.0 );
  }

  Node[] getMembers () {
    return ends;
  }

  void update (ArrayList peers) {
    for (int i=0; i<peers.size(); i++) {
	Dipole d = (Dipole) peers.get(i);
	if (d == this) continue;

	ends[0].update(d.getMembers());
	ends[1].update(d.getMembers());
    }

    ends[0].move();
    ends[1].move();

    constrainNodes();

    render();
  }

  void constrainNodes () {
    PVector distance;
    float delta;

    distance = PVector.sub(ends[0].position, ends[1].position);
    delta = (distance.mag() - 40) / distance.mag();
    distance.mult(0.5 * delta);
    ends[0].position.add(distance);
    ends[1].position.add(distance);

    distance = PVector.sub(position, ends[0].position);
    delta = (distance.mag() - 20) / distance.mag();
    distance.mult(delta);
    ends[0].position.add(distance);

    distance = PVector.sub(position, ends[1].position);
    delta = (distance.mag() - 20) / distance.mag();
    distance.mult(delta);
    ends[1].position.add(distance);
  }

  void render () {
    stroke(60, 60, 60);
    line(ends[0].position.x, ends[0].position.y, ends[1].position.x, ends[1].position.y);
    
    ends[0].render();
    ends[1].render();
  }
} 

Code:

class Node {

  public PVector position;
  PVector velocity;
  PVector acceleration;

  public float charge;
  float mass;
  float damp;

  ////// Constructors //////
  Node (PVector p, float q) {
    position = p;
    velocity = new PVector( 0, 0);
    acceleration = new PVector( 0, 0);

    charge = q;
    mass = 1.0;
    damp = 0.9;
  }
  //////////////////////////

  void update (Node[] peers) {

    //Magnetic Replusion/Attraction Algorithm
    PVector radius;
    PVector direction;
    PVector force;

    for (int i=0; i<peers.length; i++) {
	Node peer = peers[i];

	radius = PVector.sub(position, peer.position);
	direction = PVector.sub(position, peer.position);
	direction.normalize();
	force = PVector.mult( direction, ( charge * peer.charge / (radius.mag() * radius.mag()) ) + ( abs(charge) * abs(peer.charge) / pow(radius.mag(), 12) ) );
	force.mult( 60.0 / mass);
	acceleration = force;

	if (radius.mag() > 0.01)
	  velocity.add(acceleration);
    }

    velocity.mult(damp);
  }

  void move () {
    position.add(velocity);
  }

  void render () {
    stroke( map(charge, -1.0, 1.0, 0, 255), 0, map(charge, -1.0, 1.0, 255, 0) );
    rect(position.x-5, position.y-5, 10, 10);
  }
} 

I guess this might be stating the obvious but the problem looks like it's in the constrainNodes method.  I did a println on distance and whilst it's presumably intended to maintain a constant distance between two nodes, it's not working.  In fact in some cases it was showing NaN (in fact this causes the 'magnet' to disappear).

Code:

  void constrainNodes () {
    PVector distance;
    float delta;

    distance = PVector.sub(ends[0].position, ends[1].position);
    delta = (distance.mag() - 40) / distance.mag();
    distance.mult(0.5 * delta);
    println(distance);
    ...