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

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

#define pie 3.14159
#define c 2.99792458*pow(10.00, 8.0)
#define charge 1.602176565*pow(10.00, -19.0)
#define coulombConst  8.987551787*pow(10.00, 9.0)
#define massOfProton  1.67262177*pow(10.00, -27.00)
#define epsilon  3.8
#define lattice 4.9133*pow(10.0, -10.0)

                          //function to slow the particle
void slwp(double (&vel)[4], int& mvg)
{
	//printf("slwp called \n");
	double hyps,shyps, theta, phi, radiusxy, r ;
																	//lattice is the lattice spacing of quartz
																  //epsilon is the electric permittivity of fuzed quarts
	/*
	printf("before adjustment\n");
	for (int i = 0 ; i<3; i++) 
	{
		printf("vel[ %d ] = %f \n", i, vel[i]);
	}
	 */
	hyps = 0; theta = 0; phi = 0; radiusxy = 0; shyps =0;
	r = lattice*0.5*0.32835;           //0.32835 is the average value of x^2 from 0 to 1
	
	for (int i = 0; i < 3; i++) 				//it takes the velocity, gets the radial component whose max is 1
	{																		//coulomb repulsion calculated to subtract from the velocity
		hyps += pow(vel[i], 2.00);				//
	}																		// 
	radiusxy = sqrt(pow(vel[0], 2.00) + pow(vel[1], 2.00));  
	shyps = sqrt(hyps);
	phi = acos(vel[2]/shyps);          //spherica coordinates: phi is off positive z axis.
	theta = acos(vel[0]/radiusxy);     //theta is around the x-y plane
	//vf^2 = vi^2 - (2/m)*(kq^2/r^2epsilon)*d 
	hyps = hyps - (2.0*lattice*(pow(charge, 2.00))*coulombConst)/(5000.0*c*massOfProton*epsilon*pow(r, 2.00));          //Coulomb Repulsion of two unitary charges
	shyps = sqrt(hyps);
	vel[2] = cos(phi)*shyps;
	vel[1] = sin(theta)*sin(phi)*shyps;
	vel[0] = cos(theta)*sin(phi)*shyps;
  /* 
	for (int i = 0; i<4; i++) 
	{
		printf("vel [ %d ] = %f\n", i, vel[i]);
	} 
	 */
	
	mvg = 0;														//determines if the particle is moving fast enough to be considered
	
	for (int i = 0 ; i<3; i++) 
	{
		if (vel[i] < 0.01)           //if mvg = three, particle has ~ stopped. 
		{
			mvg ++;
		}
		//printf("vel[ %d ] = %f \n", i, vel[i]);
	}
	if(hyps < 0)
	{
			mvg = 3;
	}
	
	
	//Frank Tamm formula
	//dE/(dwdx) = (e^2)*w/(c^2)*[1-1/((Beta^2)*(n(w))^2]
}
/*
 The magnitude is reduced by the coulomb force (0.003657065764756304). 
 an average of descrete points is taken on a graph of r^2 from 0 to 1. the zero to 1 represents the
 middle between two atoms in the lattice (zero) to the next closest atom (1). which is ~0.3
 E = 0.5mv^2 was used along with kqq/r^2 to find the drop in energy from the particle as it traveled.
 coulomb force was multiplied by the distance from one atom to another. for energy lost from work.
 the epsilon in k was 3.8 since it passes through matter.
 The force of the coulomb repulsion was incrediably large and I didn't know what to do, so I divided it by 
 the speed of light and then again by 500
 */



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




                             //track the photons emitted
