clear
clc
close all

% Read in the data
a = csvread('Run 5 Staack EP.csv',9,0);
KE = a(:,9);
x =  a(:,13);
y =  a(:,14);
z =  a(:,15);

%Initialize variables
counter = 0;
icounter = 0;
total = size(z,1);
ke_sum = 0;
kcount = 1;
mass_flow_rate = 2.2612e21;
Kb = 1.38e-23; % J/K
q = -1.602e-19; % C
Me = 9.109e-31; % kg
Mi = 6.633e-26; % kg
Na = 6.022141e23;       % particles/mol
AL_density = 2700;      % kg/m^3
SI_density = 2330;      % kg/m^3
AL_Molar = 0.0269815/Na;% kg/particles
SI_Molar = 0.06008/Na;  % kg/particles
Area = 23.9; % m^3
Impact_KE = 0;
SI_Yield = 32.84; % sputter energy threshold
AL_Yield = 21.5;  % sputter energy threshold
AL_atom_to_meters = AL_Molar/(AL_density*Area); % converts # of atoms to meters eroded
SI_atom_to_meters = SI_Molar/(SI_density*Area); % converts # of atoms to meters eroded
Time = 60*60*24*30; % converts seconds to months
Scale_up = mass_flow_rate/total; % scales simulated ions to total mass flow rate

% Determine what the radius is from x and y
for g=1:total
    if x(g)>=0 && y(g)>=0 
        rc(g) = sqrt(x(g)^2+y(g)^2);
    else
        rc(g) = -1*sqrt(x(g)^2+y(g)^2);
    end
end

r = rc(:);

% Binning
for i=1:total
    if z(i)<=1e-5
        counter = counter+1;
        if abs(rc(i))<=5
            icounter = icounter + 1;
            Impact_KE(icounter) = KE(i)/1.602176E-19;
        else 
            icounter = icounter;
        end
    else
        counter = counter;
    end
end



for j=1:icounter
    Ve(j) = sqrt(2*(Impact_KE(icounter)/6.24e18)/Me);
    Te(j) = 1/3*(Me*(Ve(j).^2)/Kb);
    VBohm = sqrt(Kb*Te(j)/Mi);
    Potential(j)=-Mi/q*((2*Impact_KE(icounter)/(6.24e18*Mi))-VBohm^2);
    Ji(j)=(icounter/25.4)*VBohm;
end



% backfield = (counter/size(z,1))*100
% impact = (icounter/size(z,1))*100


% for k=1:size(Impact_KE',1)  % Erosion of Silicon
%     if Impact_KE(k)<32.84 % Just above the erosion threshold for Silicon
%         YieldSI(k) = 0;
%     elseif Impact_KE(k)<37 && Impact_KE(k)>=32.84 % Below the erosion threshold so no erosion
%         YieldSI(k) = 1e-5;
%     else
%         YieldSI(k) = Sputter_yieldSI(round(Impact_KE(k)));
%     end
% end 
% 
% for k=1:size(Impact_KE',1)  % Erosion of Silicon
%     if Impact_KE(k)<21.5 % Just above the erosion threshold for Silicon
%         YieldAL(k) = 0;
%     elseif Impact_KE(k)<27 && Impact_KE(k)>=21.55 % Below the erosion threshold so no erosion
%         YieldAL(k) = 1e-5;
%     else
%         YieldAL(k) = Sputter_yield(round(Impact_KE(k)));
%     end
% end 

for k=1:size(Impact_KE',1) 
    if Impact_KE(k)<=21.5497
        YieldAL2(k)=0;
    else    
%         SnAL = 0.5*log(1+1.2288*6.28194e4)/(6.28194e4+(0.1728*sqrt(6.28194e4))+(0.008*(6.28194e4^0.1504)));
%         YieldAL2(k)=4.7928*SnAL*(((Impact_KE(k)/21.5497)-1)^2.0810)/(0.4713+((Impact_KE(k)/21.5497)-1)^2.0810);
        YieldAL2(k)= 2.7E-7*(Impact_KE(k)^3)-1.57E-5*(Impact_KE(k)^2)+0.00035*Impact_KE(k)-0.003;
    end
end

for k=1:size(Impact_KE',1) 
    if Impact_KE(k)<=32.838
        YieldSI2(k)=0;
    else    
%         SnSI = 0.5*log(1+1.2288*6.67979e4)/(6.67979e4+(0.1728*sqrt(6.67979e4))+(0.008*(6.67979e4^0.1504)));
%         YieldSI2(k)=3.2299*SnSI*(((Impact_KE(k)/32.838)-1)^1.5284)/(0.2770+((Impact_KE(k)/32.838)-1)^1.5284);
        YieldSI2(k)= -6.44e-9*(Impact_KE(k)^4)+1.03e-6*(Impact_KE(k)^3)-5.5e-5*Impact_KE(k)^2+0.00121*Impact_KE(k)-0.00945;
    end
end

Total_AL_Erosion = sum(YieldAL2); % in # of AL atoms eroded
Total_SI_Erosion = sum(YieldSI2); % in # of SI atoms eroded

Total_erosionAL = 1e9 * Time * Total_AL_Erosion * Scale_up * AL_atom_to_meters % nm/mo
Total_erosionSI = 1e9 * Time * Total_SI_Erosion * Scale_up * SI_atom_to_meters % nm/mo

mean(YieldAL2)
