/* This is a single body test simulation */

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

#define G 6.7*pow(10,-11)
#define m_s 2*pow(10,5)		/* Mass of the sun, in units of 10^25 kg */
#define t_max pow(10,8)		/* stop time */
#define h 1			/* stepsize */


double xdot(double **arrayptr, double delta, short j, unsigned int i);  	/* an easy reference for the first order equations */ 
void RungeKutta4(double**, unsigned int i);		     	/* iteration for forth order (4 step) method */


int main()
{

const unsigned int n = t_max/h;			/* array size */

double** data;				/* this holds the datas */

unsigned int i;						/* loop index */

FILE *fp;

data = malloc(4*sizeof(double*));
for(i=0; i<n; ++i)
	data[i] = malloc( n*sizeof(*data[i]) );

data[0][0] = 0;					/* initial earth x-position */
data[1][0] = 1.5 * pow(10,10);			/* initial earth y-position */
data[2][0] = 3 * pow(10,-6);			/* initial earth x-velocity */
data[3][0] = 0;					/* initial earth y-velocity */


for(i=0; i<n; ++i)
	RungeKutta4(data, i);

fp = fopen("single_body.txt", "w+");
fwrite(data, sizeof(double), 4*n*sizeof(double)/sizeof(double), fp);

return 0;
}


double xdot(double **data0, double delta, short j, unsigned int i)	/* returns a first order ODE from the system of ++j equations */
{
	if(j<2)
		return (data0[j+2][i] + delta);
	else 					/* |r|^3 is calculated twice, but it's a bit messy if you don't */
		return G * m_s * (-(data0[j-2][i] + delta) )/pow(sqrt( pow(-(data0[j][i] + delta),2) + 
		pow(-(data0[j][i] + delta), 2) ), 3);
}

void RungeKutta4(double **datai, unsigned int i)	/* a single 'four'-step iteration  */
{
	short j;
	double s1, s2, s3, s4;
	
	for(j = 0; j<4; ++j)
	{
		s1 = xdot(datai, 0, j, i);
		s2 = xdot(datai, h*s1/2, j, i);
		s3 = xdot(datai, h*s2/2, j, i);
		s4 = xdot(datai, h*s3, j, i);
		
		datai[j][i+1] = datai[j][i] + h*(s1 + 2*s2 + 2*s3 +s4)/6;
	}
}

data = malloc(4*sizeof(double*));
for(i=0; i<n; ++i)
  data[i] = malloc( n*sizeof(*data[i]) );