function [eigVals, phi_q] = getShifted1DPSWFs(t_arr, a, b, Wc, W, h, nmax, Nmax)
% Shifted PSWF parameters
% 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;

%% Solve standard 1D-PSWF psi(x) with spectral order nmax 
    Pbar = getNormalizedLegendrePolynomials(nmax); 
    x = -1+h/2:h:1-h/2;

    % Construct A matrix. 
    A = zeros(nmax+1); 
    for ii = 1:nmax+1
        k = ii-1; 
        A(ii, ii) = k*(k+1)+(2*k*(k+1)-1)/((2*k+3)*(2*k-1))*c^2;
    end
    
    for ii = 1:nmax-1
        k = ii-1;
        A(ii, ii+2) = (k+2)*(k+1)/((2*k+3)*sqrt((2*k+1)*(2*k+5)))*c^2; 
    end
    
    for ii = 1:nmax-1
        k = ii-1; 
        A(ii+2, ii) = (k+2)*(k+1)/((2*k+3)*sqrt((2*k+1)*(2*k+5)))*c^2; 
    end
    
    Aeven = A(1:2:nmax+1, 1:2:nmax+1); 
    Aodd = A(2:2:nmax+1, 2:2:nmax+1); 
    
    % Assume the diagonals of Ds are in increasing order. 
    [Veven, ~] = eig(Aeven); 
    [Vodd, ~] = eig(Aodd); 
    
    
    % Compute lambda0. 
    
    y0 = evalPSWF(Pbar, Veven, Vodd, 0, x); 
    % y1 = evalPSWF(Pbar, Veven, Vodd, 1, x); 
    % y2 = evalPSWF(Pbar, Veven, Vodd, 2, x); 
    % y3 = evalPSWF(Pbar, Veven, Vodd, 3, x); 
    % 
    % figure(1); 
    % plot(x, y0); hold on; 
    % plot(x, y1); 
    % plot(x, y2);
    % plot(x, y3);
    
    yout = numericalBandlimitedOperatorApply(y0, c, x, h); 
    
    % plot(x, y0); hold on; grid on; 
    % plot(x, real(yout)); 
    
    lambda      = zeros(nmax+1, 1); 
    lambda(1)   = numericalInnerProd(y0, real(yout), h); 
    
    % Further compute lambdas. 
    yprev       = y0; 
    yderprev    = evalDerPSWF(Pbar, Veven, Vodd, 0, x); 

    for idx = 2:nmax+1
        ycur    = evalPSWF(Pbar, Veven, Vodd, idx-1, x); 
        ydercur = evalDerPSWF(Pbar, Veven, Vodd, idx-1, x); 
    
        z1       = numericalInnerProd(ycur, yderprev, h);
        z2       = numericalInnerProd(ydercur, yprev, h); 
        
        ratio   = z1/z2;
        lambda(idx) = 1i*lambda(idx-1)*sqrt(abs(ratio)); 
    
        yprev = ycur; 
        yderprev = ydercur; 
    end
    
    eigVals = c/(2*pi)*(abs(lambda).^2); 

%% Re-scale the eigen-functions 
    % assert(all(t_arr >= a)); 
    % assert(all(t_arr <= b)); 
    xq = (2*t_arr-(b+a))/(b-a); % Note: there may be values outside [-1,1]. 
    Nq = length(xq); 
    phi_q = zeros(Nmax, Nq); 

    for idx = 1:Nmax
        psi_val = evalPSWF(Pbar, Veven, Vodd, idx-1, x);
        
        % Apply operator Fc to extend the PSWF values outside [-1,1] 
        phi_val = h * psi_val * exp(1i*c*(x.')*xq); 
        phi_val = phi_val / lambda(idx); 

        % Scale the phi values to satisfy the L2 norm constraint on the
        % whole real line R. This is equivalent to setting energy in [-1,1]
        % to be eigVals(idx). 
        phi_q(idx, :) = real(phi_val) * sqrt(eigVals(idx)) * sqrt(2/(b-a)); 

        % TODO: Merge the previous two operations 
    end

end

function Pbar = getNormalizedLegendrePolynomials(n)

    Pbar = zeros(n+1); 

    Pbar(1,1) = 1;  % P0(x) = 1
    Pbar(2,2) = 1;  % P1(x) = x

    for ord = 2:n
        Pbar(ord+1, :) = (2*ord-1)/(ord)*[0, Pbar(ord, 1:n)] - (ord-1)/ord*Pbar(ord-1,:); 
    end

    for idx = 1:n+1
        Pbar(idx, :) = Pbar(idx, :)*sqrt(1/2+idx-1); 
    end

end

function y = evalPoly(p, x)
    ord = length(p)-1;
    y = zeros(size(x));
    t = ones(size(x));
    for k = 0:ord
        y = y + p(k+1)*t;
        t = t .* x; 
    end
end

function y = evalPSWF(Pbar, Veven, Vodd, ordPSWF, x)
    [nmax, ~] = size(Pbar); 
    nmax = nmax - 1;  

    y = zeros(size(x)); 
    
    if mod(ordPSWF, 2) == 0
        v = Veven(:, ordPSWF/2+1); 
        Ps = Pbar(1:2:end, :); 
    else
        v = Vodd(:, (ordPSWF+1)/2); 
        Ps = Pbar(2:2:end, :); 
    end

    for idx = 1:length(v)
        y = y + v(idx)*evalPoly(Ps(idx, :), x); 
    end

end

function y = evalDerPSWF(Pbar, Veven, Vodd, ordPSWF, x)
    [nmax, ~] = size(Pbar); 
    nmax = nmax - 1;  

    y = zeros(size(x)); 
    
    if mod(ordPSWF, 2) == 0
        v = Veven(:, ordPSWF/2+1); 
        Ps = Pbar(1:2:end, :); 
    else
        v = Vodd(:, (ordPSWF+1)/2); 
        Ps = Pbar(2:2:end, :); 
    end

    dPs = Ps(:,2:end).*(1:nmax); 

    for idx = 1:length(v)
        y = y + v(idx)*evalPoly(dPs(idx, :), x); 
    end

end

function I = numericalInnerProd(x, y, h)
    I = sum(x.*y)*h;
end

function y = numericalBandlimitedOperatorApply(phi, c, x, h)
    L = length(x); 
    y = zeros(1, L);

    for idx = 1:L
        y(idx) = numericalInnerProd(exp(1i*c*x(idx)*x), phi, h); 
    end

end




