/**
 * @file nbody.c
 * @author Sam Moore (20503628) 2012
 * @purpose N-Body simulator - Definition of simulation functions; single threaded version
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include <stdbool.h>
#include <stdlib.h>
#include <unistd.h>

#include "nbody.h"

RUNSTATE runstate = RUN;

/**
 * @function Body_Print
 * @purpose Print the body to a file (or stdout)
 * @param a - Body to print'
 * @param out - FILE* to print to
 */
inline void Body_Print(Body * a, FILE * out)
{
	//fprintf(out,"Body %p M=%f X=%f Y=%f Z=%f Fx=%f Fy=%f Fz=%f Vx=%f Vy=%f Vz=%f\n", 
           //(void*)a, a->mass, a->x[0], a->x[1], a->x[2], a->F[0], a->F[1], a->F[2], a->v[0], a->v[1], a->v[2]);
	fprintf(out, "%f %f %f %f %f %f %f %f %f %f\n", 
		a->mass, a->x[0], a->x[1], a->x[2], a->v[0], a->v[1], a->v[2], a->F[0], a->F[1], a->F[2]);
}

/**
 * @function Body_Force
 * @purpose Calculates the force on a single Body due to all bodies in a System
 * @param a - The Body
 * @param b - The System
 */
inline void Body_Forces(Body * a, System * s) 
{


	for (unsigned i = 0; i < DIMENSIONS; ++i) //Set Force to zero
		a->F[i] = 0;

	for (unsigned index = 0; index < s->N; ++index) // Step through all bodies in the System
	{
		Body * b = s->body+index; //Current body
		if (b == a)
			continue; //Otherwise we would have infinite force
		
		Body_Force(a, b);
	}
}

/**
 * @function Body_Force
 * @purpose Compute force on body a due to body b
 * @param a - Body a
 * @param b - Body b
 */
inline void Body_Force(Body * a, Body * b)
{
	float distance;
	float con;
	float gd;
	//Calculate distance between a and b
	distance = 0.0;
	for (unsigned i = 0; i < DIMENSIONS; ++i)
		distance += square(b->x[i] - a->x[i]);
	distance = sqrt(distance);
	con = G * a->mass * b->mass / square(distance);
	gd = con / distance;
	for (unsigned i = 0; i < DIMENSIONS; ++i) //Add force from b to a's net force
		a->F[i] += gd * (b->x[i] - a->x[i]);
}

/**
 * @function Body_Velocity
 * @purpose Compute velocity of a body
 * @param a - The Body
 */
inline void Body_Velocity(Body * a) 
{
	for (unsigned i = 0; i < DIMENSIONS; ++i)
		a->v[i] += a->F[i] / a->mass * DELTA_T;
}

/**
 * @function Body_Position
 * @purpose Compute position of a body
 * @param a - The Body
 */
inline void Body_Position(Body * a) 
{
	for (unsigned i = 0; i < DIMENSIONS; ++i)
		a->x[i] += a->v[i] * DELTA_T;
}

/**
 * @function System_Compute
 * @purpose Compute forces and then positions for bodies in System
 * NOTE: Only used in the single threaded version of the program
 */
inline void System_Compute(System * s)
{
	System_Forces(s, s);
	System_Positions(s);
}

/**
 * @function System_Forces
 * @purpose Compute forces on each object in System s1 due to all bodies in System s2
 * @param s1, s2 - The two systems
 * NOTE: In single threaded version, called with s1 == s2 == &universe
 *	 In multi threaded version, called with s2 == &universe, but s1 is constructed for each thread
 *		In nested multi-threaded version, s2 is constructed for the nested threads as well
 */
inline void System_Forces(System * s1, System * s2) 
{
	//printf("Compute forces for %p - %d bodies...\n", (void*)s1, s1->N);
	for (unsigned i = 0; i < s1->N; ++i)
	{
		Body_Forces(s1->body+i, s2);
		Body_Velocity(s1->body+i);
	}
	//printf("Compute positions for %p - Done!\n", (void*)s1);
}

/**
 * @function System_Positions
 * @purpose Compute positions of all objects in System
 * @param s - The system
 */
inline void System_Positions(System * s)
{
	//printf("Compute positions for %p - %d bodies...\n", (void*)s, s->N);
	for (unsigned i = 0; i < s->N; ++i)
		Body_Position(s->body+i);
	//printf("Compute positions for %p - Done!\n", (void*)s);
	
}




/**
 * @function System_Init
 * @purpose Initialise a System from an input file
 * @param s - The System
 * @param fileName - The input file
 */
inline void System_Init(System * s, const char *fileName) 
{
	char line[LINE_SIZE];
	char * token;
	FILE * file;

	file = fopen(fileName, "rt");

	if (file == NULL)
	{
		fprintf(stderr, "%s\n", fileName);
		perror("Couldn't open file");
		exit(EXIT_FAILURE);
	}

	s->N = atoi(fgets(line, LINE_SIZE, file));
	s->body = (Body*) calloc((size_t)s->N, sizeof(Body));
	s->steps = 0;

	//printf("----------------------Initial field-------------------------------\n");

	for (unsigned i = 0; i < s->N; ++i)
	{
		if (fgets(line, LINE_SIZE, file) != NULL)
		{
			Body * a = s->body+i;
			token = strtok(line, ",; \t");
			a->mass = atof(token);
			
			for (unsigned j = 0; j < DIMENSIONS; ++j)
			{
				token = strtok(NULL, ",; \t");
				a->x[j] = atof(token);
			}
			for (unsigned j = 0; j < DIMENSIONS; ++j)
			{
				token = strtok(NULL, ",; \t");
				a->v[j] = atof(token);
			}

			// Ignore any additional tokens
			while (token != NULL)
				token = strtok(NULL, ",; \t");
			//Body_Print(a);
		}
	}

	//printf("--------------------------------------------------------------\n");
	
	fclose(file);
}

