close all; clear; clc;

%%
N = 100;
options = optimset('Display','none');
for j = 1:10
    R1 = rand; R2 = rand; % Random reflectivities
    c = rand*3; l = rand*2; % Random speed of light and cavity length
    I = NaN(1,N); I(1) = 1;
    for k = 2:100
        I(k) = I(k-1)*R1*R2;
    end
    t = (0:(numel(I)-1)).*(l/c);
    
    % Fit with an exponential
    f = @(tau) mean((exp(-t/tau)-I).^2);
    tau = fminsearch(f,1,options);
    
    % Plot
    figure(1); clf;
    hold on;
    plot(t,I,'o');
    plot(t,exp(-t/tau));
    %a = gca; a.YScale = 'log';
    grid on; box on; axis square;
    xlim([min(t) max(t)]);
    drawnow;
    
    % Print info
    fprintf('actual tau = %.5f\n',tau);
    fprintf(' fitted tau = %.5f\n',-1/log(R1*R2).*(l/c));
    fprintf('      delta = % .3g\n',(-1/log(R1*R2).*(l/c))-tau);
    fprintf('\n');
end