/*
 * posfunc.cpp
 * CherenProj
 *
 *  Created by Alexander Sage on 7/25/13.
 *  Supported by Stony Brook University. No rights reserved.
 *
 */

#include <stdio.h>
#include <iostream>
#include <math.h>
#include <time.h>
using namespace std;

#define pie 3.14159

// should return the 'norm' random number between 0 and 1
void mkrnd(double& norm)
{
	double x;
	                  // particles can't travel faster than c ( c=1 )
	x =  rand() ;        
	norm = x / RAND_MAX;
	
	//printf("norm = %f\n", norm);
	
}	


			//	------------------------------------------------------------------------



void posfunc( double (&pos)[4], double (&vel)[4])
{
	//if we want to assume particles come from a small neighborhood around a postion (0,0,0)
	// we will generate a random gaussian distribution around this point
	//the velocity is a randomly generated vector within a cone whose apex is at the origin, 
	// and whose base is a square
	
	// pos is an array that holds the position and time
	//vel is an array that holds velocity and time

	double gdist1, gdist2, norm, vrad, vtheta, vphi, tempvar;
									//gdist1 and gdist2 are placeholder variables to break the formula up into parts
									//norm from mkrnd 
									//vrad, vtheta, vphi are the velocity components of the proton in sphericle coordinates

	pos[3] = 0;        		 // initialize time of array start zero
	vel[3]= 0;
	
	
	for (int i = 0; i < 3; i++) 				//the position of the particles is decided by a randomly
  {																		//generated gaussian distribution
  	
		mkrnd(norm);
		gdist1 =  -0.2*log(norm);
		mkrnd(norm);
		gdist2 =  cos(2*pie* norm);
		pos[i] =  pow(gdist1, 0.5)*gdist2;    	     // Box-Muller method gaussian distribution
		                                    				 //defines elements of array 
	//printf("this is posfunc() pos[ %d ] = %f _  %f _ %f \n", i, pos[i], norm, genhld);
	}
	//printf("pos[3] is %f", pos[3]);
	
	
	mkrnd(norm);
	vrad = norm*0.25+0.75;             //particle will have velocity 0.25 < v < 1
																		 //printf("vrad = %f\n", vrad);
	mkrnd(norm);
	vphi = norm*0.25*pie + 0.375*pie;				//these lines give random location on half the unit sphere
	mkrnd(norm);
	vtheta = norm*0.25*pie + 0.375*pie;	// (0, 0, 1) in cartesian would have vrad =1 and vtheta =0
	vel[0] = vrad*cos(vphi)*sin(vtheta); 		       //vtheta is off z axis and vphi is in xy plane
	vel[1] = vrad*sin(vphi)*sin(vtheta);					 //y is restricted to positive random number
	vel[2] = vrad*cos(vtheta);						// between zero and 1 since the detector is in              
																			// the positive y direction           
																			//x and z restricted to a square whose center is (0, y, 0)
																			//		
																			//printf("vrad = %f , vtheta = %f , vphi = %f , \n 
																			//vel[0] = %f , vel[1] = %f , vel[2] = %f  \n", 
																			//vrad, vtheta, vphi, vel[0], vel[1], vel[2]);
	tempvar = pow(vrad , 2.00);
	if (tempvar >= 1) 
	{
		printf("vrad = %f\n", vrad);
		printf("vrad^2 = %f\n", tempvar);
		printf("\n");
		printf("\n");
		printf("\n");	 printf("Error");
		printf("\n");
		printf("\n");
		printf("\n");
	}
	
	
	
		
}