close all; clear; clc;

%% Make random mirror and starting points, and show the results

odefun = @(t,y) [y(2);0];
options = odeset('RelTol',1e-6,'AbsTol',1e-6,'MaxStep',0.01);
K = exp(1i*linspace(0,2*pi,1e2+1))*1.5; K(end) = [];
N = 100;
randpm = @() 2*rand-1; % Uniform random number in the +/-1 range

c = 1; % Speed of light
v = 1+0.*1i; v = v./norm(v)*c; % Initial velocity

while true
    y0 = [0+0.5*randpm().*1i v]; % Starting position
    M1 = resonator_alignment_make_mirror(2*randpm(),1);
    M2 = resonator_alignment_make_mirror(2*randpm(),2);
    
    % Rotate M1 by angle theta
    % For large angles, no stable resonator configuration will occur
    theta = 0; % 0.1;
    M1 = ((M1+1).*exp(1i*(theta))-1);
    
    options.Events = @(t,p) resonator_alignment_eventfun(t,p,M1,M2);
    
    Y = y0(1); T = 0;
    ystart = y0;
    for j = 1:N
        [t,y] = ode23(odefun,[0 10],ystart,options);
        ye = y(end,:);
        ystart = resonator_alignment_impactmap(ye,M1,M2);
        
        Y = [Y; y(end,1)]; T = [T; t(end)+T(end)]; %#ok<AGROW>
        if abs(ye(1))>=1.5 % Escaped the cavity
            break
        end
    end
    
    if j > 2 % Animate the result if there were at least two reflections
        figure(1); clf;
        plot(real(y0(1)),imag(y0(1)),'ko'); hold on;
        plot(M1,'k-','linewidth',2); plot(M2,'k-','linewidth',2);
        plot(K,'k--'); % Plot killplane
        axis equal; grid on;
        xlim([-1 1]*2); ylim([-1 1]*2);
        
        cm = viridis(numel(T));
        for ti = 2:numel(T)
            plot(Y(ti-1:ti),'color',cm(ti,:),'linewidth',1);
            if ~mod(ti,round(numel(T)/20)); drawnow; end
        end
    end
end