// http://www.cs.unc.edu/~andrewz/comp203/hw2/index.html
/*
gcc mygravity7.c -lm
./a.out > gravity7.out

This version: 12 particles, (0,0) mass 100000000
set xrange [-1:10]
set yrange [-1:10]

40 time steps. Plot with Gnuplot
plot 'gravity7.out' using 1:2, 'gravity7.out' using 3:4, \
'gravity7.out' using 5:6, 'gravity7.out' using 7:8, \ 
'gravity7.out' using 9:10, 'gravity7.out' using 11:12;

*/

#include <stdio.h>

//============================================================================
// mygravity5.cpp
//============================================================================

// http://www.cs.unc.edu/~andrewz/comp203/hw2/gravity/gravity.cpp

#include <stdlib.h>
#include <math.h>

#define G 6.673e-11          // GRAVITATIONAL CONSTANT
#define PI 3.1415926535898   // CONSTANT PI

struct Particle
{
  double px, py;
  double vx, vy;
  double Mass;
  double ForceAccumX, ForceAccumY;
};

void ComputeGravForces_Fast(struct Particle Particles[], int NumParticles);
void UpdateParticles_MidPoint(struct Particle Particles[], int NumParticles, 
		float DeltaT);

void CheckMomentumConservation(struct Particle Particles[], int NumParticles);

void InitConfig1(struct Particle Particles[], int NumParticles);

void ComputeGravForce(struct Particle A, struct Particle B, double *ForceX, double *ForceY)
{
  double ABx = B.px - A.px;
  double ABy = B.py - A.py;
  double DistAB2 = ABx*ABx+ABy*ABy;
  double DistAB = sqrt(DistAB2);
  double DistAB3 = DistAB2 * DistAB;
  double Scale = (G * A.Mass * B.Mass) / DistAB3;
  *ForceX = Scale * ABx;
  *ForceY = Scale * ABy;
}

void ComputeGravForces_Fast(struct Particle *Particles, int NumParticles)
{
  int i,j;
  double ForceX, ForceY;


  for (i=0; i<NumParticles; i++) {
    Particles[i].ForceAccumX=0;
    Particles[i].ForceAccumY=0;
  }

  for (i=0; i<(NumParticles-1); i++)
    for (j=i+1; j<NumParticles; j++) {
      ComputeGravForce(Particles[i],Particles[j],&ForceX,&ForceY);
      Particles[i].ForceAccumX+=ForceX;
      Particles[i].ForceAccumY+=ForceY;
      Particles[j].ForceAccumX-=ForceX;
      Particles[j].ForceAccumY-=ForceY;
    }
}

void UpdateParticles_MidPoint(struct Particle *Particles, int NumParticles, 
		 float DeltaT)
{
  double Scale;
  double dVx;
  double dVy;
  int i;
  	 
  for (i=0; i<NumParticles; i++)  {
    Scale = DeltaT / Particles[i].Mass;
    dVx = Particles[i].ForceAccumX * Scale;
    dVy = Particles[i].ForceAccumY * Scale;
    Particles[i].px += (Particles[i].vx + dVx*0.5) * DeltaT;
    Particles[i].py += (Particles[i].vy + dVy*0.5) * DeltaT;
    Particles[i].vx += dVx;
    Particles[i].vy += dVy;
//    printf("Part. %d Mass: %.1f  pos: %f, %f and vel: %f,%f\n", 
//		    i, Particles[i].Mass, Particles[i].px,
//		Particles[i].py, Particles[i].vx, Particles[i].vy);
    printf("\t%f\t%f", Particles[i].px, Particles[i].py);

  }
  printf("\n");
}

void CheckMomentumConservation(struct Particle *Particles, int NumParticles)
{
  double TotalMomentumX=0, TotalMomentumY=0;
  int i;
  
  for (i=0; i<NumParticles; i++)
  {
    TotalMomentumX += (Particles[i].Mass * Particles[i].vx);
    TotalMomentumY += (Particles[i].Mass * Particles[i].vy);
  }  
  printf("%g %g\n",TotalMomentumX,TotalMomentumY);
}

