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

% Berechnung aller ITDs
ITDs = nan(numel(azimuths),numel(elevations));
for i = 1:numel(azimuths)
    for j = 1:numel(elevations)
        xc = xcorr( squeeze(hrir_l(i,j,:)), squeeze(hrir_r(i,j,:)), 'coeff');
        [corr_max,g] = max(xc); 
        ITDcorrsamples = size(hrir_l,3)-g;
        ITDcorr = ITDcorrsamples/samplingrate;
        ITDs(i,j) = -ITDcorr;
    end
end

figure
hold on
contourf(azimuths, elevations, 1e6*ITDs', -300:10:300)
contour(azimuths, elevations, 1e6*ITDs', -1000:50:1000, 'LineWidth', 3, 'Color', 'k')
contour(azimuths, elevations, 1e6*ITDs', -1000:1000:1000, 'LineWidth', 5, 'Color', 'r', 'linestyle', ':')
[~, max_I] = max(ITDs(:));
[max_I_row, max_I_col] = ind2sub(size(ITDs), 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(ITDs(:));
[min_I_row, min_I_col] = ind2sub(size(ITDs), 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('ITD')