% 10.34 - Fall 2006
% Hw Set #8 - Problem #2
% Rob Ashcraft - Oct. 25, 2006


% reaction in a spherical catalyst particle

function p2_main_nonuniform

clear;

global Ds_gel S_bulk rate_0 Km

Nr_vec = [101];
    
Sh = 10;

R_vec = [0.0005, 0.001, 0.005];

%R_vec = logspace(log10(2e-4), log10(0.02), 21)'
%R_vec = linspace(2e-4, 0.02, 21)';


S_bulk = 1000;  % moles/m^3

rhoE = 0.01*1e6;  % mgE/m^3 gel
Ds_gel = 1e-10;  % m^2/s
Ds_liq = 7e-10;  % m^2/s

Vm = 200/60/1e6;  % mole/s-mgE
Km = 200;  % mole/m^3

rate_0 = Vm*rhoE;  % mole/m^3-s

S_vec = zeros(max(Nr_vec),length(R_vec));


for ct = 1:length(R_vec)
    
    Nr = Nr_vec(1);
    
    R = R_vec(ct);

    k_mass(ct) = Sh*Ds_liq/(2*R);    % m/s

    V_total(ct) = 4/3*3.1415*R^3;
    dV = V_total(ct)/(Nr-1);
    r(1,ct) = 0;
    for i=2:Nr 
        r(i,ct) = (3*dV/(4*3.1415) + r(i-1,ct)^3)^(1/3);  % not scaled to problem
        dr(i-1,ct) = r(i,ct) - r(i-1,ct);   % not scaled to problem
    end
%     for i=1:Nr
%         if(i == 1)
%             r_temp(i,ct) = 0;
%             dr_temp(i,ct) = 0.01;
%         else
%             r_temp(i,ct) = r_temp(i-1,ct) + dr_temp(i-1,ct);  % not scaled to problem
%             dr_temp(i,ct) = dr_temp(i-1,ct)/1.01;   % not scaled to problem
%         end
%     end
% 
%     r(:,ct) = r_temp(:,ct)*(R/r_temp(end,ct));  % scale it to the appropriate length in the problem
%     dr(:,ct) = dr_temp(:,ct)*(R/r_temp(end,ct));
    
    dr(1,ct)/dr(end,ct)
    
    options=optimset('display','off','tolfun',1e-12,'tolx',1e-12,'MaxFunEvals',100000,'MaxIter',1000);
        
    if(ct == 1)
        var_guess = linspace(0,1000,Nr)';
    else
        var_guess = S_guess(:,ct-1);
    end
    S_guess(1:Nr,ct) = fsolve(@fin_diff_eqn_const,var_guess,options,Nr,r(:,ct),dr(:,ct),k_mass(ct));

    options=optimset('display','off','tolfun',1e-12,'tolx',1e-12,'MaxFunEvals',100000,'MaxIter',1000);
    var_0 = S_guess(:,ct);
    S_vec(1:Nr,ct) = fsolve(@fin_diff_eqn,var_0,options,Nr,r(:,ct),dr(:,ct),k_mass(ct));
    resid(:,i)=fin_diff_eqn(S_vec(:,ct),Nr,r(:,ct),dr(:,ct),k_mass(ct));
    resid(1)
    resid(end)
end

max(S_vec)

for i=1:length(R_vec)
    cons_rate(i) = k_mass(i)*(S_bulk - S_vec(end,i))*4*3.1415*R_vec(i)^2;    % moles/sec
    cons_rate_in(i) = Ds_gel*(S_vec(end,i) - S_vec(end-1,i))/dr(end-1,i)*4*3.1415*R_vec(i)^2;    % moles/sec
    mass_E(i) = 4/3*3.1415*R_vec(i)^3 * rhoE;  % mgE per bead
    mole_per_kgE(i) = cons_rate(i)/mass_E(i)*1e6;     % moles/mgE-sec
    mole_per_kgE_in(i) = cons_rate_in(i)/mass_E(i)*1e6;     % moles/mgE-sec
    
    disp(['Consumption Rate Out for R = ',num2str(R_vec(i)),' m:  ',num2str(mole_per_kgE(i)),' moles S per kg E']);
    disp(['Consumption Rate In for R = ',num2str(R_vec(i)),' m:  ',num2str(mole_per_kgE_in(i)),' moles S per kg E']);

    r_scaled(:,i) = r(:,i)/R_vec(i);
