close all; clear all; clc;


hr      = 0.001;                        % numerical step
c       = 12;
r       = hr/2:hr:1-hr/2;               % Mid-point grid

ht      = 0.005;
theta   = 0:ht:2*pi;

% Compute for N=0.
N       = 0;
nmax    = 30;

[lambda, varphi_val] = getOrderN2DPSWFs(N, c, r, hr, nmax);

beta    = lambda/sqrt(c);
alpha   = 2*pi*beta;
R0      = varphi_val./sqrt(2*pi*r);     % orthonormal in D.
T       = ones(size(theta));
concentrationSeq = (c*abs(alpha)/(2*pi)).^2;

save('./N0.mat', 'concentrationSeq', 'alpha', 'R0', 'T', 'hr', 'ht', 'c', 'r', 'nmax');
% Compute for N>=1.

% Positive/negative angular frequency separate.

Nmax = 15;
concentrationSeqN = zeros(Nmax, nmax+1);
for N = 1:Nmax
    [lambda2, varphi2_val] = getOrderN2DPSWFs(N, c, r, hr, nmax);

    beta2   = lambda2/sqrt(c);
    alpha_p = 2*pi*exp(1i*pi/2*N)*beta2;
    alpha_n = 2*pi*exp(-1i*pi/2*N)*beta2;

    RN      = varphi2_val./sqrt(2*pi*r);
    TNp     = exp(1i*N*theta);
    TNn     = exp(-1i*N*theta);

    concentrationSeqN(N, :) = (c*abs(alpha_p)/(2*pi)).^2;
    save(sprintf('./N%d.mat', N), 'alpha_n', 'alpha_p', 'RN', 'TNp', 'TNn', 'hr', 'ht', 'c', 'r', 'nmax');
    fprintf('N=%d complete\n', N);
end

cseq = [concentrationSeq; concentrationSeqN];
save('./Concentrations.mat', 'cseq');

