/**
 * @file nbody.c
 * @purpose Barnes Hut 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

#include "tree.h"



omp_lock_t runstate_lock;

#ifdef USE_THREADS
omp_lock_t nested_lock;
#endif //USE_THREADS

Node * parent = NULL;

/**
 * @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)
{	
	
	omp_init_lock(&runstate_lock);
	#ifdef USE_THREADS
	omp_init_lock(&nested_lock);
	#endif //USE_THREADS

	#ifdef DUAL_UNIVERSE
	if (!options.pedantic_graphics)
	{
		fprintf(stderr, "Warning: Using pedantic graphics to ensure alternate reality matches our own (as closely as possible).\n");
		options.pedantic_graphics = true;
	}
	#endif //DUAL_UNIVERSE

	omp_set_nested(1);
	if (omp_get_nested() == 0)
	{
		fprintf(stderr, "Couldn't enable nested parallelism\n.");
		exit(EXIT_FAILURE);
	}	

	if (options.num_threads == 0)
	{
		options.num_threads = omp_get_max_threads();
		if (options.draw_graphics)
			options.num_threads -= 1;
	}
	if (options.nested_threads == 0)
	{
		options.nested_threads = omp_get_max_threads() - options.num_threads;
		if (options.nested_threads == 0)
			options.nested_threads = 1;
	}

	#ifdef DUAL_UNIVERSE
		System_Init(&alternate_universe, options.input);
	#endif //DUAL_UNIVERSE

	if (options.draw_graphics)
	{	
		Graphics_Run(argc, argv);
		if (options.pedantic_graphics == false)
		{
			


			#pragma omp parallel num_threads(2)
			{
				// This can't be done with sections!
				// Because glut is useless, and can only be dealt with in the main thread
				// I just hope that thread "0" is always the main thread!
				if (omp_get_thread_num() == 0) // #pragma omp section
				{
					//printf("Thread %d gets graphics\n", omp_get_thread_num());
					glutMainLoop();
					QuitProgram(false); // Program gets to here if the window is quit
				} 
				else // #pragma omp section
				{
					//printf("Thread %d gets computations\n", omp_get_thread_num());
					Compute();
				}

			}
			return;
		}
	}
	Compute();

	omp_destroy_lock(&runstate_lock);
	#ifdef USE_THREADS
	omp_destroy_lock(&nested_lock);
	#endif //USE_THREADS
	


}

/**
 * @function Compute
 * @purpose The main computation loop
 * 	Uses threading, but only on the first level
 *	I tried to introduce threading to some of the recursive functions in the tree...
 *	But failed
 */
void Compute(void) 
{
	omp_set_nested(1);
	if (omp_get_nested() == 0)
	{
		fprintf(stderr, "Couldn't enable nested parallelism\n.");
		exit(EXIT_FAILURE);
	}	

	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);
	}

	parent = Node_Create(NULL); 
	
	Generate_Tree(&universe, parent);




#pragma omp parallel
{
	while (!ExitCondition())
	{
		#pragma omp single
		{
			Node_CalculateMass(parent);
		}

		#pragma omp for // The easiest way to parallelise it
		for (unsigned a = 0; a < universe.N; ++a)
		{
			for (unsigned i = 0; i < DIMENSIONS; ++i)
				universe.body[a].F[i] = 0;

			Node_Forces(universe.body+a, parent, options.theta);
			//printf("Forces on %p done\n", (void*)(universe.body+a));
			Body_Velocity(universe.body+a);	
		}

		// put a nowait here so that the first thread out can start updating the tree
		#pragma omp for nowait 
		for (unsigned a = 0; a < universe.N; ++a)
		{

			Body_Position(universe.body+a);
		}

		#pragma omp single // Ideally the first thread out should do this
		{
			Update_Tree(&universe, parent);
		}

		// I haven't noticed any difference in performance.
		// I'm just showing off.
		#pragma omp barrier

		// Only the master may do stuff with graphics
		#pragma omp master
		{
			if (options.draw_graphics && options.pedantic_graphics)
			{
				glutMainLoopEvent();
				Graphics_Display();
			}
		}

		// Someone should do this stuff too
		#pragma omp single
		{
			universe.steps += 1;
			if (options.verbosity != 0 && universe.steps % options.verbosity == 0)
				DisplayStatistics();
		}
		
	}
}
	Node_Destroy(parent);
	free(parent);

	QuitProgram(false);
}


void QuitProgram(bool error)
{
	if (runstate == QUIT || runstate == QUIT_ERROR)
		return;
	omp_set_lock(&runstate_lock);
	runstate = (error) ? QUIT_ERROR : QUIT;
	omp_unset_lock(&runstate_lock);
}

#ifdef DUAL_UNIVERSE
void BeforeDraw()
{
	// Before drawing the main universe... compute and draw the alternate universe
	// Not threaded, because I am too lazy
	System_Compute(&alternate_universe);
	
	for (int i = 0; i < alternate_universe.N; ++i) 
	{
	
		Body * b = alternate_universe.body+i;
		Body * original = universe.body+i;
		//glColor3f(0.0f, b->mass/1e11*100, 0.0f);
		glColor3f(1.0f, 0.0f, 1.0f);
		glPushMatrix(); // to save the current matrix
		glTranslated(scale*b->x[0], scale*b->x[1], scale*b->x[2]);
		glutSolidSphere (BALL_SIZE, 10, 10);
		glPopMatrix(); // restore the previous matrix

		// Draw a line from the "alternate" body to the real one
		glBegin(GL_LINES);
		glVertex3f(scale*b->x[0], scale*b->x[1], scale*b->x[2]);
		glVertex3f(scale*original->x[0], scale*original->x[1], scale*original->x[2]);
		glEnd();
	}
}
#endif //DUAL_UNIVERSE

#ifndef AfterDraw
void AfterDraw()
{
	/* // Uncomment to draw cubes
	glPolygonMode( GL_FRONT_AND_BACK, GL_LINE );
	Draw_Node(parent);
	glFlush();
	glPolygonMode( GL_FRONT_AND_BACK, GL_FILL);
	*/
}
#endif //AfterDraw