Tim Waters

clear all; %clears all variables
close all; %closes all figures

eps = .0001 %can also customize; input('Enter the tick-size (epsilon) ');
n1 = 1 %input('Enter the index of refraction in region containing A ');
n2 = 1.5 %input('Enter the index of refraction in region containing B ');
mirror_axis = 0:eps:1; %cut (0,1) into grid points
c = mirror_axis'; %creates a column vector
a = ones(length(c),1)*rand; %ones(rows,columns) returns a matrix of 1s
b = ones(length(c),1)*rand; %rand gives a random scalar in (0,1)

%Compute the optical path length, opl
d_ac = sqrt(a.^2 + c.^2);
d_cb = sqrt((1-c).^2 + b.^2);
opl_ac=n1*d_ac;
opl_cb=n2*d_cb;
opl = opl_ac + opl_cb;

opl_min = min(opl);
i.html">i=1;
while opl(i.html">i)~= opl_min
    i.html">i = i.html">i+1;
end
% i now contains the location of opl_min in the vector opl
%************************plot results *************************************
m = a(1)+b(1); %slope of straiht line path connecting a to -b
plot(c,zeros(size(c)),'r:') %plots interface axis y=0

hold on; %allows you to add more things to the plot

y_ab = -m*c + a(1);
plot(c,y_ab) %plots straight line path

m_a = a(1)/c(i.html">i);
y_ac = -m_a*c(1:i.html">i)+ a(1);
plot(c(1:i.html">i),y_ac) %plots computed shortest path from A to C

m_b = b(1)/(1-c(i.html">i));
y_cb = -m_b*c(i.html">i:length(c)) + m_b*c(i.html">i);
plot(c(i.html">i:length(c)),y_cb) %plots computed shortest path from A to C

text(0,a(1),'Point A'); %labels via text(x,y,'string')
text(1,-b(1),'Point -B');
plot(c(i.html">i),0,'kx') %plots point c with a black x over it
text(c(i.html">i),0,sprintf('Point C = %0.5g ',c(i.html">i)));

d_ab = sqrt(1+(b(1)+a(1))^2); %dist between pts A and -B
i_ab=floor((a(1)/m)/eps); %location of straight line path
title({'Straight Line Path and Least Time Path for n1/n2 interface';
   sprintf('Least time path dist = %0.5g ',(d_ac(i.html">i) + d_cb(i.html">i)));
  sprintf('Straight line path dist = %0.5g ',d_ab);
  sprintf('Least time path opl = %0.5g ',opl_min);
  sprintf('Straight line path opl = %0.5g ',opl(i_ab))})

%Check Answer with Snell's Law
n1/n2
c_ratio=(1-c(i.html">i))/c(i.html">i);
answr = c_ratio*d_ac(i.html">i)/d_cb(i.html">i);
xlabel({'Snells Law says n_i/n_f = (1-c)d_{ac}/(c*d_{cb})';
    sprintf('Thus compare n_i/n_f = %0.5g ',n1/n2);
    sprintf('with (1-c)d_ac/(c*d_cb) = %0.5g ',answr)});
hold off