
function K = matrix_K_mod_Ener_Aug(Aux,C,C1,C2,phi)
% To be multiplied with dt later while computing the new phi


global Invariant ini_Val Material Constraints dt

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

if isfield(Material,'Lmd')
    Lmd = Material.Lmd;
else
    Lmd = 0;
end

% For the (alpha/2)*log(beta*phi)  term added to the energy
if isfield(Material,'alpha')
    alpha = Material.alpha;
else
    alpha = 0;
end


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

%
if (alpha >0)
    tmp4 = alpha./phi./phi + alpha./(Material.phi_max-phi).^2;
else
    tmp4=0;
end

% tmp4 = alpha*beta^2*sech(beta*(Material.phi_max-phi)).^2.*tanh(beta*(Material.phi_max-phi));
% tmp6 = -gamma*delta^2.*sech(delta*(phi+0.05)).^2.*tanh(delta*(phi+0.05));


K=zeros(ini_Val.ncontrol);
connect = Invariant.connect;
aux =  (C(:,1)./ Aux.au) .* (RT_a0 * Material.phi_max./(Material.phi_max - phi) ...
    + (Material.k0*C0.^2 + Chi + tmp4).*phi).* Invariant.gweight*dt;

%aux2 =  ((C1(:,1).*Aux.au - C(:,1).*Aux.dau)./(Aux.au).^2).^2 .* Lmd .* Invariant.gweight;


for iel=1:size(connect,1)
    ig=(1:ini_Val.gaussno);
    igglob=(iel-1)*ini_Val.gaussno+ig;
    K_loc=Invariant.dNsh(igglob,:)'*diag(aux(igglob))*Invariant.dNsh(igglob,:);
    K(connect(iel,:),connect(iel,:)) = ...
        K(connect(iel,:),connect(iel,:)) + K_loc;
end