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


% reaction in a spherical catalyst particle

function p2_main_unscaled

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_Nr = 1:length(Nr_vec)
    for ct = 1:length(R_vec)

        Nr = Nr_vec(ct_Nr);

        R = R_vec(ct);

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

        dr(ct) = R/(Nr-1);
        for i=1:Nr 
            r(i,ct) = 0 + (i-1)*dr(ct);
        end

        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,dr(ct),k_mass(ct));
        
        options=optimset('display','iter','tolfun',1e-14,'tolx',1e-14,'MaxFunEvals',100000,'MaxIter',1000);
        var_0 = S_guess(:,ct);
        S_vec(1:Nr,ct) = fsolve(@fin_diff_eqn,var_0,options,Nr,dr(ct),k_mass(ct));

        resid(:,i)=fin_diff_eqn(S_vec(:,ct),Nr,dr(ct),k_mass(ct));
        resid(1)
        resid(end)
        %Ds_gel*(S_vec(Nr,ct) - S_vec(Nr-1,ct))/dr(ct) - (k_mass(ct)*(S_bulk - S_vec(Nr,ct)))
    end

    max(S_vec);

    for i=1:length(R_vec)
        cons_rate(i,ct_Nr) = k_mass(i)*(S_bulk - S_vec(end,i))*4*3.1415*R_vec(i)^2;    % moles/sec
        cons_rate_in(i,ct_Nr) = Ds_gel*(S_vec(end,i) - S_vec(end-1,i))/dr(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,ct_Nr) = cons_rate(i,ct_Nr)/mass_E(i)*1e6;     % moles/mgE-sec
        mole_per_kgE_in(i,ct_Nr) = cons_rate_in(i,ct_Nr)/mass_E(i)*1e6;     % moles/mgE-sec

        %disp(['Consumption Rate (Outer) for R = ',num2str(R_vec(i)),' m:  ',num2str(mole_per_kgE(i,ct_Nr)),' moles S per kg E']);
        disp('  ');
        disp(['Results for R = ',num2str(R_vec(i)),' meters:']);
        disp(['[S] at the bead surface:  ',num2str(S_vec(end,i)),' moles/m^3']);
        disp(['Consumption Rate:  ',num2str(cons_rate_in(i,ct_Nr)),' moles/sec']);
        disp(['Consumption Rate per kg of Enzyme:  ',num2str(mole_per_kgE_in(i,ct_Nr)),' moles/sec S per kg E']);

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

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

figure;
plot(r(:,1)/R_vec(1),S_vec(:,1),'-',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(r_scaled,S_vec)

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

global Ds_gel S_bulk rate_0 Km

S = var;

for i=1:Nr
    
    r(i) = 0 + (i-1)*dr;
    
    if(i == 1)
        resid(i) = (S(i) - 4/3*S(i+1) + 1/3*S(i+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 - (k_mass/Ds_gel*(S_bulk - S(i)));
        resid(i) = (Ds_gel*(S(i) - S(i-1))/dr - (k_mass*(S_bulk - S(i))))*1e9;
    else
        resid(i) = r(i)/dr^2*(S(i+1) - 2*S(i) + S(i-1))...
            + 2/(2*dr)*(S(i+1) - S(i-1))...
            - r(i)*rate_0/Ds_gel*S(i)/(Km+S(i));
    end
end

return;

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

global Ds_gel S_bulk rate_0 Km

S = var;

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

return;