class particleSystem {
  int numOfParticles;
  int numOfSprings;

  Particle   particles      [];
  Particle   ptemp          []; 
  Spring     s              [];
  Derivative deriv          [];
  float      gravitational = 0.0f;
  float      drag          = 0.1f;

  int ptCount;

  particleSystem(int numOfParticles_, int numOfSprings_ ) {
    numOfParticles = numOfParticles_;
    numOfSprings   = numOfSprings_;

    particles      = new Particle[numOfParticles];
    ptemp          = new Particle[numOfParticles];
    deriv          = new Derivative[numOfParticles];
    s              = new Spring[numOfSprings];
    

    for (int i = 0; i < numOfParticles; i++) {
      particles[i]               = new Particle(); 
      particles[i].mass          = 10;
      particles[i].fixed         = false;
      
      Vector3D pos = new  Vector3D((float)Math.random()*400.0f,(float)Math.random()*400,(float)Math.random()*0);
      Vector3D vel = new Vector3D(1,0,0);     
      
      particles[i].setPropertiesAndTurnOn(pos,vel);
      ptemp[i]                   = new Particle();
      deriv[i]                   = new Derivative();

    }
    for (int i = 0; i < numOfSprings; i++) {
      s[i]             = new Spring();
      s[i].k           = 0.1;
      s[i].damp        = 0.01;
      s[i].restLength  = 400;
    }
    s[0].from  = 0; 
    s[0].to  = 1;
    s[1].from  = 0; 
    s[1].to  = 2;
    s[2].from  = 0; 
    s[2].to  = 3;
    s[3].from  = 4; 
    s[3].to  = 1;
    s[4].from  = 4; 
    s[4].to  = 2;
    s[5].from  = 4; 
    s[5].to  = 3;
    s[6].from  = 1; 
    s[6].to  = 2;
    s[7].from  = 2; 
    s[7].to  = 3;
    s[8].from  = 3; 
    s[8].to  = 1;
  }








  void idle(float dt) {
    updateParticles(particles, numOfParticles, s, numOfSprings, dt, 1);
    for (int i=0; i<numOfParticles; i++) {
      particles[i].draw();
      
      }
    for (int i = 0; i < numOfSprings; i++) {
      int p1 = s[i].from;
      int p2 = s[i].to;
      stroke(0);
      line(particles[p1].pos.x,particles[p1].pos.y,particles[p1].pos.z, particles[p2].pos.x,particles[p2].pos.y,particles[p2].pos.z);



      /* 
       updateParticles(Particle particles[], int numOfParticles, Spring s[],
       int numberOfSpring, float dt, int method)
       
       */

    }
  }






  void calculateForces (Particle particles[], int numOfParticles, Spring s[], 
  int numOfSprings) {

    int i, p1, p2;
    Vector3D down = new Vector3D(0.0,0.0,-1.0);
    Vector3D zero = new Vector3D(0.0,0.0,0.0);
    float dx,dy,dz, len;
    Vector3D force = new Vector3D(0,0,0);

    for(i=0; i<numOfParticles; i++) {
      particles[i].force = zero;
    //  if(particles[i].fixed) continue;

      /* Gravitation */
      down.mult(gravitational);
      particles[i].force.x += gravitational * particles[i].mass * down.x;
      particles[i].force.y += gravitational * particles[i].mass * down.y;
      particles[i].force.z += gravitational * particles[i].mass * down.z;

      /* Drag */

      particles[i].force.x -= drag * particles[i].vel.x;
      particles[i].force.y -= drag * particles[i].vel.y;
      particles[i].force.z -= drag * particles[i].vel.z;

    }


    /* Handle the spring interaction */
    for (i=0; i<numOfSprings; i++) {
      p1 = s[i].from;
      p2 = s[i].to;
      dx = particles[p1].pos.x - particles[p2].pos.x;
      dy = particles[p1].pos.y - particles[p2].pos.y;
      dz = particles[p1].pos.z - particles[p2].pos.z;
      len = sqrt(dx*dx + dy*dy + dz*dz);

      force.x  =  s[i].k     * (len - s[i].restLength);
      force.x += s[i].damp  * (particles[p1].vel.x - particles[p2].vel.x) * dx / len;
      force.x *= -dx / len;

      force.y  =  s[i].k     * (len - s[i].restLength);
      force.y += s[i].damp  * (particles[p1].vel.y - particles[p2].vel.y) * dy / len;
      force.y *= -dy / len;

      force.z  =  s[i].k     * (len - s[i].restLength);
      force.z += s[i].damp  * (particles[p1].vel.z - particles[p2].vel.z) * dz / len;
      force.z *= -dz / len;

      if (!particles[p1].fixed) {
        particles[p1].force.x += force.x;
        particles[p1].force.y += force.y;
        particles[p1].force.z += force.z;
      }
      if (!particles[p2].fixed) {
        particles[p2].force.x -= force.x;
        particles[p2].force.y -= force.y;
        particles[p2].force.z -= force.z;
      }
    }
  }

