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