function [Psi_k] = Psi_update(K,NumdS,dS,rho_k,G_k,theta_k,sigma_2)

Psi_k=zeros(3,K);

for k=1:K
    temp_fenzi=zeros(3,1);
    for n=1:NumdS
        temp_fenzi=temp_fenzi+G_k(:,:,k,n)*theta_k(:,k,n);
    end
    temp_fenzi=temp_fenzi*dS;
    
    temp_fenmu=zeros(3,3);
    for j=1:K
        temp0=zeros(3,1);
        for n=1:NumdS
            temp0=temp0+G_k(:,:,k,n)*theta_k(:,j,n);
        end
        temp0=temp0*dS;      
        temp_fenmu=temp_fenmu+temp0*temp0';
    end
    temp_fenmu=temp_fenmu+sigma_2*eye(3);
    
    Psi_k(:,k)=sqrt(1+rho_k(k))*(temp_fenmu)^(-1)*temp_fenzi;
    
    Psi_k(:,k)=Psi_k(:,k)/norm(Psi_k(:,k));
end


end