end
k_mass;

%plot(r,dr);

figure; plot(R_vec,mole_per_kgE,'o');

figure;
plot(r(:,1)/R_vec(1),S_vec(:,1),'-o',r(:,2)/R_vec(2),S_vec(:,2),'-.',r(:,3)/R_vec(3),S_vec(:,3),'--');
xlabel('r/R  (scaled radius)'); ylabel('Concentration  (mole/m^3)');
legend('R = 0.5 mm','R = 1.0 mm','R = 5.0 mm');

% figure;
% plot(S_vec)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function resid = fin_diff_eqn(var,Nr,r,dr,k_mass)

global Ds_gel S_bulk rate_0 Km

S = var;

for i=1:Nr
    
    if(i == 1)
        dr12 = dr(i);
        dr13 = dr(i) + dr(i+1);
        resid(i) = (S(i)*(dr13^2 - dr12^2) - (S(i+1)*dr13^2 - S(i+2)*dr12^2))*1e9;
        
    elseif(i == Nr)  % BC @ r=R: Dgel*dS/dr = k_mass*(S_bulk - S(r=R))
        resid(i) = ((S(i) - S(i-1))/dr(i-1) - (k_mass/Ds_gel*(S_bulk - S(i))))*1e9;
        %resid(i) = (S(i) - 1000);
        
    else
%          resid(i) = r(i)/(dr(i-1)*dr(i))*(S(i+1) - 2*S(i) + S(i-1))...
%              + 2/(dr(i-1)+dr(i))*(S(i+1) - S(i-1))...
%              - r(i)*rate_0/Ds_gel*S(i)/(Km+S(i));
        resid(i) = r(i)*( (S(i+1) - S(i))/dr(i) - (S(i) - S(i-1))/dr(i-1) )/(dr(i)/2 + dr(i-1)/2)...
            + 2*( (S(i+1) - S(i))*dr(i-1)/dr(i) + (S(i) - S(i-1))*dr(i)/dr(i-1) )/(dr(i-1)+dr(i))...
            - r(i)*rate_0/Ds_gel*S(i)/(Km+S(i));
    end
end

return;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function resid = fin_diff_eqn_const(var,Nr,r,dr,k_mass)

global Ds_gel S_bulk rate_0 Km

S = var;

for i=1:Nr
    
    if(i == 1)
        resid(i) = S(i) - S(i+1);
        
    elseif(i == Nr)  % BC @ r=R: Dgel*dS/dr = k_mass*(S_bulk - S(r=R))
        resid(i) = (S(i) - 1000);
        
    else
%          resid(i) = r(i)/(dr(i-1)*dr(i))*(S(i+1) - 2*S(i) + S(i-1))...
%              + 2/(dr(i-1)+dr(i))*(S(i+1) - S(i-1))...
%              - r(i)*rate_0/Ds_gel*S(i)/(Km+S(i));
        resid(i) = r(i)*( (S(i+1) - S(i))/dr(i) - (S(i) - S(i-1))/dr(i-1) )/(dr(i)/2 + dr(i-1)/2)...
            + 2*( (S(i+1) - S(i))*dr(i-1)/dr(i) + (S(i) - S(i-1))*dr(i)/dr(i-1) )/(dr(i-1)+dr(i))...
            - r(i)*rate_0/Ds_gel*S(i)/(Km+S(i));
    end
end

return;