% ========================================== %
%
% This file is to study the beamforming properties of the 2D-PSWFs.
%
% ========================================== %

clear; clc; close all;

c0 = 3e8;
fc = 3.5e9;
Na = 8;
lambda = c0/fc;
Lt = Na*(lambda/2);
At = pi*(Lt/2)^2;
kc = 2*pi*fc/c0;

Grid.hr     = 0.001;
Grid.r      = Grid.hr/2:Grid.hr:1-Grid.hr/2;    % Mid-point grid
Grid.theta  = (0:1000-1)/1000*(2*pi);           % Gauss quadrature points
Grid.ht     = Grid.theta(2);
[rs, thetas] = meshgrid(Grid.r, Grid.theta);
Grid.rs     = rs;
Grid.thetas = thetas;



c       = 2*pi*(0.2*kc)*Lt;

% Compute the first nmax radial eigen-functions for radial order N.
N       = 0;
nmax    = 30;   % spectral order
[lambda2, varphi2_val] = getOrderN2DPSWFs(N, c, Grid.r, Grid.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*Grid.r);
TNp     = exp(1i*N*Grid.theta);
TNn     = exp(-1i*N*Grid.theta);
concentration = (c*abs(alpha_p)/(2*pi)).^2;


% stem(concentration);

%% Compute the beam pattern of this thing
radialOrder = 1;

% DoF estimation
% estDoF = (pi*(0.2*kc)^2*pi*Lt^2)/(2*pi)^2; % This may be incorrect.
[azimuthAngles, elevationAngles] = meshgrid(-pi/2.1:10*Grid.ht:pi/2.1, -pi/2.1:10*Grid.ht:pi/2.1);
[N_ele, N_azi] = size(azimuthAngles);
beamGain = zeros(N_ele, N_azi);

parfor n_azi = 1:N_azi
    for n_ele = 1:N_ele
        phi_azi = azimuthAngles(n_ele, n_azi);
        theta_ele = elevationAngles(n_ele, n_azi);
        fy = kc*cos(theta_ele)*sin(phi_azi);
        fz = kc*sin(theta_ele);

        % Compute int(\psi, exp(-i(y*fy + z*fz))) dA
        % dA = r dr dtheta

        ys = (Lt/2)* (Grid.rs.*cos(Grid.thetas));
        zs = (Lt/2)* (Grid.rs.*sin(Grid.thetas));

        tmp =  (exp(-1i*(fy*ys+fz*zs)).*RN(radialOrder+1, :)).*(Grid.rs).*exp(1i*N*2*pi*Grid.thetas);

        beamGain(n_ele, n_azi) = (Grid.hr*Grid.ht)*sum(tmp, 'all');
    end
    fprintf('%d/%d\n', n_azi, N_azi);
end

% Plot the 3D shape of this beam

[X, Y, Z] = sph2cart(azimuthAngles, elevationAngles, abs(beamGain));


%% Visualizations

fig2 = figure(2);
surf(X, Y, Z, abs(beamGain));
shading interp;
colorbar;
xlabel('X');
ylabel('Y');
zlabel('Z');
axis equal;
title('3D波束图');

