close all; clear; clc;

%%
Rs = fliplr([0.9 0.8 0.5 0.2 0.1]);
l = 1;
for Ri = 1:numel(Rs)
    R = Rs(Ri); % Cavity length and reflectivity
    fc = 5/l; % 5 oscillations per cavity length
    fsr = 1/(2*l); % Free spectral range
    N = max(20,round(log(0.001)/log(R))); % Number of reflections
    X = bsxfun(@plus,linspace(0,l,1e3),2*l*(0:(N-1))'); % x's for all reflections
    fs = linspace((fc-1.5*fsr),(fc+1.5*fsr),1e3); % Frequencies
    Y = sin(bsxfun(@times,X,2*pi*permute(fs,[3,1,2]))); % All the reflections
    Y = bsxfun(@times,Y,R.^(0:(N-1))'); % Loss of R on each reflection
    % Peak-to-peak amplitude of the sum(Y) - numerical Airy distribution
    airy_calc(Ri,:) = (1-R)^2*squeeze(2*mean(sum(Y).^2,2));
    frel = (fs-fc)/fsr; % Frequency relative to the central frequency
    airy(Ri,:) = ((1-R)^2)./((1-R)^2+4*R*sin(pi*frel).^2); % Analytical Airy distr.
end
%% Plot
figure(1); clf;
subplot(1,2,2); hold on;
cm = flipud(viridis(numel(Rs)+1));
for j = numel(Rs):-1:1
    color = cm(numel(Rs)-j+1,:);
    plot(frel,airy_calc(j,:),'-','color',color,'linewidth',1);
end
ylim([0 1]);
a = gca; a.YColor = 'k';

lambda = 1./fs;
distance = @(x) abs(x-round(x));
h = plot(frel,distance(2*l./lambda),'--','color',[1 1 1]*0.3);
uistack(h,'bottom');
I = c2(frel,-1.15);
plot(frel(I),airy(:,I),'ro');
a = gca; a.YColor = 'k';

grid on; box on; axis square;
xlim([-1.5 1.5]);
xlabel('(fs-fc)/fsr');
ylabel('T');

a = subplot(1,2,1); hold on;
D = Y(1:10,:,I)';
cm = lines;
x = linspace(0,l,1e3);
for j = 1:size(D,2)
    if mod(j,2)
        color = cm(1,:); region = x<0.98;
    else
        color = cm(2,:); region = x>0.02;
    end
    xx = x(region);
    plot(xx,D(region,j)-2*(j-1),'color',color,'linewidth',1);
    patch([xx(1) xx xx(end)],[0 D(region,j)' 0]-2*(j-1),color,'edgecolor','none','facealpha',0.2);
end
plot([0 0],[1 -19],'k-','linewidth',1.5); plot([1 1],[1 -19],'k-','linewidth',1.5);
xlim([-0.1 l+0.1]); ylim([-20 2]);
a.YTick = -18:2:0;
a.YTickLabel = 9:-1:0;
ylabel('Reflections'); xlabel('l');
grid on; box on; axis square;