// http://www.cs.unc.edu/~andrewz/comp203/hw1/index.html

/*
 gcc gravitytest1.c -lm
 ./a.out > gravitytest1.out
 gnuplot
 gnuplot> set yrange[-1:7];
 gnuplot> set xrange[-1:7];
 gnuplot> plot 'gravitytest1.out' using 1:2, 'gravitytest1.out' using 3:4, \
> 'gravitytest1.out' using 5:6;

*/

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

//#define NUM_PARTICLES  1000
//#define NUM_ITERATIONS 300

#define NUM_PARTICLES  3
#define NUM_ITERATIONS 100

typedef struct ParticleStruct {
  double m;              // particle mass
  double x;              // x-coordinate
  double y;              // y-coordinate
  double vx;             // x-velocity
  double vy;             // y-velocity
  double fx;             // force accumulation on x
  double fy;             // force accumulation on y
} Particle;


/*********************************************************************/
/**   Andrew Zaferakis                                              **/
/**   UNC Chapel Hill                                               **/
/**   COMP 203 Spring 2000                                          **/
/**   Programming Homework 1                                        **/
/**   N-body particle system                                        **/
/**   File:  nbody_particle.cpp                                     **/
/*********************************************************************/


#define PI        3.1415927
#define G         6.673E-11
#define EPSILON   0.0001
#define cube(a)   a*a*a
#define square(a) a*a

void UpdateSystemFast(Particle P[], int numParticles, double dt)
{
  double xdist, ydist, delta_vx, delta_vy, scalar;
  int i, j;

  for (i = 0; i < numParticles; i++)
    P[i].fx = P[i].fy = 0;

  for (i = 0; i < numParticles; i++) {
    for (j = i+1; j < numParticles; j++) {
      xdist    = P[j].x - P[i].x;
      ydist    = P[j].y - P[i].y;
      scalar   = G * P[i].m * P[j].m / (cube(sqrt(square(xdist) + square(ydist))) + EPSILON);
      P[i].fx += scalar * xdist;
      P[i].fy += scalar * ydist;
      P[j].fx -= scalar * xdist;    // must reverse the subtraction for jth
      P[j].fy -= scalar * ydist;
    }
  }

  for (i = 0; i < numParticles; i++) {
    delta_vx = P[i].fx * dt / P[i].m;    // Finalize the ith body, not the jth
    delta_vy = P[i].fy * dt / P[i].m;
    P[i].x  += (P[i].vx + delta_vx / 2) * dt;
    P[i].y  += (P[i].vy + delta_vy / 2) * dt;
    P[i].vx += delta_vx;
    P[i].vy += delta_vy;
  }
}


void UpdateSystemSlow(Particle P[], int numParticles, double dt)
{
  double xdist, ydist, delta_vx, delta_vy, scalar;
  int i, j;

  for (i = 0; i < numParticles; i++) {
    P[i].fx = P[i].fy = 0;
    for (j = 0; j < numParticles; j++) {
      xdist    = P[j].x - P[i].x;
      ydist    = P[j].y - P[i].y;
      scalar   = G * P[i].m * P[j].m / (cube(sqrt(square(xdist) + square(ydist))) + EPSILON);
      P[i].fx += scalar * xdist;
      P[i].fy += scalar * ydist;
    }
  }

  for (i = 0; i < numParticles; i++) {
    delta_vx = P[i].fx * dt / P[i].m;
    delta_vy = P[i].fy * dt / P[i].m;
    P[i].x  += (P[i].vx + delta_vx / 2) * dt;
    P[i].y  += (P[i].vy + delta_vy / 2) * dt;
    P[i].vx += delta_vx;
    P[i].vy += delta_vy;
    printf("\t%f\t%f", P[i].x, P[i].y);
  }
  printf("\n");
}

void InitializeSystem(Particle P[], int numParticles, int config)
{
  int i;
  double tmp;
  if (config == 1)

  for (i = 0; i < numParticles; i++) {
      P[i].m  = 1.0;
      P[i].x  = P[i].y = ((double)i) / ((double)numParticles);
      P[i].vx = P[i].vy = P[i].fx = P[i].fy = 0;
    }
  else if (config == 2)
    for (i = 0; i < numParticles; i++) {
      tmp = ((double)i) / ((double)numParticles);
      P[i].m  = ((double)(numParticles - i))/ ((double)numParticles);
      P[i].x  = 0.5 * (1 + tmp) * cos(2*PI*tmp);
      P[i].y  = 0.5 * (1 + tmp) * sin(2*PI*tmp);
      P[i].vx = P[i].vy = P[i].fx = P[i].fy = 0;
    }
    
   P[0].m = 100000000; 
   P[0].x = 0.0;
   P[0].y = 0.0;
   P[0].vx = 0.0;
   P[0].vy = 0.0;
  
   P[1].m = 100000000;
   P[1].x = 5.0;
   P[1].y = 0.0;
   P[1].vx = 0.0;
   P[1].vy = 0.0;

   P[2].m = 100000000;  
   P[2].x = 2.5;
   P[2].y = 4.33;
   P[2].vx = 0.0;
   P[2].vy = 0.0;
    
}

int CheckMomentum(Particle P[], int numParticles)
{
  double momentum_x = 0;
  double momentum_y = 0;
  int i;

  for (i = 0; i < numParticles; i++) {
    momentum_x += P[i].m * P[i].vx;
    momentum_y += P[i].m * P[i].vy;
  }
  return ((momentum_x <= 0.01) && (momentum_x >= -0.01) && \
	  (momentum_y <= 0.01) && (momentum_y >= -0.01))
    ? 1 : 0;
}

/*********************************************************************/
/**   Main                                                          **/
/*********************************************************************/
int main(int argc, char *argv[])
{
  int i;
  Particle *P = (Particle*)malloc(sizeof(Particle)*NUM_PARTICLES);

  InitializeSystem(P, NUM_PARTICLES, 2);

  // Check the system and report what we find
 
  fprintf(stderr, "Initializing:  Particles = %d\n", NUM_PARTICLES);
  fprintf(stderr, "               N-body using %i bodies and %i iterations of timestep 1\n", NUM_PARTICLES, NUM_ITERATIONS);

  // Begin the timetrials!
  fprintf(stderr, "Starting:      Number of Iterations = %i\n", NUM_ITERATIONS);
  for (i = 0; i < NUM_ITERATIONS; i++)
     UpdateSystemSlow(P, NUM_PARTICLES, 1);
}