% ============================================= % 
%
% This script aims to construct the PSWF codebook 
% and study their spatial/angular response 
% 
% ============================================= % 
clear all; close all; clc; 

fc  = 28e9; 
c   = 3e8; 
lambda = c/fc; 

angularOrder = 0; 
radialOrder = 0; 

load(sprintf('PSWFdata/N%d.mat', angularOrder)); 

Rfunc   = R0(radialOrder+1, :); 
rs      = r; 
thetas  = 0:ht:2*pi;  

% Meshgrid for computing the 2D PSWF-shaped beam 
% Compute the angular beam pattern using Fourier transform. 

% Compute Fourier transforms of psi(r, \theta)



[azimuthAngles, elevationAngles] = meshgrid(-pi/2:5*ht:pi/2, -pi/3:5*ht:pi/3); 

[N_ele, N_azi] = size(azimuthAngles); 
beamGain = zeros(N_ele, N_azi); 

for 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 = cos(theta_ele)*sin(phi_azi); 
        fz = sin(theta_ele); 
        
        % Compute int(\psi, exp(-i(y*fy + z*fz))) dA
        % dA = r dr dtheta
        tmp = 0; 
        for idx_theta = 1:length(thetas)
            y = rs*cos(thetas(idx_theta)); 
            z = rs*sin(thetas(idx_theta)); 

            tmp = tmp + hr*ht*sum(exp(-1i*(fy*y+fz*z)).*Rfunc.*rs)*exp(1i*angularOrder*2*pi*thetas(idx_theta)); 
        end
        beamGain(n_ele, n_azi) = tmp;  
    end
end

% Plot the 3D shape of this beam 

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


%% Visualizations 

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

