/**
 * @file nbody.c
 * @purpose OpenMP version of N-Body simulator, implementation
 * @author Sam Moore (20503628) - 2012
 */

#include "nbody.h"
#include "../single-thread/nbody.c" // Include original code
#include "graphics.h" // For Graphics_Run function only


#ifdef OVER_ENGINEERED
System * sub_system = NULL; 
#endif //OVER_ENGINEERED

/**
 * @function Simulation_Run
 * @purpose Start the simulation
 * @param argc, argv - Passed to Graphics_Run if needed, otherwise unused
 */
void Simulation_Run(int argc, char ** argv)
{	
	
	if (options.num_threads == 0)
		options.num_threads = omp_get_max_threads();

	if (options.draw_graphics)
	{	
		Graphics_Run(argc, argv);
		if (options.pedantic_graphics == false)
		{
			options.num_threads += 1;
			omp_set_nested(1);
			if (omp_get_nested() == 0)
			{
				fprintf(stderr, "Couldn't enable nested parallelism\n.");
				exit(EXIT_FAILURE);
			}	


			#pragma omp parallel sections num_threads(2)
			{
				#pragma omp section
				{
					Compute();
				}
				#pragma omp section
				{
					glutMainLoop();
				}
			}
			return;
		}
	}
	#ifdef OVER_ENGINEERED
	Compute();
	#else
	Compute2();
	#endif //OVER_ENGINEERED



}


void Compute(void) 
{
	omp_set_num_threads(options.num_threads);
	if (omp_get_max_threads() != options.num_threads)
	{
		fprintf(stderr, "Couldn't use %d threads!\n", options.num_threads);
		exit(EXIT_FAILURE);
	}

	#ifdef OVER_ENGINEERED
	sub_system = Split_System(&universe, options.num_threads); // I could use omp for loops instead, but I like this way
	#endif //OVER_ENGINEERED

	#pragma omp parallel
	{
		#ifdef OVER_ENGINEERED
		int index = omp_get_thread_num();
		#endif //OVER_ENGINEERED
		while (!ExitCondition())
		{

			#ifdef OVER_ENGINEERED
				System_Forces(sub_system+index, &universe); 
				#pragma omp barrier // Look! An explicit barrier!
				System_Positions(sub_system+index);
				#pragma omp barrier
			#else
				#pragma omp for
				for (unsigned a = 0; a < universe.N; ++a)
				{
					for (unsigned i = 0; i < DIMENSIONS; ++i)
						(universe.body+a)->F[i] = 0;
					//printf("Thread %d computes force for body %d\n", omp_get_thread_num(), a);
					omp_set_num_threads(options.nested_threads);
					#pragma omp parallel for 
					for (unsigned b = 0; b < universe.N; ++b)
					{
						if (b != a)
							Body_Force(universe.body+a, universe.body+b);
					}
				}
		
				#pragma omp for
				for (unsigned a = 0; a < universe.N; ++a)
				{
					//printf("Thread %d computes position for body %d\n", omp_get_thread_num(), a);
					Body_Velocity(universe.body+a);
					Body_Position(universe.body+a);
				}

			#endif //OVER_ENGINEERED

			#pragma omp master
			{
				// The if statement has a cost, but if can be removed easily in an optimised version.
				// Only the main thread is allowed to mess with glut, because glut is evil (not thread safe).
				// Although this function may be nested, that will only happen if the if would be false.
				// ("if the if would be false"; I am awesome at explaining code)
				if (options.draw_graphics && options.pedantic_graphics)	
				{
					glutMainLoopEvent();
					Graphics_Display();
				}
			}
			#pragma omp single
			{
				// Only one thread should execute this, and we don't care which one.
				// So much easier than the pthread version!
				universe.steps += 1;
				if (options.verbosity != 0 && universe.steps % options.verbosity == 0)
					DisplayStatistics();
			}	
				
		}
	}
}