%using 241Am source of known activity to test if the code works
% theoritical calculation is in cm
% detector diameter is 3.89 cm and the height is 4.72cm
% the base of the cylinder is at (0,0,0)
%calculates the total hits and the solid angle
% %draw Simple Cylinder and checks for hits for an individual trajectory
r=1.9500;    % inner radius of the beaker in cm
thickness=1.2;   % thickness of the sample

Hits=0;             % keeps track of number of gamma that did not scattered in the detector
i=1;
n=100000;
%totalSc=0;
% totalSj=0;
totalDistance=0;
PotentialHits=0;                % keeps track of how many times it fulfilled condition 3
k=0;                % keeps track of number of gammas that did not scattered in the sample
NotCompton=0;
NotPhotoelectric=0;
Indetector=0;


while i<=n
    %getting source points
    theta1=2*180*rand();            % problem theta is never going to be zero or 360
    z1=(randi([140,260],1)+rand())/100;
    x1=1.37886;
    y1=1.37886;
    
    %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== 90.0000) && (ceil(theta*10000)/10000 == 90.0000 ) % comparing up to four sig fig
        if z1>=1.78 && z1 <=6.50
            xa=x1;
            xb=x1;
            xSource=xa-x1;     % always zero
            %check for the absolute value of x
            if abs(xa)== 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(xa)<r
                % intersects the circle at two points
                Condition=1;
                %since xaa=x1 the y of the source and the cylinder coordinate is
                %the same
                ya1=(r^2-(xa)^2)^.5;
                yb1=-(r^2-(xb)^2)^.5;
                
                ySourceCoa1=ya1-y1;             % converting into source coordinate
                ySourceCoa2=yb1-y1;
                
                %check for direction
                
                % for this condition the aPhiPrime1 and aPhiPrime3 will be positive
                % 90 no matter what.
                PhiPrime1= atand(ySourceCoa1/xSource)*1.0000001; %finds the Phi in the source coordinate and multipling by 1.0000001 to make it a floating point number
                PhiPrime3= atand(ySourceCoa2/xSource)*1.0000001;
                
                %for PhiPrime1
                if xa>=0 && ySourceCoa1>=0
                    PhiPrime1=PhiPrime1+0;
                elseif xa<0 && ySourceCoa1>=0
                    PhiPrime1=PhiPrime1 +180;
                elseif xa<0 && ySourceCoa1<0
                    PhiPrime1=PhiPrime1+180;
                else
                    PhiPrime1=PhiPrime1+360;
                end
                
                %for PhiPrime3
                if xa>=0 && ySourceCoa2>=0
                    PhiPrime3=PhiPrime3+0;
                elseif xa<0 && ySourceCoa2>=0
                    PhiPrime3=PhiPrime3+180;
                elseif xa<0 && ySourceCoa2<0
                    PhiPrime3=PhiPrime3+180;
                else
                    PhiPrime3=PhiPrime3+360;
                end
                
                % rounding to four sig digits
                Phi=ceil(Phi*10000)/10000;
                PhiPrime1=ceil(PhiPrime1*10000)/10000;
                PhiPrime3=ceil(PhiPrime3*10000)/10000;
                
                %check if the PhiPrime == Phi
                if (PhiPrime1==Phi) && (PhiPrime3==Phi)
                    Condition=3;
                    PhiPrime2=0;
                    PhiPrime4=0;
                end
            else
                % it does not intersect the circle at all
                Condition=0;
            end
        end
        
    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;
            if xSourceCo1==0 
                xSourceCo1=xa;
            end 
            if xSourceCo2==0 
               xSourceCo2=xb;
            end 
            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>=0 && zCylinderCo1<= 4.72) || (zCylinderCo2>=0 && zCylinderCo2<=4.72)       %checking if the z value found is greater then the height of the cylinder in cm
                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;
                
                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 1
                
                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 2
                
                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 3
                
                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 4
                
                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
                
                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;
                
                %checking if the trajectory is in the right direction
                if (PhiPrime1==Phi) || (PhiPrime2==Phi) || (PhiPrime3==Phi) || (PhiPrime4==Phi)
                    Condition=3;
                    
                    
                end
            end
        end
    end
    % finding the two right PhiPrime to get the values of y and
    % x
    % there is two correct value for PhiPrime=Phi. One for xa and
    % other for xb
    if Condition==3
        PotentialHits=PotentialHits+1;                            % keeps track of how many times it fulfilled condition 3
        
        % PhiPrime1==PhiPrime3 and PhiPrime2==PhiPrime4 so
        % checking for anyone of them is fine
        if Phi==PhiPrime1
            yf1=ya1;
            zf1=zCylinderCo1;
        end
        if Phi==PhiPrime2
            yf1=yb1;
            zf1=zCylinderCo1;
        end
        if Phi==PhiPrime3
            yf2=yb1;
            zf2=zCylinderCo2;
        end
        if Phi==PhiPrime4
            yf2=yb2;
            zf2=zCylinderCo2;
        end
        
            %attenuation in the detector
            % finding the distance travelled by the trajectory in the
            % detector
            
            distance =((xa-xb)^2 + (yf1-yf2)^2 + (zf1-zf2)^2)^.5; % distance travelled inside the detector
            totalDistance=totalDistance+1;
            %attenuation due to compton scattering
            Sc=-(log(1-rand()))*1.4973;       %4.12=1/attenuation length(units=cm). For detector
            Si=-(log(1-rand()))*0.1057;  % 243.76 = 1/attenuation length(units=cm). For detector
          
            if Sc<Si %Photoelectric effect
                
                if distance < Sc && distance > Si
                    Photo=1;
                end 
                if distance < Sc && distance <Si
                    Photo=0;
                    Compton=0;
                end 
                if distance > Sc 
                    Photo=1;
                    Compton=0;
               
            
            
            
                k= k+1; % keeps track of the number of trajectory that did not get scattered in the detector
                %attenuation due to photoelectric effect
                
                %totalSi=totalSi+Si;         % to calculate the average Si
                %totalDistance=totalDistance+distance;
                
                %Checking if the ray satisfied both conditions
                if Si<distance
                    Hits= Hits+1;           %number of gammas that did not scattered in the detector
                end
            end
    end
    i=i+1;
end
PotentialHits
NotCompton
NotPhotoelectric
k
Hits
AverageDistance=totalDistance/n
%AvgtotalSc=totalSc/Indetector
% SamplePerCompton=((PotentialHits-NotCompton)/PotentialHits)*100
% SamplePerPhotoelectric=((NotCompton-NotPhotoelectric)/PotentialHits)*100
% DetectorPerCompton=((NotPhotoelectric-k)/NotPhotoelectric)*100
% DetectorPerPhotoelectric=(Hits/NotPhotoelectric)*100
%meanSi=totalSi/k;
%meanSj=totalSj/j;
%meanSc=totalSc/c;
%AvegDistance=totalDistance/k;
Efficiency=Hits/n
SolidAngle=PotentialHits/n
%fprintf('____________ \ntotal successful hits = %d \nTotal efficiency = %f \nAverage distance travelled by the trajectory inside the detector = %f\nMean scatter length in detector = %f\nMean Scatter length in sample = %f\nPercentage of gamma rays absorbed in the sample= %d\nPercentage of gamma rays absorbed the detector = %d\nSolid angle = %f\n',Hits, Efficiency,AvegDistance,meanSi,meanSj,((j-NotPhotoelectric)/j)*100,(Hits/NotPhotoelectric)*100,j/n)