%calculates the total hits and the solid angle % %draw Simple Cylinder and checks for hits for an individual trajectory r=1; % inner radius of the beaker r1=1.001; % outer radius of the beaker Hits=0; i=1; n=100000; while i<=n %getting random points theta1=2*180*rand(); % problem theta is never going to be zero or 360 z1=(randi(10000)+rand())*(-1)^(randi(2)); R1=r+(r1-r)*(rand())^.5; %creates random uniform distribution x1=R1*cosd(theta1); y1=R1*sind(theta1); %getting random trajectory R2=10; % R2 is the length of the trajectory theta= acosd(1-2*rand()); %uniform theta distribution in degrees. REMEBER this is in FLOATING POINT Phi= 2*180*rand(); %uniform phi distribution in degrees. REMEBER this is in FLOATING POINT. also Phi will never be zero or 360 x2p= R2*sind(theta)*cosd(Phi); % measured in the source coordinate y2p= R2*sind(theta)*sind(Phi); %measured in the source coordinate z2p= R2*cosd(theta); %measured in the source coordinate x2=x2p+x1; % measured in the cylinder coordinate y2=y2p+y1; % measured in the cylinder coordinate z2=z2p+z1; % measured in the cylinder coordinate %plotting % %plotting point % plot3([x1,x2],[y1,y2],[z1,z2]); % hold on % % %plots a cylinder in 3d which can be rotated in 3d % cylinder(1); % a cylinder with radius 1 % colormap([1,1,1]); % hold on % daspect([1 1 1]) %streches the sphere % axis ([-3 3 -3 3]) % xlabel('x-axis') % ylabel('y-axis') % zlabel('z-axis') % grid on % rotate3d on %counting hits Condition=0; %initialization %checking for exception conditions where Phi blows up if (ceil(Phi*10000)/10000== 180.00001/2) && (ceil(theta*10000)/10000 == 180.00001/2 ) % comparing up to four sig fig xaa=x1; xSource=xaa-x1; % always zero %check for the absolute value of x if abs(x1)== r %tangent condition so does not matter which direction the %trajectory goes. But we are not going to count this as a hit. Condition=0; elseif abs(x1)0 && ySourceCoa1>0 PhiPrime1=aPhiPrime1; elseif xaa<0 && ySourceCoa1>0 PhiPrime1=aPhiPrime1 +180; elseif xaa<0 && ySourceCoa1<0 PhiPrime1=aPhiPrime1+180; else PhiPrime1=aPhiPrime1+360; end %for aPhiPrime2 if xaa>=0 && ySourceCoa2>0 PhiPrime2=aPhiPrime2; elseif xaa<0 && ySourceCoa2>0 PhiPrime2=aPhiPrime2+180; elseif xaa<=0 && ySourceCoa2<0 PhiPrime2=aPhiPrime2+180; else PhiPrime2=aPhiPrime2+360; end % rounding to four sig digits Phi=ceil(Phi*10000)/10000; PhiPrime1=ceil(PhiPrime1*10000)/10000; PhiPrime2=ceil(PhiPrime2*10000)/10000; %check if the PhiPrime == Phi if (PhiPrime1==Phi) || (PhiPrime2==Phi) Condition=3; end else % it does not intersect the circle at all Condition=0; end elseif(Phi== (180.0001/2))&& (theta== 0 ) %never touches the detector Condition=0; else %getting x and y projection %r= radius of the cylinder a=1+(tand(Phi))^2; b=2*tand(Phi)*(y1-x1*tand(Phi)); c=(x1^2)*(tand(Phi))^2-2*x1*y1*tand(Phi)+(y1^2)-r^2; xa= (-b+((b^2)-4*a*c)^.5)/(2*a); xb= (-b-((b^2)-4*a*c)^.5)/(2*a); if (isreal(xa)||isreal(xb))==1 Condition=1; %Checking for condition 2 (if it hits the detector) xSourceCo1=xa-x1; %finding new value of x prime in the source coordinate xSourceCo2=xb-x1; zSourceCo1=(xSourceCo1 * cosd(theta))/(sind(theta)*cosd(Phi)); % finding value of z prime in the source coordinate zSourceCo2=(xSourceCo2 * cosd(theta))/(sind(theta)*cosd(Phi)); zCylinderCo1=zSourceCo1+z1; %getting the value of z prime in the cylinder coordinate zCylinderCo2=zSourceCo2+z1; if (zCylinderCo1>-10000 && zCylinderCo1<10000) || (zCylinderCo2>-10000 && zCylinderCo2<10000) %checking if the z value found is greater then the height of the cylinder Condition=2; %checking for condition 3 (direction of the trajectory) %getting all the values of y ya1=-((r^2)-(xa^2))^.5; % finds the y when the trajectory intersects the cylinder ya2=-ya1; yb1=((r^2)-(xb^2))^.5; yb2=-yb1; if (isreal(ya1)||isreal(yb1))==0 fprintf(' Phi = %f theta = %f x1 = %f y1 = %f z1= %f xa= %f xb= %f ',Phi, theta, x1,y1,z1, xa,xb) end ySourceCo11=ya1-y1; % converting into source coordinate ySourceCo12=ya2-y1; ySourceCo21=yb1-y1; ySourceCo22=yb2-y1; %getting all the values of PhiPrime PhiPrime11= atand(ySourceCo11/xSourceCo1); %finds the Phi in the source coordinate PhiPrime12= atand(ySourceCo12/xSourceCo1); PhiPrime21= atand(ySourceCo21/xSourceCo2); PhiPrime22= atand(ySourceCo22/xSourceCo2); %checking for the diffrent quadrent condition %for PhiPrime 11 if xSourceCo1>=0 && ySourceCo11>=0 PhiPrime1=PhiPrime11; elseif xSourceCo1<0 && ySourceCo11>=0 PhiPrime1=PhiPrime11+180; elseif xSourceCo1<0 && ySourceCo11<0 PhiPrime1=PhiPrime11+180; else PhiPrime1=PhiPrime11+360; end %for PhiPrime 12 if xSourceCo1>=0 && ySourceCo12>=0 PhiPrime2=PhiPrime12; elseif xSourceCo1<0 && ySourceCo12>=0 PhiPrime2=PhiPrime12+180; elseif xSourceCo1<0 && ySourceCo12<0 PhiPrime2=PhiPrime12+180; else PhiPrime2=PhiPrime12+360; end %for PhiPrime 21 if xSourceCo2>=0 && ySourceCo21>=0 PhiPrime3=PhiPrime21; elseif xSourceCo2<0 && ySourceCo21>=0 PhiPrime3=PhiPrime21+180; elseif xSourceCo2<0 && ySourceCo21<0 PhiPrime3=PhiPrime21+180; else PhiPrime3=PhiPrime21+360; end %for PhiPrime 22 if xSourceCo2>=0 && ySourceCo22>=0 PhiPrime4=PhiPrime22; elseif xSourceCo2<0 && ySourceCo22>=0 PhiPrime4=PhiPrime22+180; elseif xSourceCo2<0 && ySourceCo22<0 PhiPrime4=PhiPrime22+180; else PhiPrime4=PhiPrime22+360; end %fprintf('PhiPrime1 = %f, PhiPrime2= %f, PhiPrime3= %f, and PhiPrime4= %f/n Phi =%f',PhiPrime1,PhiPrime2,PhiPrime3,PhiPrime4, Phi); Phi=ceil(Phi*10000)/10000; PhiPrime1=ceil(PhiPrime1*10000)/10000; PhiPrime2=ceil(PhiPrime2*10000)/10000; PhiPrime3=ceil(PhiPrime3*10000)/10000; PhiPrime4=ceil(PhiPrime4*10000)/10000; if (PhiPrime1==Phi) || (PhiPrime2==Phi) || (PhiPrime3==Phi) || (PhiPrime4==Phi) Condition=3; end end end end %Checking if the ray satisfied both conditions if Condition ==3 Hits= Hits+1; end i=i+1; end TotalHits= Hits; SolidAngle=TotalHits/n; fprintf('total hits = %d and Solid angle = %d\n',TotalHits, SolidAngle)