function [eigVals, phi_q] = get1DPSWFs_multipleFreqInterval(t_arr, a, b, freq_intervals, h, Nmax)
% 1D-PSWF parameters (numerical method with SVD+)
% Domain 1: A = [a,b] 
% Domain 2: [W_c - W/2, W_c + W/2] 
% 
% (c, nmax, h) 
% Nmax: max. PSWF order, should be less than nmax 
% Deprecated Model: Operator Qc with integral kernel sin(c*(x-y))/(pi*(x-y)), acting on [-1,1]

    % c = pi*W*(b-a)/2;
    Ni = length(freq_intervals); 
    % format: { [fL1, fH1], [fL1, fH1], ... } 

    x_arr = a+h/2:h:b-h/2;  % equal-spacing sampling 
    Nx = length(x_arr); 

%% Construct the discretized integral kernel K on A = [a,b] specified by the freq-domain concentration intervals. 
    K = zeros(Nx); 
    
    for idx_f = 1:Ni
        interval = freq_intervals{idx_f}; 
        W = interval(2) - interval(1); 
        Wc = mean(interval); 

        for idx_1 = 1:Nx
            for idx_2 = idx_1:Nx
                Delta_x = x_arr(idx_1) - x_arr(idx_2); 
                K(idx_1, idx_2) = K(idx_1, idx_2) + h*W*sinc(W*Delta_x)*exp(1j*2*pi*Wc*Delta_x);  
            end
        end
    end
    
    K = (K+K')-diag(diag(K)); 
    [V, D] = eig(K); 
    
    [sortedEigenvalues, indices] = sort(diag(D), "descend"); 
    eigVals = sortedEigenvalues(1:Nmax); 
    V = V(:, indices); 
    
    Nt = length(t_arr); % Query points 
    phi_q = zeros(Nmax, Nt); 
    
%% Evaluate the PSWF functions on [a,b] and outside [a,b] 
    K_cont = zeros(Nt, Nx); % Analytic continuation of the integral kernels 
    for idx_f = 1:Ni
        interval = freq_intervals{idx_f}; 
        W = interval(2) - interval(1); 
        Wc = mean(interval); 

        for idx_1 = 1:Nt
            for idx_2 = 1:Nx
                Delta_x = t_arr(idx_1) - x_arr(idx_2); 
                K_cont(idx_1, idx_2) = K_cont(idx_1, idx_2) + h*W*sinc(W*Delta_x)*exp(1j*2*pi*Wc*Delta_x); 
            end
        end
    end
    
    % Use the formula: \Pi_2 (\phi_n) = \sqrt{\lambda_n} \psi_n to evaluate
    % the values of \psi_n 
    for idx = 1:Nmax 
        % phi_q(idx, :) = interp1(x_arr, (V(:, idx)).'/sqrt(h), t_arr, "linear", "extrap");  % This is restricted. 
        phi_q(idx, :) = (K_cont * V(:, idx)/sqrt(h)).' / sqrt(eigVals(idx));  
    end
    
    
end


