function MIMO_params = setup(method, fc, L, specifiedRegions, spacing, debugging)
    c = physconst('lightspeed'); 
    lambda = c/fc; 

    % Check for missing inputs. 
    if ~exist('debug', 'var')
        debugging = false;
    end

    MIMO_params.lambda = lambda; 
    MIMO_params.fc = fc; 
    if ~exist('L', 'var')
        L = 0.8;
    end

    if ~exist('spacing', 'var')
        spacing = lambda/2; 
    end
    
    if ~exist('specifiedRegions', 'var')
        Gamma = 0.05; 
    
        RegionsTx = cell(1,1);
        RegionsTx{1} = [-0.4, 0, 0.15, Gamma];   % format: (nkx, nky, nkradius, gammaCorrelationFactor). 
        RegionsTx{2} = [0.3, 0, 0.15, Gamma];
        RegionsTx{3} = [0.7, 0, 0.15, Gamma]; 
        
        RegionsRx = cell(1,1); 
        RegionsRx{1} = [-0.85, 0, 0.15, Gamma]; 
        RegionsRx{2} = [0.1, 0, 0.15, Gamma];
        RegionsRx{3} = [0.52, 0, 0.15, Gamma];

    else
        
        Gamma = specifiedRegions.Gamma; 
        RegionsTx = specifiedRegions.RegionsTx; 
        RegionsRx = specifiedRegions.RegionsRx; 
    end

    
    if method == "SymmetricUPA"
        UPA.Lx = L;
        UPA.Ly = L;
        UPA.dx = spacing; 
        UPA.dy = spacing; 
        UPA.x = -UPA.Lx/2:UPA.dx:UPA.Lx/2; 
        UPA.y = -UPA.Ly/2:UPA.dy:UPA.Ly/2;
        UPA.Nx = length(UPA.x);
        UPA.Ny = length(UPA.y); 
        UPA.N = UPA.Nx*UPA.Ny; 
        
        dsep = 3; 
        theta_incl = 0; 
        
        UPA.coord = zeros(UPA.Nx, UPA.Ny, 3); 
        for idx = 1:UPA.Nx
            for idy = 1:UPA.Ny
                UPA.coord(idx, idy, :) = [UPA.x(idx), UPA.y(idy), 0]; 
            end
        end
        
        upa(1) = UPA; 
        
        for idx = 1:UPA.Nx
            for idy = 1:UPA.Ny
                UPA.coord(idx, idy, :) = [UPA.x(idx), UPA.y(idy)*cos(theta_incl), dsep+UPA.y(idy)*sin(theta_incl)]; 
            end
        end
        
        upa(2) = UPA;
        
        MIMO_params.TxArray = upa(1); 
        MIMO_params.RxArray = upa(2); 
        
        % Construct nk-domain grids with EM-compliance 
        [Gx, Gy] = meshgrid((-1+1/UPA.Nx):(2/UPA.Nx):(1-1/UPA.Nx), (-1+1/UPA.Ny):(2/UPA.Ny):(1-1/UPA.Ny)); 
        sel = (Gx.^2+Gy.^2) <= 1;
        sel = sel(:); 
        Gx = Gx(:);
        Gy = Gy(:); 
        
        Grids = zeros(sum(sel), 2);
        cnt = 1; 
        for idx = 1:(UPA.Nx*UPA.Ny)
            if sel(idx)
                Grids(cnt, :) = [Gx(idx), Gy(idx)]; 
                cnt = cnt + 1;
            end
        end
        
        MIMO_params.nkGrids = Grids; % normalized k-domain grids. 
        [Ng, ~] = size(Grids); 
        dictTx = zeros(MIMO_params.TxArray.Nx*MIMO_params.TxArray.Ny, Ng); 
        dictRx = zeros(MIMO_params.RxArray.Nx*MIMO_params.RxArray.Ny, Ng); 
        
        for idx = 1:Ng
            t = MIMO_params.nkGrids(idx, :);
            nkx = t(1); nky = t(2); nkz = sqrt(1-nkx^2-nky^2);  % normalized k out. 
            nk = [nkx, nky, nkz].'; 
            
            ct = reshape(MIMO_params.TxArray.coord, [MIMO_params.TxArray.Nx*MIMO_params.TxArray.Ny, 3]); 
            cr = reshape(MIMO_params.RxArray.coord, [MIMO_params.RxArray.Nx*MIMO_params.RxArray.Ny, 3]); 
            pst = ct * (2*pi*nk/lambda);
            psr = cr * (2*pi*nk/lambda); 
        
            dictTx(:, idx) = exp(-1i*pst(:));
            dictRx(:, idx) = exp(-1i*psr(:)); 
        end
        
        MIMO_params.TxArray.dict = dictTx; 
        MIMO_params.RxArray.dict = dictRx; 
        MIMO_params.TxArray.RegionsTx = RegionsTx; 
        MIMO_params.RxArray.RegionsRx = RegionsRx; 

    
    elseif method == "SymmetricULAwithFixedK"
        MIMO_params.isSIMO = false; 
        MIMO_params.isULA = true; 
        ULA.Lx = L;
        ULA.dx = spacing; 
        ULA.x = -ULA.Lx/2:ULA.dx:ULA.Lx/2; 
        ULA.Nx = length(ULA.x);
        ULA.N = ULA.Nx; 
        
        dsep = 0; 
        theta_incl = 0; 
        
        ULA.coord = zeros(ULA.Nx, 3); 
        for idx = 1:ULA.Nx
            ULA.coord(idx, :) = [ULA.x(idx), 0, 0]; 
        end
        
        ula(1) = ULA; 
        
        for idx = 1:ULA.Nx
            ULA.coord(idx, :) = [ULA.x(idx), 0, dsep+ULA.x(idx)*sin(theta_incl)]; 
        end
        
        ula(2) = ULA;
        
        MIMO_params.TxArray = ula(1); 
        MIMO_params.RxArray = ula(2); 
        
        K = 2^10;           % Aligned with the TIT'14 paper 
        Gx = (-1+1/K:(2/K):(1-1/K)); 
        
        Grids = zeros(length(Gx), 2);
        cnt = 1; 
        for idx = 1:K
            Grids(cnt, :) = [Gx(idx), 0]; 
            cnt = cnt + 1;
        end
        
        MIMO_params.nkGrids = Grids; % normalized k-domain grids. 
        [Ng, ~] = size(Grids); 
        dictTx = zeros(MIMO_params.TxArray.Nx, Ng); 
        dictRx = zeros(MIMO_params.RxArray.Nx, Ng); 
        
        for idx = 1:Ng
            t = MIMO_params.nkGrids(idx, :);
            nkx = t(1); nky = t(2); nkz = sqrt(1-nkx^2-nky^2);  % normalized k out. 
            nk = [nkx, nky, nkz].'; 
            
            ct = reshape(MIMO_params.TxArray.coord, [MIMO_params.TxArray.Nx, 3]); 
            cr = reshape(MIMO_params.RxArray.coord, [MIMO_params.RxArray.Nx, 3]); 
            pst = ct * (2*pi*nk/lambda);
            psr = cr * (2*pi*nk/lambda); 
        
            dictTx(:, idx) = exp(-1i*pst(:));
            dictRx(:, idx) = exp(-1i*psr(:)); 
        end
        
        MIMO_params.TxArray.dict = dictTx; 
        MIMO_params.RxArray.dict = dictRx; 
        MIMO_params.TxArray.RegionsTx = RegionsTx; 
        MIMO_params.RxArray.RegionsRx = RegionsRx; 

        if debugging
            fprintf('Configuration: MIMO ULA with %d Antennas\n', MIMO_params.RxArray.N); 
        end

    elseif method == "SIMOULAwithFixedK"
        MIMO_params.isULA = true; 
        MIMO_params.isSIMO = true; 

        ULA.Lx = L;
        ULA.dx = spacing; 
        ULA.x = -ULA.Lx/2:ULA.dx:ULA.Lx/2; 
        ULA.Nx = length(ULA.x);
        ULA.N = ULA.Nx; 
        
        ULA.coord = zeros(ULA.Nx, 3); 
        for idx = 1:ULA.Nx
            ULA.coord(idx, :) = [ULA.x(idx), 0, 0]; 
        end
       
        MIMO_params.RxArray = ULA; 
        
        K = 2^10; 
        Gx = (-1+1/K:(2/K):(1-1/K)); 
        
        Grids = zeros(length(Gx), 2);
        cnt = 1; 
        for idx = 1:K
            Grids(cnt, :) = [Gx(idx), 0]; 
            cnt = cnt + 1;
        end
        
        MIMO_params.nkGrids = Grids; % normalized k-domain grids. 
        [Ng, ~] = size(Grids); 
        dictRx = zeros(MIMO_params.RxArray.Nx, Ng); 
        
        for idx = 1:Ng
            t = MIMO_params.nkGrids(idx, :);
            nkx = t(1); nky = t(2); nkz = sqrt(1-nkx^2-nky^2);  % normalized k out. 
            nk = [nkx, nky, nkz].'; 
            
            cr = reshape(MIMO_params.RxArray.coord, [MIMO_params.RxArray.Nx, 3]); 
            psr = cr * (2*pi*nk/lambda); 
        
            dictRx(:, idx) = exp(-1i*psr(:)); 
        end
        
        MIMO_params.RxArray.dict = dictRx; 
        MIMO_params.RxArray.RegionsRx = RegionsRx; 
        
        if debugging
            fprintf('Configuration: SIMO BS-ULA with %d Antennas\n', MIMO_params.RxArray.N); 
        end

    else
        error('Method %s not defined.', method); 
    end

    MIMO_params.K = K;                  % Wavenumber domain discretization number
    MIMO_params.Gamma = Gamma;  


    MIMO_params.FreqShiftFactor = 0; 
    MIMO_params.FreqShift       = false; 
    MIMO_params.debug           = false; 
    MIMO_params.P               = 10;
    MIMO_params.P_channelEst    = 10; 
    MIMO_params.Npilot          = 300; 

    % Methods
    MIMO_params.Estimators.BT_PSWF.eta = 0.5; 

    % Other utils
    if exist('results/CDL_corrMat.mat', 'file')
        load('results/CDL_corrMat.mat', 'corrMat'); 
        MIMO_params.CDL_corrMat = corrMat; 
    end

end
