function [f,g] = Incremental_Mech_Rayleighian_S(x)
global Invariant ini_Val Material Constraints ALG C_1 dt C_psi
%Initialization
C0=Constraints.C0;
p = ini_Val.p;

if isfield(Constraints,'f0')
    f0 = Constraints.f0;
else
    f0 = 0;
end

% To check for fixed pressure or fixed volume ensemble between membrane and
% support
if isfield(Constraints,'iVol_s')
    iVol_s = Constraints.iVol_s;
    if(iVol_s ==0)
        Pr = Constraints.Pr_S;
    end
else
    iVol_s = 0;
    Pr = 0;
end


if isfield(Constraints,'f0z') %Point load applied at the tip of bud
    f0z = Constraints.f0z;
else
    f0z = 0;
end


[C,C1,C2]=interpolate(Invariant,x);
Aux = auxiliary_calculations(Invariant,p,C0,C,C1,C2);


[cE,dcE_dCP] = Curvature_functional(Aux);

%[cE,dcE_dCP] = Curvature_functional_Reg(Aux,C);




% potential for interaction membrane - support
if isfield(Material,'a') && isfield(Material,'b') &&  isfield(Material,'h')
    [tU,dU_dCP] = incremental_potential(Aux);
else
    tU = 0 ;
    dU_dCP = 0;
end


fE = pi* x(end-1)^2*f0;
tmp = zeros(size(dcE_dCP));
tmp(1,end) = 2*pi*x(end-1)*f0;
dfE_dCP = tmp;

fzE = x(2)*f0z;
tmpz = zeros(size(dcE_dCP));
tmpz(2,1) = x(2)*f0z;
dfzE_dCP = tmpz;

%Elastic_outer = Material.K_outer/2*(x(end-1)-Material.R_ini)^2;
%Elastic_outer = 0;
%if x(2) < 6
%    Elastic_outer = Material.K_outer/2*(x(2)-6)^2;
%end

%tE = cE + tU -fE-fzE + Elastic_outer;
tE = cE + tU -fE-fzE;
dE_dCP = dcE_dCP + dU_dCP - dfE_dCP - dfzE_dCP;

if (iVol_s==0)
    [tV,dV_dCP] = Volume_functional_support(Aux);
    tE = tE - Pr.*tV;
    dE_dCP = dE_dCP - Pr.*dV_dCP;
end


if isfield(Material,'k_dc') && isfield(Constraints,'D0')
    tE=tE+DE;
    dE_dCP = dE_dCP + dDE_dCP;
end

[tD,dD_dCP] = Incremental_mem_diss(Aux,C,C_1,C1);
[D_s,dD_dCP_s] = Incremental_drag_diss(Aux,C,C_1);


%[tD,dD_dCP,dD_dCPdot] = Wdissipation_functional(Aux,(x-x_1)/dt);


f=(tE + tD + D_s);


%f=tE + dt * tD;
%dD_dCP=reshape(dD_dCP,length(x),1);
%D_W = matrix_W(Aux,C,C1,C2);
%dD_dCPdot=D_W*(x-x_1)/dt;

g=(dE_dCP +  dD_dCP + dD_dCP_s);

if isfield(Material,'xi_Cor_perp')
    if Material.xi_Cor_perp~=0
        [tD_c,dD_dCP_c] = Incremental_cortex_diss(Aux,C,C_1,C1);
        f = f + tD_c;
        g = g + dD_dCP_c;
    end
end
%g(end-1) = g(end-1)  + Material.K_outer * (x(end-1)-Material.R_ini);
%if x(2) < 6
%    g(2) = g(2) +  Material.K_outer*(x(2)-6);
%end

%g=dE_dCP;
g=reshape(g,length(x),1);
%g = g + dt*dD_dCP + dD_dCPdot;


