/*
 * geom.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>
using namespace std;

#define pie 3.14159

//geometrical shape of LQBar
void geom(double (&pos)[4], bool& htbl, bool& leftvac)
{
	double dtByf, dtByb, dtBxl, dtBxr, dtBzb, dtBzt;
	double dtTyf, dtTyb, dtTxl, dtTxr, dtTzb, dtTzt;
	//the LQbar will have a small square region 'facing' the origin
	// a rectanglular region will go out/back away from the square (along y-axis)
	// then a 45 degree angle will be defined, after which some particles will be reflected
	// up a rectangular region that is perpendicular to the first (up z-axis), 
	//possibly striking the square detector
	
	// the bottom rectangular region
	dtByf = 2;
	dtByb = 3.25;
	dtBxl = -0.25;
	dtBxr = 0.25;
	dtBzb = -0.25;
	dtBzt = 0.25;
	
	//the top rectangular region
	dtTyf = 2.75;
	dtTyb = dtByb;
	dtTxl = dtBxl;
	dtTxr = dtBxr;
	dtTzb = dtBzt;
	dtTzt = dtBzt + 1;      //1.25
	
	////define [0]=x [1]=y [2]=z [3]=t
	
	//checks enter detector
	//requirements for enter detector for top and bottom
	if ((pos[1] >= dtByf) && (pos[0] >= dtBxl) && (pos[0] <= dtBxr) && (pos[1] <= dtByb)  
			&& (pos[2] >= dtBzb) && (pos[2] <= dtTzt))     //requirements for enter detector
	{
		if ( ((pos[1] <= dtTyf) && (pos[2] <= dtBzt)) || ( (pos[1] > dtTyf) ))   //bottom || top
		{htbl = true;}          //in detector
		else 
		{htbl = false;}        //outside detector
		
	}
	else 
	{htbl = false;}        //outside detector
	
	if ( (pos[0] > dtBxr) || (pos[0] < dtBxl) || (pos[2] < dtBzb) || 
			((pos[2] > dtBzt) && (pos[1] > dtTyf)) || (pos[1] > dtTyb) || (pos[2] > dtTzt) ) 
	                                           	//old if statement ((pos[0])*(pos[0]) + (pos[2])*(pos[2]) ) >= 1
	{leftvac = true;}                           // i don't know why that old condition was there
		
			//printf("x is %f \n", pos[0]);
			//printf("x^2  is %f \n", (pos[0])*(pos[0]));
			//printf("z is %f \n", pos[2]);
			//printf("z^2  is %f \n", (pos[2])*(pos[2]));
			//printf("particle absorbed into wall\n");
	
}


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


/**
double diagonal(double y)
{
	double yn = y-2.75;
	double z = y;
	double zn = z + 0.25;
	return zn;
}*/


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


//z - up; y is forward; x is side to side;
//reflect photon off walls
void rfctp (double (tmppos)[4], double (photpos)[4], double (&photvel)[4])
{
	
	double dtByf, dtByb, dtBxl, dtBxr, dtBzb, dtBzt;
	double dtTyf, dtTyb, dtTxl, dtTxr, dtTzb, dtTzt;
	
	// the bottom rectangular region
	dtByf = 2;
	dtByb = 3.25;
	dtBxl = -0.25;
	dtBxr = 0.25;
	dtBzb = -0.25;
	dtBzt = 0.25;
	
	//the top rectangular region
	dtTyf = 2.75;
	dtTyb = dtByb;
	dtTxl = dtBxl;
	dtTxr = dtBxr;
	dtTzb = dtBzt;
	dtTzt = dtBzt + 1;   //1.25
	
	
	//implement angled slide
	//have set positions
	//box to put the angle in is total; .5 along y(start 2.75); 0.5 along z(middle 0);
	
	//double z = diagonal(tmppos[1]);
	/**if (tmppos[2] <= z) 
	{
		photvel[2] = -1* photvel[2];
		photvel[1] = -1* photvel[1];
	}*/
	//this needs some work. Its close but the code doesn't stop running with it
	if ((tmppos[0] < dtBxl) || (tmppos[0] > dtBxr))     //reflects from walls in y-z plane
	{
		photvel[0] = -1* photvel[0];
	}
	//top detector
	else if((tmppos[2] > dtBzt) && (tmppos[1] > dtTyb))
	{
			photvel[1] = -1* photvel[1];    //reflects from walls in x-z plane
	}
	else if(tmppos[2] > dtTzt)
	{
		photvel[2] = -1* photvel[2];
	}
	//bottom detector
	else if((tmppos[2] < dtBzt))
	{
		if ((tmppos[1] < dtByf) || (tmppos[1] > dtByb))      //reflects from walls in x-z plane
		{
			photvel[1] = -1* photvel[1];
		}
		else if ((tmppos[2] < dtBzb) || ((tmppos[1] < dtTyf) && (tmppos[2] > dtBzt)))      //reflects from walls in x-y plane
		{
			photvel[2] = -1* photvel[2];
		}
	}
	
	
	//whole if statement is about that stupid corner
	else if ((tmppos[1] < dtTyf) && (tmppos[2] > dtBzt))  
	{
		if ((photpos[2] > dtBzt) && (photpos[1] > dtTyf)) 
		{
			photvel[1] = -1* photvel[1];				//using photon's last position inside the detector to tell
		}																			//where the photon is coming from and where it should be reflected
		else if ( (photpos[2] < dtBzt) && (photpos[1] < dtTyf) )
		{
			photvel[2] = -1* photvel[2];
		}
		else if ((photpos[2] < dtBzt) && (photpos[1] > dtTyf) )
		{
			
			if ( (dtTyf - tmppos[1]) < (tmppos[2] - dtBzt) ) 
			{
				photvel[1] = -1* photvel[1];
				printf("reflect y\n");
				photvel[1] -= 0.01;
			}
			else if ((dtTyf - tmppos[1]) > (tmppos[2] - dtBzt)) 
			{
				photvel[2] = -1* photvel[2];
				printf("reflect z\n");
				photvel[2] -= 0.01;
			}
			if ( (tmppos[1] - dtTyf) > (dtBzt - tmppos[2]) ) 
			{
				photvel[1] = -1* photvel[1];
				//printf("reflect y\n");
				photvel[1] -= 0.01;
			}
			else if ( (tmppos[1] - dtTyf) < (dtBzt - tmppos[2]) ) 
			{
				photvel[2] = -1* photvel[2];
				//printf("reflect z\n");
				photvel[2] -= 0.01;
			}
			
			else 
			{
				printf("unresolved ERROR (rfctp)\n");
			}

			/**printf("Warning: poor reflection (rfctp) at corner\n");
			for (int j = 0; j < 4; j++) 
			 {
			 printf("photpos[ %d ] = %f  |   tmppos[ %d ] = %f\n", j, photpos[j], j, tmppos[j]);
			 }
			printf("\n");*/
		}
		else 
		{
			printf("Huge Error (rfctp)\n \n");
		}


	}
	
	else if ((tmppos[2] == dtBzt) && (tmppos[1] == dtTyf))   //if it hits the corner of the two walls exactly
	{
		photvel[2] = -1* photvel[2];
		photvel[1] = -1* photvel[1];
	}
	
}

