/************************************************************************
Program name: Nbody
Description: This program does the main Nbody computations,it takes the 
number of particles whose orbits is being calculated, calls the
acceleration subroutine to calculate the subsequent position using the 
leap frog integration scheme.
Output:Returns 0 if successful
Input:none
*************************************************************************/
#include<stdio.h>
#include<file.h>
#include<math.h>
#include<omp.h>
#include<stdlib.h>
/*define the number of particle*/
#define N 2
/*define the mass of the system in grams*/
#define MASS 2.0e+3
/*define the radius of the system in cm*/
#define RADIUS 10
/*define the number of orbits of the system*/
#define ORBIT 1
/*define the time resolution of the problem*/
#define TIME 1e-2
/*Universal gravitational constant*/
#define G 6.67e-8
/*dimensions in the problem*/
#define DIM 2
double fAccel(double*,double,double*,int);
int main()
{
	double dTotTime,dParticleMass,dTimeinc;
	int iTimeSteps;
	dParticleMass=MASS/N;
	double dAccelGravity[N]={0.0};
	FILE * fp1;
	FILE * fp2;
	FILE * accfp;
	/*loop variables*/
	int i,j,k,iCount;
	dTotTime=1/(sqrt(G*MASS/pow(RADIUS,3)));
	
	/*the total time steps is how many time stamps the system will take,
 	  this will be higher if the time resolution of the system 
	 is smaller*/
	
	iTimeSteps=20;
	dTimeinc=dTotTime/(double)iTimeSteps;
	
	
	/*dummy variable used inside the loop to pass values as
	 * a vector to the acceleration subroutine*/
	
	double dPosition[DIM],dOtherPosition[DIM];
	
	/*variables to store the x and y values over varying time
	 * steps of N particles*/
	
	double dPositionkx[iTimeSteps][N];
	double dPositionky[iTimeSteps][N];
	
	/*initial x and y positions of the particles at time -1
	 * for particle 0 and particle 1, they are placed in 
	 * opposite ends of the cooridinate system*/
	
	dPositionkx[0][0]=-10;
	dPositionkx[0][1]=10;
	dPositionky[0][0]=0;
	dPositionky[0][1]=0;
	
	/*x and y positiions at time =0*/

	dPositionkx[1][0]=-9.996;
	dPositionkx[1][1]=9.996;
	dPositionky[1][0]=0.04;
	dPositionky[1][1]=-0.04;
	

	fp1=fopen("particle1.dat","w+");
	fp2=fopen("particle2.dat","w+");
	accfp=fopen("accel.dat","w+");
	
	/*the for loop that does the time integration and calculates
	particle positions for each time step*/

for(k=2;k<iTimeSteps;k++)
{/*for loop for the time step*/
	
   for(i=0;i<N;i++)
   {	
	   /*for loop that cycles through all the particles in order to 
	 * calculate the effect of gravity
	 * The next assignment is the x and y coordinates to 
	 * the dummy variable*/
	

	   dPosition[0]=dPositionkx[k-1][i];
	   dPosition[1]=dPositionky[k-1][i];
	

	for(j=0;j<N;j++)
	{	
		if(i==j)
		{ /*the if construct to check for other particles*/
			continue;
		}
	dOtherPosition[0]=dPositionkx[k-1][j];
	dOtherPosition[1]=dPositionky[k-1][j];
dAccelGravity[i]+=fAccel(dPosition,dParticleMass,dOtherPosition,DIM);
	}
	
	/*Take the acceleration due to gravity and change the position
	 * based on the previous 2 time stamp positions
	 * the time increment is 10^-2 which is the time resolution*/


dPositionkx[k][i]=2*dPositionkx[k-1][i]-dPositionkx[k-2][i]+dAccelGravity[i]*pow(TIME,2);
dPositionky[k][i]=2*dPositionky[k-1][i]-dPositionky[k-2][i]+ dAccelGravity[i]*pow((TIME),2);

printf("the acceleration of particle %d is: %e\n",i,dAccelGravity[i]);

/*this section just writes the position of the particle in
 * different files for different particles*/

if(i==1)
{
fprintf(fp1,"%e",dPositionkx[k][i]);
fprintf(fp1,"\t%e\n",dPositionky[k][i]);
fprintf(accfp,"%e",dAccelGravity[i]);
}

else
{
fprintf(fp2,"%e",dPositionkx[k][i]);
fprintf(fp2,"\t%e\n",dPositionky[k][i]);
fprintf(accfp,"\t%e\n",dAccelGravity[i]);
}
   }
	}
fclose(fp1);
fclose(fp2);
	return 0;
}