// http://www.cs.unc.edu/~andrewz/comp203/hw2/index.html
/*
gcc gravity2.c -lm
*/

#include <stdio.h>

//============================================================================
// gravity2.c
//============================================================================

// 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;
   
   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 = 0.0;
   Particles[0].py = 0.0;
   Particles[0].vx = 0.0;
   Particles[0].vx = 0.0;
  
   Particles[1].px = 5.0;
   Particles[1].py = 5.0;
   Particles[1].vx = 0.0;
   Particles[1].vx = 0.0;

}

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=2;

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