/* * 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); } }