BImage bg;
Planet[] p;

float rx, ry;
int mdownX, mdownY;

double galaxyCX, galaxyCY, galaxyCZ;

float G = .012;
int nPlanets = 400;
double initWeightSeed = 4;
double initRandSpeed = 2;
double initGalaxySpeed = 1;
int initGalaxyWidth = 40;
int initGalaxyMin = 310;
int initGalaxyMax = 490;

void setup()
{
  size(800, 800);
  bg = loadImage("sombrerogalaxy.jpg"); 
  background(bg);
  p = new Planet[nPlanets];
  initialize(nPlanets);
}

void initialize(int n)
{
  nPlanets = n;

  double initX, initY, initZ, initdX, initdY, initdZ, initW;
  int initc1, initc2;

  for(int i=0; i<nPlanets; i++)
  {
    initX = random(initGalaxyWidth);
    initY = random(initGalaxyWidth);
    initZ = random(initGalaxyWidth/2) - initGalaxyWidth/4;
    initdX = random(int(initRandSpeed)) - initRandSpeed/2;
    initdY = random(int(initRandSpeed)) - initRandSpeed/2;
    initdZ = random(int(initRandSpeed/2)) - initRandSpeed/4;
    initW = random(int(initWeightSeed))+1;
    initc1 = 255;
    initc2 = 0;

    int clusters = nPlanets/4;
    if(i<clusters)
    {
      initX += initGalaxyMin;
      initY += initGalaxyMin;
      initdX -= initGalaxySpeed;
      initdY += initGalaxySpeed/2;
    }
    else if(i<clusters*2)
    {
      initX += initGalaxyMin;
      initY += initGalaxyMax;
      initdX += initGalaxySpeed;
      initdY += initGalaxySpeed;
    }
    else if(i<clusters*3)
    {
      initX += initGalaxyMax;
      initY += initGalaxyMin;
      initdX -= initGalaxySpeed;
      initdY -= initGalaxySpeed;
    }
    else if(i<clusters*4)
    {
      initX += initGalaxyMax;
      initY += initGalaxyMax;
      initdX += initGalaxySpeed;
      initdY -= initGalaxySpeed/2;
    }

    p[i] = new Planet(initX, initY, initZ, initdX, initdY, initdZ, initW, initc1, initc2);
  }
}

void loop()
{
  translate(width/2, height/2, 0);
  rotateY(ry);
  rotateX(rx);
  translate(-width/2, -height/2, 0);

  if(!keyPressed)
    background(bg);

  noFill();
  stroke(100);
  line(300,400,0, 500,400,0);
  line(400,300,0, 400,500,0);
  line(400,400,-100, 400,400,100);
  ellipse(325,325,150,150);

  updatePlanets();
}

void mouseDragged()
{
  ry += (mouseX - pmouseX)*.01;
  rx -= (mouseY - pmouseY)*.01;
}

void keyPressed()
{
  if(key == '1')
    initialize(100);
  else if(key == '2')
    initialize(200);
  else if(key == '3')
    initialize(300);
  else if(key == '4')
    initialize(400);
}

class Planet
{
  double X, Y, Z, dX, dY, dZ, W, Xprev, Yprev, Zprev;
  int C1, C2;

  Planet(double x, double y, double z, double dx, double dy, double dz, double w, int c1, int c2)
  {
    X = x;
    Y = y;
    Z = z;
    dX = dx;
    dY = dy;
    dZ = dz;
    W = w;
    C1 = c1;
    C2 = c2;
  }

  void draw()
  {
    if(keyPressed)
    {
      stroke(C1, 0, 255-C1);
      point(int(Xprev), int(Yprev), int(Zprev));
    }
    stroke(255);
    point(int(X), int(Y), int(Z));
  }
}

void updatePlanets()
{
  double a, ax, ay, az;
  double r, rx, ry, rz, rxn, ryn, rzn;

  galaxyCX = 0;
  galaxyCY = 0;
  galaxyCZ = 0;

  for(int i=0; i<nPlanets; i++)
  {
    for(int c=0; c<nPlanets; c++)
    {
      if(c != i)
      {
        rx = ry = rz = 0;
        rx = p[c].X - p[i].X;
        ry = p[c].Y - p[i].Y;
        rz = p[c].Z - p[i].Z;
        r = rx*rx + ry*ry + rz*rz;

        rxn = rx / (Math.abs(rx) + Math.abs(ry) + Math.abs(rz));
        ryn = ry / (Math.abs(rx) + Math.abs(ry) + Math.abs(rz));
        rzn = rz / (Math.abs(rx) + Math.abs(ry) + Math.abs(rz));
        a = G*p[i].W/Math.sqrt(r);
        ax = a*rxn;
        ay = a*ryn;
        az = a*rzn;

        p[i].dX += ax;
        p[i].dY += ay;
        p[i].dZ += az;
      }
    }
    p[i].C1 = int(Math.abs(p[i].dX + p[i].dY + p[i].dZ)*50);
    p[i].Xprev = p[i].X;
    p[i].Yprev = p[i].Y;
    p[i].Zprev = p[i].Z;
    p[i].X += p[i].dX;
    p[i].Y += p[i].dY;
    p[i].Z += p[i].dZ;
    p[i].draw();

    galaxyCX += p[i].X-400;
    galaxyCY += p[i].Y-400;
    galaxyCZ += p[i].Z;
  }

  galaxyCX = int(galaxyCX / nPlanets + 400);
  galaxyCY = int(galaxyCY / nPlanets + 400);
  galaxyCZ = int(galaxyCZ / nPlanets);

  if(keyPressed && (key == 'c'))
  {
    for(int i=0; i<nPlanets; i++)
    {
      p[i].X -= galaxyCX - 400;
      p[i].Y -= galaxyCY - 400;
      p[i].Z -= galaxyCZ;
    }
  }

  stroke(0, 255, 0);
  point(int(galaxyCX), int(galaxyCY), int(galaxyCZ));
}