% To add adsorption to the membranes
% A) Sorption reaction given by Langmuir type + Curvature sensing (No
% change in time scale) + Self Interaction
% Variables
% C_s: solute concentration in the bulk
% k_A: Coefficient of rate of adsorption 
% k_D: Coefficient of rate of desorption 

function R = Residual_Sorp_New(Aux,C,C1,C2,phi_k,d_phi_k,phi_1,dd_phi_k)
global Invariant ini_Val Material Constraints dt

if isfield(Material,'k_A') && isfield(Material,'k_D') && isfield(Constraints,'C_s')
    k_A = Material.k_A;
    k_D = Material.k_D;
    C_s = Constraints.C_s;
else
    k_A = 1;
    k_D = 0;
    C_s = 1;
end
    
if isfield(Material,'Lmd')
    Lmd = Material.Lmd;
else
    Lmd = 0;
end

% To limit the value of phi between 0 and 1
if isfield(Material,'alpha')
    alpha = Material.alpha;
else
    alpha = 0;
end

if isfield(ini_Val,'iopen')
    iopen = ini_Val.iopen;
else
    iopen = 0;
end

R = zeros(ini_Val.ncontrol,1);
connect = Invariant.connect;


Chi = Material.Chi;
RT_a0 = Material.RT_a0;
C0 = Constraints.C0;

% To Compute R1, R2, R3 and R4
tmp1 = (C(:,1).* Aux.au);
tmp2 = (C(:,1)./ Aux.au);
tmp3 = ((C1(:,1).*Aux.au - C(:,1).*Aux.dau)./(Aux.au).^2);

% To bound the phi between 0 and 1
tmp4 = alpha./phi_k./phi_k + alpha./(Material.phi_max-phi_k).^2;

% Langmuir type sorption
tmp8 = (k_A*C_s*(Material.phi_max - phi_k) - k_D*phi_k);
% Curvature sensing and generation to sorption
tmp9 = RT_a0.*log((k_A*C_s*(Material.phi_max-phi_k)./(k_D*phi_k)));
tmp10 = (Material.k0*Constraints.C0*Aux.HmC0-Chi*phi_k)./tmp9;

if (nargin == 8)
    L_phi_m = (tmp3.* d_phi_k + tmp2.*dd_phi_k)./tmp1;
    aux5 = tmp3.*phi_k.*L_phi_m.*Lmd .*Invariant.gweight*dt;
    aux6 = tmp2.*d_phi_k.*L_phi_m.*Lmd .*Invariant.gweight*dt;
    aux7 = tmp2.*phi_k.*L_phi_m.*Lmd .*Invariant.gweight*dt;
    aux8 = -tmp1.*tmp8.*Material.eta_a0.*Lmd.* L_phi_m./tmp9 .*Invariant.gweight*dt;
end

aux1 = tmp1.* (Material.eta_a0*phi_k) .* Invariant.gweight;

aux2 = - tmp1.* (Material.eta_a0*phi_1) .* Invariant.gweight;

aux3_a = tmp3.*(Material.k0*C0.*phi_k.*Aux.Hu).*Invariant.gweight*dt;
aux3_b = tmp2.*(Material.k0*C0.*d_phi_k.*Aux.Hu).*Invariant.gweight*dt;
aux3_c = tmp2.*(Material.k0*C0.*phi_k.*Aux.Hu).*Invariant.gweight*dt;

aux4 =  tmp2.* (RT_a0 * Material.phi_max./(Material.phi_max - phi_k) ...
    + (Material.k0*C0.^2 + Chi+ tmp4).*phi_k).* d_phi_k .*Invariant.gweight*dt;



aux9 = -tmp1.*tmp8.*Material.eta_a0.*Invariant.gweight*dt;

aux10 = -tmp1.*tmp8.*tmp10.*Material.eta_a0.*Invariant.gweight*dt;

for iel=1:size(connect,1)
    ig=(1:ini_Val.gaussno);
    igglob=(iel-1)*ini_Val.gaussno+ig;
    R1_loc = (aux1(igglob)')*Invariant.Nshape(igglob,:);
    R2_loc = (aux2(igglob)')*Invariant.Nshape(igglob,:);
    R3_loc_a = (aux3_a(igglob)')*Invariant.dNsh(igglob,:);
    R3_loc_b = (aux3_b(igglob)')*Invariant.dNsh(igglob,:);
    R3_loc_c = (aux3_c(igglob)')*Invariant.ddNsh(igglob,:);
    R4_loc = (aux4(igglob)')*Invariant.dNsh(igglob,:);
    R_loc = R1_loc + R2_loc + R3_loc_a +R3_loc_b+ R3_loc_c + R4_loc;
    if (nargin == 8)
        R5_loc = (aux5(igglob)')*Invariant.dNsh(igglob,:);
        R6_loc = (aux6(igglob)')*Invariant.dNsh(igglob,:);
        R7_loc = (aux7(igglob)')*Invariant.ddNsh(igglob,:);
        R8_loc = (aux8(igglob)')*Invariant.Nshape(igglob,:);
        R_loc = R_loc + R5_loc + R6_loc + R7_loc + R8_loc;
    end
    R9_loc = (aux9(igglob)')*Invariant.Nshape(igglob,:);
    R10_loc = (aux10(igglob)')*Invariant.Nshape(igglob,:);
    R_loc = R_loc + R9_loc + R10_loc;
    R(connect(iel,:)) = ...
        R(connect(iel,:)) + R_loc(:);
end