/**
 * @function Universe_Cleanup
 * @purpose Cleans up the universe by freeing all memory
 *	 Also prints the bodies in the universe to a file specified in "options" if required.
 */
inline void Universe_Cleanup()
{
	//fprintf(stderr, "Called Universe_Cleanup()\n");
	if (options.output != NULL) // An output file was specified...
	{
		FILE * save = NULL;
		if (strcmp("stdout", options.output) == 0)
			save = stdout;
		else if (strcmp("stderr", options.output) == 0)
			save = stderr;
		else
			save = fopen(options.output, "w");
		if (save == NULL)
			perror("Couldn't save universe to file.");
		else
		{
			// Print the output in the same format as the initial field file
			fprintf(save, "%u\n", universe.N);
				
			for (unsigned i = 0; i < universe.N; ++i) // Print all bodies to the file
			{
				Body_Print(universe.body+i, save);			
				
								
			}
			fclose(save);
		}

	}
	free(universe.body); // Cleanup memory used for bodies
	
}

/**
 * @function DisplayStatistics
 * @purpose Display useful information about the computations done,
 *	- Total number of steps computed
 *	- Time in seconds elapsed
 *	- Clock cycles executed 
 *	- Equivelant time for a single thread to execute same number of clock cycles
 */
inline void DisplayStatistics()
{
	clock_t end = clock();
	struct timeval end_time;
	if (gettimeofday(&end_time, NULL) != 0)
	{
		perror("Couldn't get time of day.\n");
		return;
	}
	float runtime = (float)(end_time.tv_sec-options.start_time.tv_sec) + (float)(end_time.tv_usec-options.start_time.tv_usec) / 1.0e6;
	float clock_time = (float)(end - options.start_clock) / (float)(CLOCKS_PER_SEC);
	printf("%u\t%f\t%u\t%f\n", universe.steps, runtime, (unsigned)(end - options.start_clock), clock_time);
}


/**
 * @function ExitCondition
 * @purpose Helper to check whether the program is supposed to exit
 * 	Does not check for user triggered quit
 *	Checks for either a timeout, or computation of the required number of steps
 * 
 * The reason timeouts are integers is because this function is called a lot in lots of threads
 *	So using floats (and working out microseconds every time) is expensive
 *
 */
inline bool ExitCondition(void)
{
	bool result = (runstate != RUN || 
			(options.num_steps >= 0 && universe.steps >= options.num_steps) ||
			(options.timeout >= 0 && ((int)(time(NULL) - options.start_time.tv_sec) >= options.timeout)));
	//fprintf(stderr,"runstate %d\n timeout %d\n steps %d\n", (int)(runstate != RUN), (int)(options.timeout > 0.00 && ((unsigned)(time(NULL) - options.start_time.tv_sec) >= options.timeout)), (int)(options.num_steps > 0 && universe.steps > options.num_steps));
	return result;
}

/**
 * @function Split_System
 * @purpose Helper to divide one system into an array of systems
 *	Each sub system will have N = (s->N / n) bodies in it
 * @param s - The original system (typically &universe)
 * @param n - The number of sub systems in the array
 *
 * WARNING: It is the caller's responsibility to free() the returned array
 */
System * Split_System(System * s, unsigned n)
{
	System * result = (System*)(calloc(n, sizeof(System)));
	if (result == NULL)
	{
		perror("Couldn't create array of sub systems");
		exit(EXIT_FAILURE);
	}

	unsigned n_per_system = (s->N) / n;
	unsigned remainder = (s->N) % n;

	for (unsigned i = 0; i < n; ++i)	
	{
		result[i].N = n_per_system;
		if (i == n-1)
			result[i].N += remainder;
		result[i].body = (s->body) + (n_per_system * i);
		result[i].steps = 0;
	}
	return result;
}

/**
 * @function System_Random
 * @purpose Randomly generate initial body field
 * @param s - The system
 * @param n - Number of bodies
 */
void System_Random(System * s, int n)
{
	srand(time(NULL));
	s->N = (unsigned)n;
	s->body = (Body*) calloc((size_t)s->N, sizeof(Body));
	s->steps = 0;

	s->body[0].mass = 1e11;

	for (unsigned a = 1; a < s->N; ++a)
	{
		s->body[a].mass = 1e7;
		for (unsigned i = 0; i < DIMENSIONS; ++i)
		{
			s->body[a].x[i] = -10000 + rand() % 20000;
			s->body[a].v[i] = -2 + rand() % 4; 
		}
		s->body[a].v[2] = 0;
		s->body[a].x[2] = 0; //set z to zero; force in plane
		
	}
}