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


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

if isfield(Material,'ka')
    [tS,dS_dCP] = Area_functional(Aux);
    %[S_L,dS_L_dCP] = localArea_functional(Aux);
    tE = 2*pi*sum(Material.k0/2*Aux.HmC02.*Aux.ds.*Invariant.gweight)...
        +Material.ka.*tS;
    
elseif isfield(Constraints,'anisotropy') && Constraints.anisotropy ==1
    %Aux.HmC0_bis=(Aux.Hu-C0).*phi; Aux.HmC02_bis=(Aux.Hu-C0).^2.*phi;

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

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

if nargout == 2
    %Derivatives of E (local)
    if isfield(Constraints,'anisotropy') && Constraints.anisotropy == 1
        temp  = Material.k_b/2.*(Aux.MHu2);
        temp1 = Material.k0/2.*Aux.MHmC02;
        temp2 = Material.k_b/2.*2*(Aux.MHu).*Aux.Mau.*Aux.MC1;
        temp3 = Material.k0/2.*2*(Aux.MHmC0).*Aux.Mau.*Aux.MC1;
        
        dE_dr_l = (((temp2+temp3).*Aux.dH_dr_l)+((temp+temp1).*(Aux.MC1.*Aux.da_dr_l ...
            +Aux.Mau.*Invariant.Nshape))).*Aux.Mgweight;
        dE_dz_l = (((temp2+temp3).*Aux.dH_dz_l)+((temp+temp1).*(Aux.MC1.*Aux.da_dz_l))) ...
            .*Aux.Mgweight;


    else    
        temp = Material.k0/2.*Aux.MHmC02;
        temp1 = Material.k0/2.*2*(Aux.MHmC0).*Aux.Mau.*Aux.MC1;
        dE_dr_l = (((temp1.*Aux.dH_dr_l)+(temp.*(Aux.MC1.*Aux.da_dr_l ...
            +Aux.Mau.*Invariant.Nshape))).*Aux.Mgweight);
        dE_dz_l = (((temp1.*Aux.dH_dz_l)+(temp.*(Aux.MC1.*Aux.da_dz_l))) ...
            .*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;
    
    if isfield(Material,'ka')
        dE_dCP=dE_dCP+dS_dCP*Material.ka;
    end
    
    
    
    
end
