function p = get_phi(Theta,pmf, SNR_dB)

NN = 50;        % number of antennas
MM = 10;        % number of processing channels
TT = 100;       % number of snapshots


NUser = length(Theta);

if nargin < 3
    SNR_dB = [0,20,20,20,20,20,20,20,20].';
end
%SNR_dB = [0,10,10,10,10,10,10,10,10].';

sigmas = 10.^(SNR_dB/10);
sigmas2_assume = 100;

deltaTheta = 0.1;
Theta_ = -90:deltaTheta:90;


gamma = 0.001;          % step size
nloop = 100;



%% Signal generation


S = (sqrt(sigmas/2)).*(randn(NUser,TT)+1j*randn(NUser,TT));
Noise = (sqrt(1/2))*(randn(NN,TT)+1j*randn(NN,TT));
A = steer(NN,0.5,Theta);

X = A*S + Noise;

Phi = randn(MM,NN)+1j*randn(MM,NN);
PP = Phi*Phi';
Phi0 = PP^-0.5 * Phi;

%% Update of \Phi

for iloop=1:nloop

    sum1 = 0;
    sum2 = 0;
    sum3 = 0;



    for kk=1:length(Theta_)
        theta_kk = Theta_(kk);
        steer_kk = steer(NN,0.5,theta_kk);
        SI_mat = sigmas2_assume * steer_kk * steer_kk' + eye(NN);
        Cyy = Phi * SI_mat * Phi';

        sum1 = sum1 + pmf(kk)/det(Cyy);
        sum2 = sum2 + pmf(kk)/det(Cyy)*inv(Cyy)*Phi*SI_mat;
        sum3 = sum3 + pmf(kk)         *inv(Cyy)*Phi*SI_mat;
    end

    %Phi = Phi0 + gamma * (sum2/sum1 - sum3);    % in paper
    Phi = Phi0 + gamma * (sum2/sum1 - sum3 * 1e5);    % variation


    PP = Phi*Phi';     
    Phi = PP^-0.5 * Phi; 
    Phi0 = Phi;
end

Y = Phi*X;
Ryy = Y*Y';

index=0;
for ii=-90:0.1:90
    index=index+1;
    aa = steer(NN,0.5,ii);
    p(index) = real(aa'*Phi'*Phi*aa / (NN*aa'*Phi'*inv(Ryy)*Phi*aa));
end



end