function [tS,dS_dCP] = Area_functional(Aux)
global Invariant ini_Val


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

tS = 2*pi*sum((Aux.ds).*Invariant.gweight);

if nargout == 2
    %Derivatives of S (local)
    dS_dr_l = (Aux.MC1.*Aux.da_dr_l+Aux.Mau.*Invariant.Nshape).*Aux.Mgweight;
    dS_dz_l = Aux.MC1.*Aux.da_dz_l.*Aux.Mgweight;
    
    ns = length(dS_dr_l)/gaussno;
    vs = [1:1:ns];
    dS_dCPr = reshape(sum(reshape(dS_dr_l,gaussno,ns,p+1),1),ns,p+1);
    dS_dCPz = reshape(sum(reshape(dS_dz_l,gaussno,ns,p+1),1),ns,p+1);
    dS_dCP = zeros(2,ini_Val.ncontrol);
    
    for i = 0:p
        dS_dCP(1,vs(:)+i) = dS_dCP(1,vs(:)+i) + dS_dCPr(:,i+1)';
        dS_dCP(2,vs(:)+i) = dS_dCP(2,vs(:)+i) + dS_dCPz(:,i+1)';
    end
    dS_dCP=2*pi*dS_dCP;
end