function rho_k = rho_k_update(K,NumdS,dS,G_k,Psi_k,theta_k)

rho_k=zeros(K,1);

for k=1:K   
    eps=zeros(3,1);
    for n=1:NumdS
        eps=eps+G_k(:,:,k,n)*theta_k(:,k,n);
    end
    eps=real((Psi_k(:,k))'*eps*dS);
    
    rho_k(k)=(eps+eps*sqrt(eps^2+4))/2;
end


end

