function Aux = auxiliary_calculations(Invariant,p,C0,C,C1,C2,phi,C3)
global Material Constraints
if nargin==6
    phi = 1;
end


gweight=Invariant.gweight;
Nshape=Invariant.Nshape;
dNsh=Invariant.dNsh;
ddNsh=Invariant.ddNsh;

Aux.au = sqrt(C1(:,1).^2 + C1(:,2).^2);
Aux.bu = -(C2(:,1).*C1(:,2))+(C2(:,2).*C1(:,1));
Aux.dau = ((C1(:,1).*C2(:,1))+(C1(:,2).*C2(:,2)))./Aux.au;
Aux.Hu = ((Aux.bu./(Aux.au.^2))+(C1(:,2)./C(:,1)))./(Aux.au);
Aux.ds = Aux.au.*C(:,1);
Aux.dv = (C(:,1).^2).*C1(:,2)/2;
Aux.dv_supp = C(:,1).*C(:,2).*C1(:,1);
Aux.HmC0=(Aux.Hu-C0*phi); Aux.HmC02=Aux.HmC0.^2;
Aux.HmC0_bis=2.*(Aux.Hu-C0).*phi; Aux.HmC02_bis=(Aux.Hu-C0).^2.*phi;
Aux.K = Aux.bu.*C1(:,2)./(Aux.au.^4.*C(:,1));

% potential support bilayer
   %y_min=0.976491;
   y_min=0;
if isfield(Material,'a') && isfield(Material,'b') &&  isfield(Material,'h')

    Aux.dz2 = (1./(C(:,2)+y_min)).^2;
    Aux.dz3 = (1./(C(:,2)+y_min)).^3;
    Aux.dz4 = exp(-Material.a.*(C(:,2)+y_min));
end

if isfield(Material,'k_dc') && isfield(Constraints,'D0')
    Aux.Du = ((C1(:,2)./C(:,1))-(Aux.bu./(Aux.au.^2)))./(Aux.au);
    Aux.DmD0 = (Aux.Du - Constraints.D0*phi); Aux.DmD02 = Aux.DmD0.^2;
    Aux.DmD0_bis=2.*(Aux.Du-Constraints.D0).*phi; Aux.DmD02_bis=(Aux.Du-Constraints.D0).^2.*phi;
    
end


if nargin==8
    tmp1 =  C1(:,1).*C2(:,1) + C1(:,2).*C2(:,2);
    tmp2 = -C3(:,1).*C1(:,2) + C1(:,1).*C3(:,2);
    tmp3 =  C2(:,2).* C(:,1) - C1(:,2).*C1(:,1);
    Aux.Hpu = (tmp1.*(-3*Aux.bu./Aux.au.^2 - C1(:,2)./C(:,1)) ...
        + tmp2)./(Aux.au.^3) + tmp3 ./C(:,1).^2./ Aux.au;
end

Aux.MC1= C(:,1)*ones(1,p+1);
Aux.MC2= C(:,2)*ones(1,p+1);
Aux.MC12= C1(:,2)*ones(1,p+1);
Aux.MC22= C2(:,2)*ones(1,p+1);
Aux.MC21= C2(:,1)*ones(1,p+1);
Aux.MC11= C1(:,1)*ones(1,p+1);
Aux.Mau= Aux.au*ones(1,p+1);
Aux.Mdau= Aux.dau*ones(1,p+1);
Aux.Mds = Aux.au.*C(:,1)*ones(1,p+1);
Aux.Mbu=Aux.bu*ones(1,p+1);
Aux.Mau2=Aux.Mau.^2;Aux.Mau3=Aux.Mau.*Aux.Mau2;Aux.Mau4=Aux.Mau2.^2;Aux.Mdau2=Aux.Mdau.^2;
Aux.MHu=Aux.Hu*ones(1,p+1); Aux.MHu2=Aux.Hu.^2*ones(1,p+1);
Aux.MHmC0=Aux.HmC0*ones(1,p+1); Aux.MHmC02=Aux.HmC02*ones(1,p+1); 
Aux.MHmC0_bis=Aux.HmC0_bis*ones(1,p+1); Aux.MHmC02_bis=Aux.HmC02_bis*ones(1,p+1);

Aux.Mgweight=gweight*ones(1,p+1);
Aux.y_min=y_min*ones(length(Aux.MC2),p+1);
%Derivatives of b (local)
Aux.db_dr_l = -(Aux.MC12.*ddNsh)+(Aux.MC22.*dNsh);
Aux.db_dz_l = -(Aux.MC21.*dNsh)+(Aux.MC11.*ddNsh);
%Derivatives of a (local)
Aux.da_dr_l = (Aux.MC11./Aux.Mau).*dNsh;
Aux.da_dz_l = (Aux.MC12.*dNsh)./Aux.Mau;
%Derivatives of a'^2 (local)
Aux.ddau_dr_l = 2*((dNsh.*Aux.MC21+Aux.MC11.*ddNsh)-Aux.Mdau.*Aux.da_dr_l).*Aux.Mdau./Aux.Mau;
Aux.ddau_dz_l = 2*((Aux.MC22.*dNsh+Aux.MC12.*ddNsh)-Aux.Mdau.*Aux.da_dz_l).*Aux.Mdau./Aux.Mau;
%Derivatives of H (local)
Aux.dH_dr_l = (Aux.db_dr_l./Aux.Mau3)-(3*Aux.Mbu.*Aux.da_dr_l./Aux.Mau4)- ...
        Aux.MC12.*Nshape./(Aux.Mau.*Aux.MC1.^2)-Aux.MC12.*Aux.da_dr_l./(Aux.Mau2.*Aux.MC1);
Aux.dH_dz_l = (Aux.db_dz_l./Aux.Mau3)-(3*Aux.Mbu.*Aux.da_dz_l./(Aux.Mau4))+ ...
        dNsh./(Aux.Mau.*Aux.MC1)-Aux.MC12.*Aux.da_dz_l./(Aux.Mau2.*Aux.MC1);
%Derivatives of D (local)
if isfield(Material,'k_dc') && isfield(Constraints,'D0')
    Aux.MDu=Aux.Du*ones(1,p+1); 
    Aux.MDmD0=Aux.DmD0*ones(1,p+1); Aux.MDmD02=Aux.DmD02*ones(1,p+1);
    Aux.MDmD0_bis=Aux.DmD0_bis*ones(1,p+1); Aux.MDmD02_bis=Aux.DmD02_bis*ones(1,p+1);
    Aux.dD_dr_l = -(Aux.db_dr_l./Aux.Mau3)+(3*Aux.Mbu.*Aux.da_dr_l./Aux.Mau4)- ...
        Aux.MC12.*Nshape./(Aux.Mau.*Aux.MC1.^2)-Aux.MC12.*Aux.da_dr_l./(Aux.Mau2.*Aux.MC1);
    Aux.dD_dz_l = -(Aux.db_dz_l./Aux.Mau3)+(3*Aux.Mbu.*Aux.da_dz_l./(Aux.Mau4))+ ...
        dNsh./(Aux.Mau.*Aux.MC1)-Aux.MC12.*Aux.da_dz_l./(Aux.Mau2.*Aux.MC1);
end
