function [f,g] = Incremental_Mech_Rayleighian(x)
global Invariant ini_Val Material Constraints ALG phi C_1 dt
%Initialization
C0=Constraints.C0;
p = ini_Val.p;
 
if isfield(Constraints,'f0')
    f0 = Constraints.f0;
else
    f0 = 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,phi);
 
 
 
[cE,dcE_dCP] = Curvature_functional(Aux);

if isfield(Material,'k_dc') && isfield(Constraints,'D0')
    [DE,dDE_dCP] = Deviator_functional(Aux);
end


% 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;
 
 
tE = cE + tU -fE-fzE;
dE_dCP = dcE_dCP + dU_dCP - dfE_dCP - dfzE_dCP; 

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);
 
%[tD,dD_dCP,dD_dCPdot] = Wdissipation_functional(Aux,(x-x_1)/dt);
 
 
f=(tE + tD);
%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);
%g=dE_dCP;
g=reshape(g,length(x),1);
%g = g + dt*dD_dCP + dD_dCPdot;