% Run this script either after HRIR_Berechnung or after loading one of the
% results files that were saved by HRIR_Berechnung

% Berechnung aller ITDs
ILDs = nan(numel(azimuths),numel(elevations));
for i = 1:numel(azimuths)
    for j = 1:numel(elevations)
        hrtf_l = fft(squeeze(hrir_l(i,j,:)));
        hrtf_r = fft(squeeze(hrir_r(i,j,:)));
        sum_l = abs(hrtf_l)' * abs(hrtf_l);
        sum_r = abs(hrtf_r)' * abs(hrtf_r);
        ILDs(i,j) = 10 * log10(sum_r) - 10 * log10(sum_l);
    end
end

figure
hold on
contourf(azimuths, elevations, ILDs', -20:1:20)
contour(azimuths, elevations, ILDs', -40:2:40, 'LineWidth', 3, 'Color', 'k')
contour(azimuths, elevations, ILDs', -1000:1000:1000, 'LineWidth', 5, 'Color', 'r', 'linestyle', ':')
[~, max_I] = max(ILDs(:));
[max_I_row, max_I_col] = ind2sub(size(ILDs), max_I);
line([0, 0] + azimuths(max_I_row), [-3, +3] + elevations(max_I_col), 'linewidth', 1, 'color', 'r')
line([-3, +3] + azimuths(max_I_row), [0, 0] + elevations(max_I_col), 'linewidth', 1, 'color', 'r')
viscircles([azimuths(max_I_row), elevations(max_I_col)], 3, 'linewidth', 1);
% plot(azimuths(max_I_row), elevations(max_I_col), '+r');
% plot(azimuths(max_I_row), elevations(max_I_col), 'or');
[~, min_I] = min(ILDs(:));
[min_I_row, min_I_col] = ind2sub(size(ILDs), min_I);
line([-3, +3] + azimuths(min_I_row), [0, 0] + elevations(min_I_col), 'linewidth', 1, 'color', 'r')
viscircles([azimuths(min_I_row), elevations(min_I_col)], 3, 'linewidth', 1);
colorbar
set(gca, 'xtick', (-180:30:180))
set(gca, 'ytick', (-180:30:180))
set(gca, 'xlim', [-180 180])
set(gca, 'ylim', [-90 90])
set(gca, 'DataAspectRatio', [1 1 1])
grid
title('ILD')