clear;

%% create a radiator plane in the XY plane
W = 1;                  % width [m]
d = .2e-2;              % element size [m]
ds = d^2;               % element are [m^2]
x = -W/2 : d : W/2;
y = -W/2 : d : W/2;
[Y, X] = meshgrid(y, x);    % center coordinates

%% create a semicircle in the XZ plane
phivec = -pi/2 : pi/180 : pi/2;
N = length(phivec);
R = 1;
x_h = R * sin(phivec);
y_h = zeros(1, N);
z_h = R * cos(phivec);


%% create a semicircle in the YZ plane
y_v = R * sin(phivec);
x_v = zeros(1, N);
z_v = R * cos(phivec);

%% create a velocity excitation on the plane
v = sqrt(X.^2 + Y.^2) < .3;     % select elements that are close to center

%% plot the system
figure;
hold on;
pcolor(X, Y, double(v));
shading flat;
plot3(x_h, y_h, z_h, '.red');
plot3(x_v, y_v, z_v, '.blue');
axis equal;
xlabel('x');
ylabel('y');
zlabel('z');

%% analysis parameters
freq = 1e3;     % [Hz]
om = 2*pi*freq; % angular frequency [rad/s]
c = 340;        % [m/s]
rho = 1.2;      % [kg/m^3]
k = om/c;       % wave number [rad/m]

%% compute Rayleigh integral
dirchar_h = nan(1,N);
dirchar_v = nan(1,N);
for i = 1 : N
    % horizontal
    r = sqrt((x_h(i) - X).^2 + (y_h(i) - Y).^2 + z_h(i).^2);
    integrand = (-1i*om*rho/2/pi)*v.*exp(-1i * k * r) ./ r;
    dirchar_h(i) = sum(sum(integrand * ds));
    % vertical
    r = sqrt((x_v(i) - X).^2 + (y_v(i) - Y).^2 + z_v(i).^2);
    integrand = (-1i*om*rho/2/pi)*v.*exp(-1i * k * r) ./ r;
    dirchar_v(i) = sum(sum(integrand * ds));
end

%% plot direction characteristics
figure;
plot(phivec, 20*log10(abs(dirchar_h)), ...
    phivec, 20*log10(abs(dirchar_v)));
legend('horizontal', 'vertical');
xlabel('Angle [rad]');
ylabel('Directivity [dB]');
