function [phi,x]=hat2phi(phihat,rho)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% function phi=hat2phi(phihat,rho,x)
%
% rho must be given as a vector (assumed diagonal)
%
% diag(rho) = phi^{1/2} phihat phi^{1/2}

[n,m]=size(rho);
if m~=1, disp('rho must be n by 1'), return, end

rhohalf=diag(sqrt(rho));
rhohalfinv=inv(rhohalf);

x=(rhohalf*sqrtm(rhohalfinv*inv(phihat)*rhohalfinv)*rhohalf);
phi=x*x;

% rhohalf=diag(sqrt(rho./diag(phihat)));
% rhohalfinv=inv(rhohalf);
% 
% x=(rhohalf*sqrtm(rhohalfinv*inv(phihat)*rhohalfinv)*rhohalf);
% phi=x*x;
