#include <iostream>
#include <cstdlib>
#include <random>
#include <math.h>
#include <TGraph.h>
#include <TCanvas.h>


double Prob_IP() {
	double b = ((double)rand()/(double)RAND_MAX) * 20.0;
	return b;
}

int Probability(double R) {
	double e = 2.71828182846;
	if ((double)rand()/(double)RAND_MAX <= (0.1693)/(pow(e, (R-6.38)/0.535) + 1)) {
		return 1;
	}
	else {
	return 0;
	}
}

double Random() {
	double a = (((double)rand()/(double)RAND_MAX)*14) - 7;
	return a;
}

int main() {
	// "At separation <B>, # of Participants, # of Collisions"
	int Num_P = 0;
	int Num_C = 0;
	int B = 6;
	int Num_N1C = 0;
	int Num_N2C = 0;
	std::cout << "Enter Separation of Nuclei (in fermi): ";
	std::cin >> B;

	
	//Constants
	double pi = 3.14159265358979;

	// Arrays for Nucleus 1
	int N1 [197] = {0};
	double X1 [197] = {0};
	double Y1 [197] = {0};
	double Z1 [197] = {0};
	
	// Arrays for Nucleus 2
	int N2 [197] = {0};
	double X2 [197] = {0};
	double Y2 [197] = {0};
	double Z2 [197] = {0};
	
	// Temp Variables for calculation
	double x = 0.0;
	double y = 0.0;
	double z = 0.0;
	double Rmag = 0.0;
	int counter = 0;
	int counter2 = 0;
	double d2 = 0.0;

	srand( time( NULL ));


	//Initializes Nucleus #1 and stores in X1 Y1 and Z1 Arrays
	for ( ; counter < 197; ) {
		x = Random();
		y = Random();
		z = Random();
		Rmag = sqrt(x*x + y*y + z*z);
		if (Rmag <= 7.0 && Probability(Rmag) == 1) {
			X1[counter] = x;
			Y1[counter] = y;
			Z1[counter] = z;
			counter += 1;
		}


	}

	//Initializes Nucleus #2 and stores in X2 Y2 and Z2 Arrays
	for ( ; counter2 < 197; ) {
		x = Random();
		y = Random();
		z = Random() - B;
		Rmag = sqrt(x*x + y*y + (z + B)*(z + B));
		if (Rmag <= 7.0 && Probability(Rmag) == 1) {
			X2[counter2] = x;
			Y2[counter2] = y;
			Z2[counter2] = z;
			counter2 += 1;
		}
	}
	
	//Compares each pair of Nucleons to see if a collision occurred, and to count the participants
	for (int i = 0; i < 197; i++) {
		for (int j = 0; j < 197; j++) {
			d2 = (pow(Y1[i] - Y2[j], 2.0) + pow(Z1[i] - Z2[j], 2.0));
			if ( d2 < 4.0/pi ) {
				N1[i] += 1;
				N2[j] += 1;
			}
		}
	}

	//For loops to count number of participants and number of collisions
	for (int k = 0; k < 197; k++) {
		Num_N1C += N1[k];
		Num_N2C += N2[k];
		if (N1[k] > 0) {
			Num_P += 1;
		}
		if (N2[k] > 0) {
			Num_P += 1;
		}
	}
	if (Num_N1C != Num_N2C) {
		std::cout << "Num_N1C = " << Num_N1C <<" and Num_N2C = " << Num_N2C << std::endl;
break;
	}
	else {
		Num_C = Num_N1C;
	}
	//std::cout << "Total Participants: " << Num_P << std::endl << "Total Collisions: " << Num_C << std::endl;


	/*TGraph *gr1 = new TGraph(197, X1, Y1);
	TCanvas *c1 = new TCanvas("c1","Graph Draw Options", 600, 800);
	gr1->Draw("p");
	c1->SaveAs("/home/jmerges/public_html/tests/c1.png");
	*/
return 0;
}