Java/Processing calculations returning "NaN" and "Infinity"

Question

I am trying to make simulation of the three-body problem in Processing/Java. I have created a Planet object and defined several functions to calculate the attraction, the forces, and the speeds of the Planets. The program runs without any errors. The hitch is that the calculation are returning Infinity and NaN outputs. I haven't found any divide by zero errors. Why is it doing this?

float simulationTime = 0.01;
float earthMass = pow(10,24)*5.972;

Planet[] bodies = {new Planet(100,200, 0, 0,earthMass/2),
                   new Planet(400,500, 0, 0,earthMass/2),
                   new Planet(800,200, 0, 0,earthMass/2)};
int size = 500;
float G = 6.674*pow(10, -11);
float gravAttract(float mass1, float mass2, float x1, float y1, float x2, float y2) {
  float force = G*mass1*mass2/pow(dist1(x1, y1, x2, y2), 2);
  //println(pow(dist1(x1, y1, x2, y2), 2));
  //println(force);
  return force;
}
float attract(Planet one,Planet two){
  float force = G*one.mass*two.mass/pow(dist1(one.xpos,one.ypos,two.xpos,two.ypos),2);
  println(one.xpos,two.xpos,one.xspeed,two.xspeed);
  return force;
}
float[] forceXY(int body){
  float[] xy = {0,0};
  for (int ii = 0; ii<bodies.length; ii = ii+1){
    if (ii == body){
      continue;
    }
    else{
      xy[0] = xy[0]+attract(bodies[body],bodies[ii])*
      (bodies[ii].xpos-bodies[body].xpos)/dist1(bodies[body].xpos,bodies[body].ypos,bodies[ii].xpos,bodies[ii].ypos);
      xy[1] = xy[1]+attract(bodies[body],bodies[ii])*
      (bodies[ii].ypos-bodies[body].ypos)/dist1(bodies[body].xpos,bodies[body].ypos,bodies[ii].xpos,bodies[ii].ypos);
      println(xy);
    }
  }
  return xy;
}
float dist1(float x1,float y1,float x2,float y2){
  float x = dist(x1,y1,x2,y2);
  x = x*149.6*pow(10,7);
  //println(x1,y1,x2,y2);
  return x;
}
class Planet {
  float xpos;
  float ypos;
  float xspeed;
  float yspeed;
  float mass;

  Planet(float a, float b, float c, float d, float e) {
    xpos = a;
    ypos = b;
    xspeed = c;
    yspeed = d;
    mass = e;
  }
  void update(float xforce, float yforce) {
    //println("xy",xpos,ypos);
    float xa = xforce/mass;
    float ya = yforce/mass;
    xspeed = xspeed + xa*simulationTime;
    yspeed = yspeed + ya*simulationTime;
    xpos = xpos + xspeed*simulationTime;
    ypos = ypos + yspeed*simulationTime;
  }
}
void setup() {
  size(1000, 1000);
  frameRate(1);
}
void draw() {
  background(0);
  
  for (int ii = 0; ii < bodies.length; ii = ii+1){
    float[] asdf = new float[2];
    asdf = forceXY(ii);
    circle(bodies[ii].xpos,bodies[ii].ypos,10);
    bodies[ii].update(asdf[0],asdf[1]);
  }
}

Answer 1

I see you're trying to do the calculations using the real physical values. The problem with that is that the astronomical masses, distances and forces are astronomical, not just literally but also figuratively. The numbers are so big that the results of the intermediate calculations cannot be represented using the limited range of the float data type.

The solution is to change the units. For example, instead of using meters, kilograms, and seconds, measure the distance in "astronomical units" (distance between Earth and Sun), mass in multiples of Earth masses, and time in years. Or use the standard "astronomical system of units". Most importantly, this changes the value of the gravitational constant G.

If you're not too attached to making a simulation that lets you predict the real positions of the planets, it's much easier to just "wing it" and set G=1, set the masses to small numbers, and play around with the velocities. You'll see it's very easy to come up with a realistic simulation.