  /* Perform one step of the solver */

  void updateParticles(Particle particles[], int numOfParticles, Spring s[],
  int numberOfSpring, float dt, int method) {
    int i;
    

    for (i=0; i<numOfParticles; i++) {
      ptemp[i].force.setXYZ(0,0,0);
    }

    switch (method) {
    case 0:   /* Euler */
      calculateForces(particles,numOfParticles,s,numOfSprings);
      calculateDerivatives(particles,numOfParticles,deriv);
      for (i=0; i<numOfParticles;i++) {
        particles[i].pos.x += deriv[i].dp.x * dt;
        particles[i].pos.y += deriv[i].dp.y * dt;
        particles[i].pos.z += deriv[i].dp.z * dt;
        particles[i].vel.x += deriv[i].dv.x * dt;
        particles[i].vel.y += deriv[i].dv.y * dt;
        particles[i].vel.z += deriv[i].dv.z * dt;
      }
      break;
    case 1:    /* Midpoint */
      calculateForces(particles,numOfParticles,s,numOfSprings);
      calculateDerivatives(particles,numOfParticles,deriv);    
      for(i=0; i< numOfParticles; i++) {
        ptemp[i] = particles[i];
        ptemp[i].pos.x += deriv[i].dp.x * dt / 2;
        ptemp[i].pos.y += deriv[i].dp.y * dt / 2;
        ptemp[i].pos.z += deriv[i].dp.z * dt / 2;
        ptemp[i].vel.x += deriv[i].dv.x * dt / 2;
        ptemp[i].vel.y += deriv[i].dv.y * dt / 2;
        ptemp[i].vel.z += deriv[i].dv.z * dt / 2;
      }
      calculateForces(ptemp,numOfParticles,s,numOfSprings);
      calculateDerivatives(ptemp,numOfParticles,deriv);          
      for (i=0; i < numOfParticles; i++) {
        particles[i].pos.x += deriv[i].dp.x *dt;
        particles[i].pos.y += deriv[i].dp.x *dt;
        particles[i].pos.z += deriv[i].dp.x *dt;
        particles[i].vel.x += deriv[i].dv.x *dt;
        particles[i].vel.x += deriv[i].dv.x *dt;
        particles[i].vel.x += deriv[i].dv.x *dt;
      }
      break;

    }
    for (i=0; i<numOfParticles; i++) {
    //  ptemp[i].force.setXYZ(0,0,0);
    } 
  } 
  
  /*
   Calculate the derivatives
   dp/dt = v
   dv/dt = f / m
*/

  void calculateDerivatives(Particle particles_[], int numOfParticles, Derivative deriv[]) {
    

    for(int i=0; i < numOfParticles; i++) {
      deriv[i].dp.x = particles_[i].vel.x;
      deriv[i].dp.y = particles_[i].vel.y;
      deriv[i].dp.z = particles_[i].vel.z;

      deriv[i].dv.x = particles_[i].force.x / particles_[i].mass;
      deriv[i].dv.y = particles_[i].force.y / particles_[i].mass;
      deriv[i].dv.z = particles_[i].force.z / particles_[i].mass;
    }
  }


}