close all; clear; clc;
seed = 28; rng(seed);

%% Calculate hypersurfaces
fprintf('Calculating L1\n');
p1s = linspace(0,1,101); p2s = linspace(0,1,101);
L1 = NaN(numel(p1s),numel(p2s));
for j = 1:numel(p1s)
    for k = 1:numel(p2s)
        L1(j,k) = window_lossfun([p1s(j),p2s(k)]);
        %L3(a1i,a2i) = window_lossfun_2([a1s(a1i),a2s(a2i)]);
    end
    if ~mod(j,round(numel(p1s)/20)); fprintf('%02d%% done\n',round(j/numel(p1s)*100)), end
end

fprintf('Calculating L2 & L3\n');
b1s = linspace(0,pi/2,101);
b2s = linspace(0,1,102); b2s(1) = [];
L2 = NaN(numel(p1s),numel(p2s));
L3 = NaN(numel(p1s),numel(p2s));
for j = 1:numel(b1s)
    for k = 1:numel(b2s)
        temp = b2s(k).*exp(1i*b1s(j));
        L2(j,k) = window_lossfun([real(temp),imag(temp)]);
        L3(j,k) = window_lossfun_2([real(temp),imag(temp)]);
    end
    if ~mod(j,round(numel(p1s)/20)); fprintf('%02d%% done\n',round(j/numel(p1s)*100)), end
end
L1 = L1'; L2 = L2'; L3 = L3';

%% Plot
figure(1); clf;
axes; hold on;
colormap(viridis);

s = subplot(1,3,1);
contourf(p1s,p2s,L1,10,'edgecolor','none');
axis square;
%s.XTickLabel = ''; s.YTickLabel = '';
xlabel('p_1'); ylabel('p_2');
grid on;

s = subplot(1,3,2);
contourf(2*b1s/pi,b2s,L2,10,'edgecolor','none');
axis square;
%s.XTickLabel = ''; s.YTickLabel = '';
xlabel('2*\theta / \pi'); ylabel('r');
grid on;

s = subplot(1,3,3);
contourf(2*b1s/pi,b2s,L3,10,'edgecolor','none');
axis square;
%s.XTickLabel = ''; s.YTickLabel = '';
xlabel('2*\theta / \pi'); ylabel('r');
grid on;