%% Add directories containing support functions.
addpath('functions')

%% Inputs
d = 11.43e-3;     % Distance between microphones.
N = 12;            % Number of microphones.
f_max = 40e3;         % Frequency of source. Might be multiple frequencies?? Model non-ideal source???
phi = d2r(90);       % Look-direction of beam.
psi = linspace(0, pi, 1000);
c = 343;          % Velocity of sound (in air).

%% Frequency response surface plot
% Resolution
angleStepSize = 0.1; % deg
freqStepSize  = 100;

% output = f(theta, freq)
[t,f] = meshgrid(deg2rad(-90:angleStepSize:90), 0:freqStepSize:f_max);

% Initialize surface plot array
output = zeros(size(f));

% Can this be done better?
for sensor_idx = 0:N-1
    output = output + exp(j*2*pi.*f*sensor_idx*d.*sin(t)./c); % Computed using complex values directly.
end

% Convert to magnitude [dB]
output = (20*log10(abs(1/N.*output)));
% Limit all negligable gain values.
output(output < -30) = -30; 

%cont = contourf(f,t,output)
subplot(1,3,1)
surf(f,t,output, 'EdgeColor','interp', 'FaceColor', 'interp', 'FaceLighting', 'gouraud')
colorbar
view(-90,90) % Default view of surface plot
title('Frequency response')
axis square

% Cross section of surface plot obtained by fixing frequency to a constant
% value ie.
subplot(132)
ca = beamPlotCart(t(1,:), output(251,:))
axis square

paHandle = subplot(1,3,3,polaraxes)
pa = beamPlotPol(paHandle, t(1,:), output(251,:))

%% Surface plots in grid formation

% Parameter combinations to compute/plot
spacing_vec = [d 2*d 3*d ];
sensors_vec = [ 8  16  24 ];

figure
i = 1;
for spacing = spacing_vec
    for sensors = sensors_vec
        subplot(length(spacing_vec),length(sensors_vec),i);
        %%
        [t,f] = meshgrid(deg2rad(-90:angleStepSize:90), 0:freqStepSize:f_max);
        output = zeros(size(f));
        for sensor_idx = 0:sensors-1
            output = output + exp(1i*2*pi.*f*sensor_idx*spacing.*sin(t)./c); % Computed using complex values directly.
        end
        output = (20*log10(abs(1/sensors.*output)));
        output(output < -30) = -30; 

        surf(f,t,output, 'EdgeColor','interp', 'FaceColor', 'interp', 'FaceLighting', 'gouraud')
        view(-90,90)
        %axis square
        ylabel('Incident angle \psi'); xlabel('Frequency'); zlabel('Gain [dB]');
        title([num2str(sensors) ' sensors, ' num2str(spacing) 'm spacing'])

        i = i+1
    end
end

%camproj('orthographic')
%delay = sum(i*l.*sin(-90:90)./c)