function [eigVals, phi_q] = get1DPSWFs_numerical(t_arr, a, b, Wc, W, h, nmax, 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;
    x_arr = a+h/2:h:b-h/2; 
    Nx = length(x_arr); 

%% Construct the integral kernel on Domain 1: A = [a,b] 
    K = zeros(Nx); 

    for idx_1 = 1:Nx
        for idx_2 = idx_1:Nx
            Delta_x = x_arr(idx_2) - x_arr(idx_1); 
            K(idx_1, idx_2) = h*W*sinc(W*Delta_x)*exp(1j*2*pi*Wc*Delta_x);  
        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); 

    phi_q = zeros(Nmax, length(t_arr)); 
    for idx = 1:Nmax
        phi_q(idx, :) = interp1(x_arr, (V(:, idx)).'/sqrt(h), t_arr, "linear", "extrap"); 
    end
    
    
end