void InitConfig1(struct Particle *Particles, int NumParticles)
{
   double t; 
   int i;
   double centerx, centery;
   
   for (i=0; i<NumParticles; i++)  {
//    double t = ((double)i) / ((double)NumParticles);
    t = (double) rand()/RAND_MAX;
    Particles[i].Mass = 100000000.0;
//    Particles[i].px = t;
//    Particles[i].py = t;
//    Particles[i].vx = 0;
//    Particles[i].vy = 0;
  }

   Particles[0].px = 1.0;
   Particles[0].py = 3.0;
   Particles[0].vx = 0.0;
   Particles[0].vy = 0.0;
  
//   Particles[1].Mass = 100;

   Particles[1].px = 4.0;
   Particles[1].py = 3.0;
   Particles[1].vx = 0.0;
   Particles[1].vy = .0;

   Particles[2].px = 2.5;
   Particles[2].py = 5.6;
   Particles[2].vx = 0.0;
   Particles[2].vy = 0.0;

   Particles[3].px = 6.0;
   Particles[3].py = 3.0;
   Particles[3].vx = 0.0;
   Particles[3].vy = 0.0;
   
   Particles[4].px = 9.0;
   Particles[4].py = 3.0;
   Particles[4].vx = 0.0;
   Particles[4].vy = 0.0;
   
   Particles[5].px = 7.5;
   Particles[5].py = 5.6;
   Particles[5].vx = 0.0;
   Particles[5].vy = 0.0;

   
/*   
   Particles[0].px = 1.0;
   Particles[0].py = 7.0;
   Particles[0].vx = 0.0;
   Particles[0].vy = 0.0;
  
//   Particles[1].Mass = 100;

   Particles[1].px = 4.0;
   Particles[1].py = 7.0;
   Particles[1].vx = 0.0;
   Particles[1].vy = .0;

   Particles[2].px = 2.5;
   Particles[2].py = 9.6;
   Particles[2].vx = 0.0;
   Particles[2].vy = 0.0;

   Particles[3].px = 6.0;
   Particles[3].py = 7.0;
   Particles[3].vx = 0.0;
   Particles[3].vy = 0.0;
   
   Particles[4].px = 9.0;
   Particles[4].py = 7.0;
   Particles[4].vx = 0.0;
   Particles[4].vy = 0.0;
   
   Particles[5].px = 7.5;
   Particles[5].py = 9.6;
   Particles[5].vx = 0.0;
   Particles[5].vy = 0.0;
   
   Particles[6].px = 6.0;
   Particles[6].py = 2.0;
   Particles[6].vx = 0.0;
   Particles[6].vy = 0.0;
   
   Particles[7].px = 9.0;
   Particles[7].py = 2.0;
   Particles[7].vx = 0.0;
   Particles[7].vy = 0.0;
  
   Particles[8].px = 7.5;
   Particles[8].py = 4.6;
   Particles[8].vx = 0.0;
   Particles[8].vy = 0.0;

   
   Particles[9].px = 1.0;
   Particles[9].py = 2.0;
   Particles[9].vx = 0.0;
   Particles[9].vy = 0.0;
   
   Particles[10].px = 4.0;
   Particles[10].py = 2.0;
   Particles[10].vx = 0.0;
   Particles[10].vy = 0.0;
   
   Particles[11].px = 3.5;
   Particles[11].py = 4.6;
   Particles[11].vx = 0.0;
   Particles[11].vy = 0.0;
*/
/*   
  Particles[NumParticles-1].Mass = 100000000.0;
   Particles[NumParticles-1].px = 7.5;
   Particles[NumParticles-1].py = 9.6;
   Particles[NumParticles-1].vx = 0.0;
   Particles[NumParticles-1].vx = 0.0;
*/
  centerx = (100000000*0.0 + 100000000*2.5 + 100000000*5.0)/300000000;
  printf("Center of mass x: %f\n", centerx);
  centery = (100000000*0.0 + 100000000*0.0 + 100000000*4.33)/300000000;
  printf("Center of mass y: %f\n", centery);
}

void RunSimulation(int NumParticles, int NumIterations, double TimeStep)
{
  int k;
  
  struct Particle *Particles = (struct Particle*)malloc(sizeof(struct Particle)*NumParticles);

  InitConfig1(Particles,NumParticles);

  for (k=0; k<NumIterations; k++) {
      ComputeGravForces_Fast(Particles,NumParticles);
      UpdateParticles_MidPoint(Particles,NumParticles,TimeStep);
  }

  free(Particles);
}

int main()
{
  int NumParticles=6;

  RunSimulation(NumParticles,40, 1);
  return 0;
}