function [lambda, varphi_val] = getOrderN2DPSWFs(N, c, r, hr, nmax)
    
    P = getJacobi(N, 0, nmax);                  % Evaluate Jacobi polynomial coefficients up to order nmax. 
    [Tval, hval] = getT(N, nmax, r, P);         % Evaluate T functions on grid r, together with normalizing constants hval. 
    BN = getB(N, hval, c, nmax);                % Evaluate B matrix. 
    
    % Compute eigenfunctions with expansion in Jacobi polynomials. 
    % for each integer N>=0, we consider the angular form of exp(iNθ). 
    % Then the radial integral equation has eigensystem indexed by n>=0. 
    
    
    [V, D] = eig(BN);
    chi = -diag(D);
    d = V.';      % each row: coefficients of PSWF(N,n) in T. 
    
    
    varphi_val = get_varphi(N, nmax, r, Tval, d); 
    
    lambda = zeros(1, nmax+1);
    intVal = zeros(1, length(r)); 
    
    for nidx = 1:nmax+1
    
        for idx = 1:length(r)
            intVal(idx) = numericalInnerProd01(besselj(N, c*r(idx)*r).*sqrt(c*r(idx)*r), varphi_val(nidx, :), hr); 
        end
        
        lambda(nidx) = numericalInnerProd01(intVal, varphi_val(nidx,:), hr);
    end


end


function P = getJacobi(a, b, n)
% Return a matrix P of size (n+1)*(n+1), containing Jacobi polynomials of
% order 0~n. 

P = zeros(n+1, n+1);

P(1,1) = 1; 
P(2,1) = (a-b)/2;
P(2,2) = (a+b+2)/2;

for ord = 2:n
    a1 = 2*ord*(ord+a+b)*(2*(ord-1)+a+b);
    a2 = -(2*(ord-1)+a+b+1)*(a^2-b^2);
    a3 = (2*(ord-1)+a+b)*(2*(ord-1)+a+b+1)*(2*(ord-1)+a+b+2); 
    a4 = 2*((ord-1)+a)*((ord-1)+b)*(2*ord+a+b); 

    P(ord+1, :) = (-a2*P(ord, :)-a4*P(ord-1, :)+a3*[0, P(ord,1:n)])/a1; 
end


end

function [Tval, hval] = getT(N, n, x, P)
% Assume P is evaluated at (N, 0). 
% up to T_{N,n}-th order. 
    L = length(x);
    Tval = zeros(n+1, L);
    hval = zeros(n+1, 1); 
    
    for nidx = 0:n
        C = nchoosek(nidx+N, nidx);
        h = sqrt(2*(2*nidx+N+1))*C;
        
        t = 1-2*x.^2; 
        y = evalPoly(P(nidx+1, :), t); 
        Tval(nidx+1, :) = h*(x.^(N+1/2)).*(y)/C;
        hval(nidx+1) = h; 
    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 drawPoly(p, x)
    y = evalPoly(p, x); 
    plot(x, y);
end


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

function BN = getB(N, hval, c, n)
% Expand the PSWFs with orthonormal basis T. 
    BN = zeros(n+1);
    
    for idx = 0:n
        kappa = (N+2*idx+1/2)*(N+2*idx+3/2); 

        if N==0 && idx == 0
            gamma0 = 1/2;
        else
            gamma0 = (2*idx*(idx+1)+N*(2*idx+N+1))/(2*idx+N)/(2*idx+N+2);
        end

        BN(idx+1, idx+1) = -(kappa + c^2*gamma0); 
        
        if idx <= n-1
            gamma_p1 = -(idx+N+1)^2/(2*idx+N+1)/(2*idx+N+2)*(hval(idx+1)/hval(idx+2)); 
            BN(idx+2, idx+1) = -c^2*gamma_p1; 
        end

        if idx >= 1
            gamma_n1 = -(idx^2/(2*idx+N)/(2*idx+N+1))*(hval(idx+1)/hval(idx)); 
            BN(idx, idx+1) = -c^2*gamma_n1; 
        end
        
    end

end

function varphi_val = get_varphi(N, n, r, Tval, d)
% evaluate up to order n. 
    L = length(r);
    varphi_val = zeros(n+1, L);

    for idx = 0:n
        varphi_val(idx+1, :) = d(idx+1, :)*Tval; 
    end

end