//what I probably should have done for reflection is to draw a line from tmppos to photpos 
//find the wall it intersects, take the point on the wall, set that point as the particle's
//new position, then invert the velocity



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



//this function finds the wavelength dependant index of refraction
//based on the Sellmeier equation. taken from refractiveindex.info

void specrel(double& index)
{	
	double lambda, tmpvar;
	lambda = 3.0 * pow(10., -1.);      
	index = sqrt( 1 + 0.6961663*pow(lambda,2.)/(pow(lambda,2.)-pow(0.0684043,2.))
							 + 0.4079426*pow(lambda,2.)/(pow(lambda,2.)-pow(0.1162414,2.)) 
							 + 0.8974794*pow(lambda,2.)/(pow(lambda,2.)-pow(9.896161,2.)) );
	//printf("index = %f \n", index);
}


//this function gets Beta based on the proton's velocity
void specrel(double (vel)[4], double& Beta)
{
	double tmpvar;
	tmpvar = pow(vel[0] , 2.00) + pow(vel[1] , 2.00) + pow(vel[2] , 2.00);
	Beta = 1/(pow( (1-tmpvar), 0.50));
	
}




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




//emitting cherenkov light
void projrad(double (vel)[4], double (&photvel)[4], double alpha, double zeta, double index, bool& tosmll)
{
	
	double  tmpvar, 
        	prad, ptheta, pphi, photrad, phottheta, photphi;

	
	tmpvar = pow(vel[0] , 2.00) + pow(vel[1] , 2.00) + pow(vel[2] , 2.00);
	prad = pow(tmpvar , 0.50);
	pphi = acos((vel[0])/prad);						//now we have velocity of particle vel[ ] , converted into spherical
	ptheta = acos((vel[2])/prad) ;        // coordinates
	
	if (zeta <= 0.5*pie) 
	{
		phottheta = ptheta - alpha * pow(cos(zeta), 2.0);    //phi is in the x-y plane, theta is off z axis
		photphi = pphi + alpha *pow(sin(zeta), 2.0);
	}
	else if ((zeta <= pie) & (zeta > 0.5*pie)) 
	{
		phottheta = ptheta + alpha * pow(cos(zeta), 2.0);
		photphi = pphi + alpha *pow(sin(zeta), 2.0);
	}
	else if ((zeta > pie) & (zeta <= 1.5*pie))    				//the velocity of the photon emitted from cherenkov
	{																										//radiation is found using the vector for velocity 
		phottheta = ptheta + alpha * pow(cos(zeta), 2.0);  //of particle emitting it, then taking every
		photphi = pphi - alpha *pow(sin(zeta), 2.0);			// vector at an angle, alpha, off of the 
	}																										//radial part of that vector in increments denoted by i
	else if (zeta > 1.5*pie) 														// these if statements change the phi and theta
	{																										//components of the photon velocity vector, 
		phottheta = ptheta - alpha * pow(cos(zeta), 2.0); // but the radial component of the photon is based  
		photphi = pphi - alpha *pow(sin(zeta), 2.0);			// on the phase velocity of light in that medium
	}
	else 
	{
		printf("Error in projection of alpha\n \n \n");
		tosmll = 1;
	}
	
	photrad = 1/index;   
	photvel[0] = photrad* cos(photphi);
	photvel[1] = photrad* sin(photphi);
	photvel[2] = photrad* cos(phottheta);
	photvel[3] = vel[3];
	
	if (prad < photrad) 
	{
		tosmll = 1;
		//printf("particle moving too slow; prad = %f\n", prad);
		//printf("phase velocity is = %f \n", photrad);
	}
	
}