function [Hk, Rut, Rur, RT, RR] = get2DColoredAngularFadingChannel(MIMO_params, RegionsTx, RegionsRx)
% Regions: subset of D. 
persistent innerRut;
persistent innerRur; 

    nkGrids = MIMO_params.nkGrids; 
    [Ngp, ~] = size(nkGrids);   % Grid points. 
    
%     figure; 
%     scatter(nkGrids(:,1), nkGrids(:,2), 'blue');  hold on;  axis equal; 

    RT = zeros(Ngp); 
    LRT = length(RegionsTx);
    RR = zeros(Ngp); 
    LRR = length(RegionsRx); 
    
    doMemorize = false;
    doCalc = true; 
    if (~isfield(MIMO_params, 'memorizeCorrelation') || ~MIMO_params.memorizeCorrelation) || (MIMO_params.isULA && MIMO_params.FreqShift)
        
    else
        if exist('innerRut', 'var') && ~isempty(innerRut) && all(size(innerRut) == size(RT))
            RT = innerRut; 
            RR = innerRur; 
            doCalc = false; 
        else
            doMemorize = true; 
        end
    end

    if doCalc
        if ~MIMO_params.isSIMO
        
            for idx_lr = 1:LRT
                reg = RegionsTx{idx_lr}; 
                sel = sum((nkGrids - reg(1:2)).^2, 2) < reg(3)^2; 
                if ~any(sel)
                    sel(randi(Ngp)) = true; 
                end
                selectedGrids = nkGrids(sel, :); 
        
                L = sum(sel); 
                Corr = zeros(L, L);
                dfFactor = MIMO_params.FreqShiftFactor*randn();    % dispersive frequency shift factor
        
                if MIMO_params.debug
                    fprintf('%.2f ', dfFactor); 
                end
        
                for idx = 1:L
                    if MIMO_params.isULA && MIMO_params.FreqShift
                        temp = selectedGrids - selectedGrids(idx, :); 
                        rhos = temp(:,1); 
                        Corr(idx, :) = getSincCorrWithRandomFrequencyShift(rhos, reg, dfFactor); 
                    else
                        rhos = sqrt(sum((selectedGrids - selectedGrids(idx, :)).^2, 2));
                        Corr(idx, :) = getSincCorr(rhos, reg); 
                    end
                    
                end
        
                RT(sel, sel) = Corr; 
        
            end
        end
        
    
        for idx_lr = 1:LRR
            reg = RegionsRx{idx_lr}; 
            sel = sum((nkGrids - reg(1:2)).^2, 2) < reg(3)^2; 
            if ~any(sel)
                sel(randi(Ngp)) = true; 
            end
            selectedGrids = nkGrids(sel, :); 
    
            L       = sum(sel); 
            Corr    = zeros(L, L);
            dfFactor = MIMO_params.FreqShiftFactor*randn();    % dispersive frequency shift factor
            if MIMO_params.debug
                fprintf('%.2f ', dfFactor); 
            end
            
            for idx = 1:L
                if MIMO_params.isULA && MIMO_params.FreqShift
                    temp = selectedGrids - selectedGrids(idx, :); 
                    rhos = temp(:,1); 
                    Corr(idx, :) = getSincCorrWithRandomFrequencyShift(rhos, reg, dfFactor); 
                else
                    rhos = sqrt(sum((selectedGrids - selectedGrids(idx, :)).^2, 2));
                    Corr(idx, :) = getSincCorr(rhos, reg); 
                end
            end
    
            RR(sel, sel) = Corr; 
        end
    end
    
    % Sample a k-domain channel kH with correlation specified by (RR, RT).
    if ~MIMO_params.isSIMO
        if ~doCalc
            Rut = innerRut; 
            Rur = innerRur; 
        else
            Rut = mysqrtm(RT); 
            Rur = mysqrtm(RR); 
            if doMemorize
                innerRut = Rut; 
                innerRur = Rur; 
            end
        end
        W = (randn(Ngp) + 1j*randn(Ngp))/sqrt(2); 
        Hk = (Rur)*W*(Rut.');  

    else
        if ~doCalc
            Rur = innerRur; 
        else
            Rur = mysqrtm(RR); 
            if doMemorize
                innerRur = Rur; 
            end
        end
        w = (randn(Ngp, 1) + 1j*randn(Ngp, 1))/sqrt(2); 
        Hk = Rur * w; 
        Rut = []; 
    end  

end


function corrs = getBesselCorr(rhos, reg)
    L = length(rhos);
    corrs = zeros(size(rhos));
    eps = 1e-5; 
    for idx = 1:L
        if rhos(idx)/reg(4) < eps
            corrs(idx) = 1;
        else
            corrs(idx) = 2*besselj(1, rhos(idx)/reg(4))./(rhos(idx)/reg(4)); 
        end
    end
end

function corrs = getSincCorr(rhos, reg)
    corrs = sinc(rhos/reg(4))/reg(4);
    % corrs = sinc(rhos/reg(4)); 
end

function corrs = getSincCorrWithRandomFrequencyShift(rhos, reg, factor)
    corrs = sinc(rhos/reg(4))/reg(4).*(exp(1i*factor*rhos));
end


function Rhalf = mysqrtm(R)
    [V, D] = eig(R);
    d = diag(D); 
    
    [~, order] = sort(d, 'descend'); 
    d = d(order); 
    V = V(:, order); 

    totE = sum(abs(d).^2); 
    eps = 1e-10; 
    sel = and((real(d)>0), ((abs(d).^2/totE) >= eps)); 
    
    Vr = V(:,sel);
    Rhalf = Vr*diag(sqrt(d(sel)))*(Vr');
end
