#include "NBody.hpp"

extern FILE *snapshotsFile, *diagnosticsFile;
extern void* me;						// Pointing to NBody class; 

/////////////////////////////////// NBody::put_snapshot() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
//
// Writes a single snapshot on the output snapshot file.
//
/////////////////////////////////// NBody::put_snapshot() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

void NBody::put_snapshot()
{
	if (snapshotsFile == NULL){
		std::cout << "Error: snapshot file not opened." << std::endl;
		((NBody*)me)->cleanup();
		exit(1);
	}

	// cout.precision(16);                     			// Full double precision;
	fprintf(snapshotsFile, "%d\n", numParticles);		// Total particle number;
	fprintf(snapshotsFile, "%f\n", curr_time_step);		// Current time;

	for (int i = 0; i < numParticles ; i++){
		
		fprintf(snapshotsFile, "%f\t", pos[i*4+3]);		// Mass of particle i;
    	          	
    	for (int k = 0; k < 3; k++)
			fprintf(snapshotsFile, "%f\t", pos[i*4+k]);	// Position of particle i;
       		  	
    	for (int k = 0; k < 3; k++)
      		fprintf(snapshotsFile, "%f\t", vel[i*4+k]);	// Velocity of particle i;
    	
		fprintf(snapshotsFile, "\n");
  	}
	fprintf(snapshotsFile, "\n");
}


/////////////////////////////////// NBody::epot_particle() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
//
// Computes the potential energy of a single particle;
//
/////////////////////////////////// NBody::epot_particle() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

float NBody::epot_particle(int currPart)
{   
   	float ePot = 0.0;
	for (int i=0; i<numParticles; i++){
	   	
		if (currPart != i){

			int idx_curr = currPart*4;	// The current particle's position index in the "pos" vector;
			int idx_other = i*4;		// The other particle's position index in the "pos" vector;

			float r[3], r2= 0.0;
			for (int k = 0; k < 3 ; k++){
				r[k]= pos[idx_other+k] - pos[idx_curr+k];
				r2+= r[k]*r[k];
			}
				
			// To avoid Division by Zero:
			if (r2 == 0) r2= 0.1e-50;

			// ePot = -M1*M2/r;
			ePot += -pos[idx_curr+3]*pos[idx_other+3]/sqrt(r2);
		}
	}

   	return ePot;
}

/////////////////////////////////// NBody::write_diagnostics() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
//
// Writes diagnostics on the diagnostics file:
// current time; number of integration steps so far;
// kinetic, potential, and total energy; absolute and
// relative energy errors since the start of the run.
//
/////////////////////////////////// NBody::write_diagnostics() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

void NBody::write_diagnostics()
{

    float Ekin = 0;       					// The Kinetic energy of the n-body system;
	for (int i = 0; i < numParticles ; i++){
		
		int idx = i*4;

		for (int k = 0; k < 3 ; k++){
			// pos[idx+3] = mass;	
			// Ekin = 1/2*m*(v^2)

            Ekin += 0.5 * pos[idx+3] * vel[idx+k] * vel[idx+k];
		}
	}

    float Epot= 0;							// The Potential energy of the n-body system;
	for (int i=0; i<numParticles; i++)
    	Epot += epot_particle(i);
    
  	Epot *= 0.5;		  					// Against double counting;

    float Etot =  Ekin+Epot;  	 			// The total energy of the actual n-body system;

	if (initFlag){                       	// At first pass, the initial energy is equal to the
		Etot_init = Etot; 					// total energy of the system. We will keep this initial
											// energy for error detections;
	}										

	if (!initFlag){							// At first pass we don't write any diagnostics;
		if (diagnosticsFile == NULL){
			std::cout << "Error: diagnostics file not opened." << std::endl;
			((NBody*)me)->cleanup();
			exit(1);
		}

		// cout.precision(16);                     	// Full double precision;
		fprintf(diagnosticsFile, "At time t = %f, after %ld steps:\n", curr_time_step, taken_steps);
		fprintf(diagnosticsFile, " E_init = %f,\n E_tot = %f,\n E_kin = %f,\n E_pot = %f\n", Etot_init, Etot, Ekin, Epot);
		fprintf(diagnosticsFile, "                Absolute energy error: E_tot - E_init = %f\n", Etot - Etot_init);
		fprintf(diagnosticsFile, "                Relative energy error: (E_tot - E_init) / E_init = %f\n", (Etot - Etot_init) / Etot_init);
	}else{
		initFlag = false;
	}
}