% Energy Values
Energy = [1.00E-03	1.50E-03	2.00E-03	3.00E-03	4.00E-03	5.00E-03	6.00E-03	8.00E-03	1.00E-02	1.50E-02	2.00E-02	3.00E-02	4.00E-02	5.00E-02	6.00E-02	8.00E-02	1.00E-01	1.50E-01	2.00E-01	3.00E-01	4.00E-01	5.00E-01	6.00E-01	8.00E-01	1.00E+00
];
% This is the C.S Ratio associated with each energy above
CS_Ratio = [3.23602E-06	1.94541E-05	6.79E-05	0.000368682	0.001149019	0.002671901	0.00520463	0.014314115	0.030398117	0.110396361	0.245944822	0.556603774	0.762839248	0.86849711	0.921875	0.966951567	0.983071342	0.995273464	0.997787611	0.999152542	1	0.999793068	0.999888143	0.999872725	1
];



% Declaring three variables to store counts towards flux values
Uncollided_Flux = 0;
Collided_Flux = 0;
Undetected_Flux = 0;

% Calculation of Distance for Slab
MFP = 1;
Mass_Atten_Coeff = 7.07e-2;
Rho = 1;
Linear_Atten_Coeff = Mass_Atten_Coeff * Rho;
% x distance = MFP/Linear_Atten_Coeff
slab_distance = (MFP)/(Linear_Atten_Coeff);

% Number of trials
Trials = 50000;

for k=1:Trials
   
% Reset Energy_Gamma to 1 for new trial
Energy_Gamma = 1.0;

% Calculation of distance travelled
RAND_Distance = rand;
Distance_Travel = (1/Linear_Atten_Coeff)*log(RAND_Distance);

% Calculation of Angle
RAND_Angle = rand;
Angle = RAND_Angle * 2 * 3.141592;

% Position of X-Coordinate
X_Coord = (Distance_Travel)*(cos(Angle));
escape = 0;
    
% Check when gamma is in the slab
    while escape == 0
        escape = 1;

        if X_Coord > 0 && X_Coord <= slab_distance;
            % Chance of CS effect given energy value
            Interp_CS_Ratio = interp1(Energy,CS_Ratio,Energy_Gamma);
            RAND_CS_Chance = rand;

            % Accounting for absorption if photo-electric effect
            if RAND_CS_Chance > Interp_CS_Ratio
                Undetected_Flux = Undetected_Flux + 1;

            % New energy value if CS
            else
                Energy_Gamma = Energy_Gamma/(1+(Energy_Gamma/0.511)*(1-cos(Angle)));

                if Energy_Gamma <= 0.001
                    Undetected_Flux = Undetected_Flux + 1;

                else
                    RAND_Angle=rand;
                    Angle = RAND_Angle*2*pi;

                    RAND_Distance=rand;
                    X_Coord = X_Coord + Distance_Travel*(cos(Angle));
                    escape = 0;

                end
            end
        end
    end
        
% Check if gamma is in the wrong direction
    if X_Coord <= 0    
        Undetected_Flux = Undetected_Flux + 1;
        
% Check if gamma is uncollided
    elseif X_Coord > slab_distance && Energy_Gamma == 1;
        Uncollided_Flux = Uncollided_Flux + 1;
        
    
% Check if gamma is collided    
    elseif X_Coord > slab_distance && Energy_Gamma < 1
        Collided_Flux = Collided_Flux + 1;
    
    end
        
end    

Buildup = (Uncollided_Flux + Collided_Flux)/Uncollided_Flux;

disp ( Buildup )