function [tE,dE_dCP] = Curvature_functional_Reg(Aux,C)
global Invariant ini_Val Material Constraints


p = ini_Val.p;
gaussno = ini_Val.gaussno;

if isfield(Material,'Rg_a') %Regularization of zero radius
    Rg_a = Material.Rg_a;
else
    Rg_a = 0;
end



tE = 2*pi*sum(Aux.HmC02.*Aux.ds.*Invariant.gweight);
tE = tE*Material.k0/2;

if(Rg_a~=0)    
    tmp_lb= 0.15;
    tmp_ub= 0.9;
    tmp_sh = 100; %specifies the sharpness of the reg potential in terms of gauss points
    
    tmp_cf = 0.5; %specifies the cutoff radius
    tmp_fn = 0.5*(tanh(tmp_sh*(Invariant.gaussp-tmp_lb)) ...
        - tanh(tmp_sh*(Invariant.gaussp - tmp_ub)))';

    tE = tE + 2*pi*Rg_a*sum(1./((C(:,1) - tmp_cf).^2).*Aux.ds.*tmp_fn.*Invariant.gweight);
end

if nargout == 2
    %Derivatives of E (local)
    temp=2*Aux.MHmC0.*Aux.Mau.*Aux.MC1;
    dE_dr_l = ((temp.*Aux.dH_dr_l)+(Aux.MHmC02.*(Aux.MC1.*Aux.da_dr_l ...
        +Aux.Mau.*Invariant.Nshape))).*Aux.Mgweight;
    dE_dz_l = ((temp.*Aux.dH_dz_l)+(Aux.MHmC02.*(Aux.MC1.*Aux.da_dz_l))) ...
        .*Aux.Mgweight;
    
    if(Rg_a~=0)
        tmp_pfn = tmp_fn*ones(1,p+1);
        tmp_MC2 = (C(:,1)-tmp_cf).^2*ones(1,p+1);
        tmp_MC3 = (C(:,1)-tmp_cf).^3*ones(1,p+1);
        
        dE_dr_l = dE_dr_l + 2/Material.k0*Rg_a*((-2*Aux.Mau.*Aux.MC1./tmp_MC3.*Invariant.Nshape) ...
            + 1./tmp_MC2.*(Aux.MC1.*Aux.da_dr_l + Aux.Mau.*Invariant.Nshape)).*tmp_pfn.*Aux.Mgweight;
        
        dE_dz_l = dE_dz_l + 2/Material.k0*Rg_a./tmp_MC2.*Aux.MC1.*Aux.da_dz_l.*tmp_pfn.*Aux.Mgweight;
    end
    ns = length(dE_dr_l)/gaussno;
    vs = [1:1:ns];
    dE_dCPr = reshape(sum(reshape(dE_dr_l,gaussno,ns,p+1),1),ns,p+1);
    dE_dCPz = reshape(sum(reshape(dE_dz_l,gaussno,ns,p+1),1),ns,p+1);
    dE_dCP = zeros(2,ini_Val.ncontrol);

    for i = 0:p
        dE_dCP(1,vs(:)+i) = dE_dCP(1,vs(:)+i) + dE_dCPr(:,i+1)';
        dE_dCP(2,vs(:)+i) = dE_dCP(2,vs(:)+i) + dE_dCPz(:,i+1)';
    end

    dE_dCP=2*pi*dE_dCP*Material.k0/2;
end
