%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [var,facs,fact]=updt_den(F,var,mat,ndim)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% density update for functional adaption
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% input:  F        -> deformation gradient
%         var      -> var = rho / rho_0 - 1 internal variable density
%         mat      -> material parameters
% output: var      -> var = rho / rho_0 - 1 internal variable density
%         facs     -> factor for stresses    
%         fact     -> factor for tangent operator    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tol = 1e-8;      

emod = mat(1);   nue  = mat(2);   rho0 = mat(3);   psi0 = mat(4);   
expm = mat(5);   expn = mat(6);   dt   = mat(7);

xmu = emod / 2.0 / (1.0+nue);
xlm = emod * nue / (1.0+nue) / ( 1.0-2.0*nue ); 

J  = det(F);
C  = F'*F;
I1 = trace(C);
psi0_neo = xlm/2 * log(J)^2  +  xmu/2 *(I1 - ndim - 2*log(J));
       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% euler backward - implicit time integration (quad.conv.@large str)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if psi0_neo > tol 
          
rho_k0 = (1 + var) * rho0;
rho_k1 = (1 + var) * rho0; 

iter = 0;   res  = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% local newton-raphson iteration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    while abs(res) > tol
      iter=iter+1;
  	  
      res =((rho_k1/rho0)^(expn-expm)*psi0_neo-psi0)*dt-rho_k1+rho_k0; 
	  dres= (expn-expm)*(rho_k1/rho0)^(expn-expm)*psi0_neo*dt/rho_k1-1;
	  
      drho=- res/dres;  
	  rho_k1 = rho_k1+drho;  
      
	  if(iter>20); disp(['*** NO LOCAL CONVERGENCE ***']); return; end;      
    end
%%% local newton-raphson iteration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
rho = rho_k1;
else
rho = rho0;
end

var = rho / rho0 - 1;
facs= (rho/rho0)^expn;

if (dt~=0)
  facr = 1/dt - (expn-expm) * (rho/rho0)^(expn-expm) / rho *psi0_neo;
  fact=          expn / rho * (rho/rho0)^(    -expm) / facr;
else
  fact = 0.0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%