function [H, Correlation] = getChannel(MIMO_params, method)
    
    if method == "ColoredAngularCorrelation"
        
        if MIMO_params.isSIMO
            [H2D, ~, halfRr, ~, Rr] = get2DColoredAngularFadingChannel(MIMO_params, [], MIMO_params.RxArray.RegionsRx); 
            
            if nargout == 2
                Correlation.halfRr = halfRr; 
            end

        else
            [H2D, halfRt, halfRr, Rt, Rr] = get2DColoredAngularFadingChannel(MIMO_params, MIMO_params.TxArray.RegionsTx, MIMO_params.RxArray.RegionsRx); 

            if nargout == 2
                Correlation.halfRt = halfRt; 
                Correlation.halfRr = halfRr; 
            end
        end

        % Plot this channel 
%         figure(1); 
%         plotChannel(H2D, MIMO_params, "kDomain", gca);
%         figure(2); 
%         plotChannel(H2D.', MIMO_params, "kDomain", gca); 

        if MIMO_params.debug
            fprintf('\n'); 
            figure;
            imagesc(abs(H2D)); axis equal; colorbar; 
            title('k-domain strength')

            figure; 
            imagesc(angle(H2D)); axis equal; 
            title('k-domain phase'); 
            colorbar; 
        end
        
        Dr = MIMO_params.RxArray.dict; 
        

        if MIMO_params.isSIMO
            H = (Dr)*H2D*(2/MIMO_params.K); 

            if isfield(MIMO_params, 'normalizeChannel') && MIMO_params.normalizeChannel
                CovH = (2/MIMO_params.K)^2*Dr*Rr*Dr'; 
                totChEnergy = real(trace(CovH)); 

                H = H/sqrt(totChEnergy); 
                Correlation.CovH = CovH/totChEnergy;  
            end
            
        else
            Dt = MIMO_params.TxArray.dict; 
            H = (Dr)*H2D*(Dt')*(2/MIMO_params.K)^2; 
            
            if isfield(MIMO_params, 'normalizeChannel') && MIMO_params.normalizeChannel
                A = conj(Dt)*Rt*(Dt.'); 
                B = Dr*Rr*Dr'; 

                totChEnergy = real(trace(A))*real(trace(B))*(2/MIMO_params.K)^4; 
                H = H/sqrt(totChEnergy); 

                if isfield(MIMO_params, 'forceGetCovariance') && MIMO_params.forceGetCovariance
                    % return covariance matrix of vec(H). 
                    Correlation.CovH = (kron(A, B)*(2/MIMO_params.K)^4)/totChEnergy; 
                end
            end
            
        end

        if MIMO_params.debug
            figure;
            imagesc(abs(H)); axis equal; colorbar; 
            title('Channel Matrix'); 
        end

    elseif method == "CDL"
        % Initialize the CDL-MIMO model. 
        cdl_model = nrCDLChannel('ChannelFiltering', false);
        cdl_model.DelayProfile              = 'CDL-A'; 
        cdl_model.MaximumDopplerShift       = 300;
        cdl_model.CarrierFrequency          = MIMO_params.fc;
        cdl_model.SampleRate                = 10e2;     
        cdl_model.SampleDensity             = Inf;
        cdl_model.NumTimeSamples            = 1;
        cdl_model.TransmitAntennaArray.Size = [MIMO_params.TxArray.N, 1, 1, 1, 1];   % [M N P Mg Ng]
        cdl_model.ReceiveAntennaArray.Size  = [MIMO_params.RxArray.N, 1, 1, 1, 1]; 
        cdl_model.Seed                      = randi(10000);
        
        [pathGains, ~] = cdl_model();
        NRx = MIMO_params.RxArray.N; NTx = MIMO_params.TxArray.N;
        H = reshape(sum(pathGains(1, :, :, :), 2), [NTx, NRx]).';
        
        if isfield(MIMO_params, 'CDL_corrMat')
            covH = MIMO_params.CDL_corrMat; 
            H = H/sqrt(trace(covH));  
            covH = covH/trace(covH); 

            Correlation.CovH = covH; 
            
        else
            error('CDL Correlation characteristics unknown.'); 
        end

    else
        error('Method Error');
    end
    
    % H = H/norm(H, 'fro'); 

end