int tphoton(double (vel)[4], double (pos)[4])
{
	double photpos[4], tmppos[4], photvel[4], alpha, zeta = 0, index, time, Beta;
	bool photbl = 1, leftvac = 0, tosmll = 0;
	int mult = 0, limit;
	
	
	specrel(index);                    //gets index of refraction
	specrel(vel, Beta);                //gets Beta
	alpha = acos(1/(index*Beta));      //defines Cherenkov angle
	
	//printf("tphoton called\n");
	
	if ((index*Beta) <= 1 ) 
  {
		printf("Huge error in Beta or index \n \n \n");
	}
	
	//printf("alpha = %f\n", alpha);

	
	for (int i = 0; i < 10; i++) 
	{
		//printf("\n loop is working %d  \n", i);
		
		zeta = 0.1*2*pie*i;
		projrad(vel, photvel, alpha, zeta, index, tosmll);   //projects the radial of the velocity
																												 //vector to the cone of Cherenkov Radiation
																												 //to get a propagation vector for the photon
		
		if (tosmll)                //if particle moving too slow to emit radiation, loop is broken
		{
			//printf("this should break loop\n");
			break;
		}
		
		for (int k = 0; k < 4; k++) 
		{
		 photpos[k] = pos[k] ;          //photon starts at same position as the particle when it is emitted
		 tmppos[k] = pos[k] ;          //assign placeholder array
		}
		time = 0;
		limit = 0 ;
		for (int v = 0; v < 50000; v++) //start photon movement
		{
			for (int m = 0; m < 3; m++) 
			{
				tmppos[m] += photvel[m]*0.01;    //update place holder for position of photon
			}
			time += 0.01;
			//printf("time  = %f\n", time);
			geom(tmppos, photbl, leftvac);   //check if it left detector
			limit = 0 ;
			while (!photbl) 
			{																		//if photon left the detector geom
				rfctp(tmppos, photpos, photvel);         //reflect photon off wall
				//insert function for transmission coefficients here
				for (int p = 0; p < 3; p++)             //update temporary position
				{
					tmppos[p] += photvel[p]*0.01;
				}
				time += 0.01;
			
				if (limit >= 2) 
				{
				printf("limit = %d (tphoton)\n", limit);
					for (int j = 0; j < 4; j++) 
					{
					printf("photpos[ %d ] = %f\n", j, photpos[j]);
					}
				for (int j = 0; j < 4; j++) 
				 {
				 printf("tmppos[ %d ] = %f\n", j, tmppos[j]);
				 }
				for (int j = 0; j < 4; j++) 
					{
					printf("photvel[ %d ] = %f\n", j, photvel[j]);
					}
				printf("\n");
				}
				geom(tmppos, photbl, leftvac);          //check if sent back into detector properly
				limit ++;
				if (limit >= 4)                        //if the loop has to keep repeating, something 				
				{
					printf("Error in reflection (tphoton)\n \n \n ");    //is wrong with the reflection function
					break;
				}
			}
			if (photbl) 
			{
				for (int h = 0; h < 3; h++)        
				{                                   
					photpos[h] = tmppos[h];         //if still in detector set position equal to placeholder
				}																		
			}
			else 
			{
				printf("Error in while(!photbl); photbl = %d\n \n \n", photbl);
				for (int i = 0; i<4; i++) 
				{
					printf("vel [ %d ] = %f\n", i, vel[i]);
				}
				printf("Beta = %f\n", Beta);
			}
			if (photpos[2] > 1)              
			{
				mult ++;																	//assign region for photo-multiplier to pick up photons
				/*for (int j = 0; j < 4; j++) 
				{
					printf("photpos[ %d ] = %f\n", j, photpos[j]);
				}*/
				//printf("photon inside detector = %d\n", photbl);
				break;
			}
			if (time >= 1500)                       //stop photon from bouncing around forever and never 
			{																				//hitting the detector
				printf("time went over 150 (tphoton)\n");
				break;
			}
		 }   //end movement
				 //printf("photon movement over; time = %f \n", time);
		/*if (i == 9) 
		{
			printf("%d photon hit photo-multiplier\n", mult);
			//printf("cone loop finished\n");
		}*/
	}       //end cone loop
	return mult;
	
}