function [h_approx, logL] = fit_GPRbes(y, simPar)

    sigma2  = simPar.sigma2_n; 
    N_BS    = simPar.N; 

    % thetas          = ((-(2*N_BS-1)/2:(2*N_BS-1)/2)*(1/N_BS));
    Nig         = 10; %?
    thetas      = (linspace(-1+1/Nig, 1-1/Nig, Nig)); 
    [A, ~]          = iget_A(N_BS, thetas);
    [~, idx_max]    = max(abs(A'*y));%?FFT
    theta_max       = asin(thetas(idx_max)); %?
    mu              = 20; %?
    
    lr_theta    = 5e-4;
    lr_mu       = 4e-1; 
    decay       = 0.7; %lr decay

    for iter = 1:simPar.channelEstimator.GPR.N_iter
        [R, dR_dtheta, dR_dmu]  = iget_R(N_BS, theta_max, mu);
        K                       = R + sigma2*eye(N_BS);
        alpha                   = K\y; 
        logL                    = real(-y'*alpha - log(det(K)) - N_BS*log(pi));%???
        dlogL_dtheta            = real(alpha'*dR_dtheta*alpha-trace(K\dR_dtheta));
        
        % opt for theta问一下
        flag = true;
        while flag
            tnew = theta_max + lr_theta*dlogL_dtheta;
            [R, ~, dR_dmu] = iget_R(N_BS, tnew, mu);
            K = R + sigma2*eye(N_BS);
            alpha = K\y; 
            logLnew = real(-y'*alpha - log(det(K)) - N_BS*log(pi));
            if logLnew < logL
                lr_theta = lr_theta*decay; 
            else
                flag = false; 
            end
        end
        theta_max   = tnew; 
        logL        = logLnew; 
        dlogL_dmu   = real(alpha'*dR_dmu*alpha-trace(K\dR_dmu));

        % opt for mu
        flag = true;
        while flag
            mnew = mu + lr_mu*dlogL_dmu;
            [R, ~, ~] = iget_R(N_BS, theta_max, mnew);
            K = R + sigma2*eye(N_BS);
            alpha = K\y; 
            logLnew = real(-y'*alpha - log(det(K)) - N_BS*log(pi));
            if logLnew < logL
                lr_mu = lr_mu*decay; 
            else
                flag = false; 
            end
        end
        mu      = mnew; 
        logL    = logLnew; 
        
        % disp([iter, theta_max, mu, logL]);
    end

    h_approx = R*alpha; 
    
end


function [R, dR_dt, dR_dmu] = iget_R(n, theta, mu)
    % Assume ULA with lambda/2-spacing. 
    R_in = zeros(n);
    dR_dt = zeros(n);
    dR_dmu = zeros(n); 

    for ii = 1:n
        for jj = 1:n

            z = sqrt((+1i*pi*(ii-jj)+mu*sin(theta)).^2 + (mu*cos(theta)).^2); 
            R_in(ii, jj) = besseli(0, z);
            dR_dt(ii, jj) = besseli(1, z)*(1/(2*z))*...
                (2*(+1i*pi*(ii-jj)+mu*sin(theta))*mu*cos(theta) + 2*(mu*cos(theta))*(-mu*sin(theta)));
            
            dR_dmu(ii, jj) = besseli(1, z)*(1/(2*z))*...
                (2*(+1i*pi*(ii-jj)+mu*sin(theta))*sin(theta) + 2*(mu*cos(theta)^2));
        end
    end
    I0 = besseli(0, mu); 
    I1 = besseli(1, mu); 

    R = R_in / I0;
    dR_dt = dR_dt / I0; 
    dR_dmu = dR_dmu/I0 + R_in*(-1/(I0^2))*I1; 
end

function [A, dA] = iget_A(N_BS, thetas)
% This function returns the far-field dictionary matrix defined by thetas. 
    S = length(thetas);
    A = zeros(N_BS, S);
    dA = zeros(N_BS, S);
    
    seq = (0:N_BS-1).';
    for idx = 1:S
        A(:,idx) = exp(1j*pi*thetas(idx)*seq)/sqrt(N_BS);
        dA(:,idx) = A(:,idx).*(1j*pi*seq);
    end
end


