close all; clear; clc;

%% Show that Snell's law minimises the travel time between two points
c = 1; n0 = 1; n1 = 3;

for j = 1:10
    
    % Pick start and end points
    x0 = 2*rand-1 + rand*1i; x1 = 2*rand-1 - rand*1i;
    %x0 = 2*rand-1 + 1i; x1 = 2*rand-1 - 1i;
    %x1 = 0 - 1i;
    
    xi = @(angle) norm(x0)/tan(pi/2-angle)+real(x0);
    T = @(angle) norm(xi(angle)-x0)/(c/n0) + norm(x1-xi(angle))/(c/n1);
    alpha = fminsearch(T,0);
    
    figure(1); clf
    axes; hold on;
    plot(complex(0),'k+');
    plot([-1 1],[0 0],'k--');
    plot([x0 xi(alpha)],'-o');
    plot([xi(alpha) x1],'-o');
    plot([x0 x1],'k:');
    axis equal square;
    xlim([-1 1]); ylim([-1 1]);
    grid on; box on;
    
    theta0 = pi/2-angle(x0-xi(alpha));
    theta1 = pi/2-angle(xi(alpha)-x1);
    
    fprintf('n0*sin(theta0) = %.5f\n', n0*sin(theta0));
    fprintf('n1*sin(theta1) = %.5f\n', n1*sin(theta1));
    drawnow; pause(0.1);
end

%% Repeat for arbitrary surface shape
c = 1; n0 = 1; n1 = 3;

% Pick a surface shape
%m = @(x) zeros(size(x));
%m = @(x) x.^2;
m = @(x) 0.2*cos((1.5)*pi*x);

for j = 1:10
    % Pick start and end points
    p0 = 2*rand-1 + 1i; p2 = 2*rand-1 - 1i;
    
    x1 = @(angle) norm(p0)/tan(pi/2-angle)+real(p0);
    p1 = @(angle) x1(angle) + 1i*m(x1(angle));
    T = @(angle) norm(p1(angle)-p0)/(c/n0) + norm(p2-p1(angle))/(c/n1);
    alpha = fminsearch(T,0);
    
    figure(2); clf
    axes; hold on;
    x = linspace(-1,1,1e3);
    plot(x,m(x),'k--')
    plot(complex(0),'k+');
    %plot([x0 0],'-o');
    %plot([x1 0],'-o');
    %plot(x0-(x0.*exp(1i*theta)),'o');
    plot([p0 p1(alpha)],'-o');
    plot([p1(alpha) p2],'-o');
    plot([p0 p2],'k:');
    axis equal square;
    xlim([-1 1]); ylim([-1 1]);
    grid on; box on;
    drawnow; pause(0.1)
end