Snells_law_wavefronts.gif
No higher resolution available.
Titel
Summary
This diagram was created with
MATLAB
.
| Description | ||||
| Date | (UTC) | |||
| Source | Own work, created with MATLAB . ( Source code below.) | |||
| Author | Oleg Alexandrov | |||
|
Permission
( Reusing this file ) |
|
|||
| Other versions |
|
Source code
% Illustration of Snell's law
function main()
% indexes of refraction
n1=1.0;
n2=1.5;
sign = -1;% is the source up or down?
O=[0, -1*sign];
k=500;
% KSmrq's colors
red = [0.867 0.06 0.14];
blue = [0, 129, 205]/256;
green = [0, 200, 70]/256;
yellow = [254, 194, 0]/256;
white = 0.99*[1, 1, 1];
black = [0, 0, 0];
gray = 0.5*white;
color1=red;
color2=blue;
color3=gray;
lw = 3;
plot_line=0;
Theta=linspace(0, 2*pi, k);
V=0*Theta; W=0*Theta;
S0=7;
spacing=0.45;
p=floor(S0/spacing);
S=linspace(0, S0, p+1);
spacing=S(2)-S(1);
num_frames = 10;
for frame_iter=1:num_frames
figure(1); clf; hold on; axis equal; axis off;
% plot the interface between diellectrics
L=1.2*S0;
plot([-L, L], [0, 0], 'color', color3, 'linewidth', lw);
% plot a ray
plot_line=1;
s=L;
theta=pi/3; wfr(s, theta, n1, n2, O, sign, plot_line, color1, lw);
% plot the wafefronts
plot_line=0;
for i=1:p
s=S(i)+spacing*(frame_iter-1)/num_frames;
for j=1:k
theta=Theta(j);
[V(j), W(j)]=wfr(s, theta, n1, n2, O, sign, plot_line, color1, lw);
end
plot(V, W, 'color', color2, 'linewidth', lw);
end
% dummy points to enlarge the bounding box
plot(0, S0+2.5*spacing, '*', 'color', white);
plot(0, -(S0+2.5*spacing)/n2, '*', 'color', white);
% to know where to crop later
Lx=3.2; Ly=Lx; shift = 1;
plot([-Lx, Lx, Lx, -Lx -Lx], ...
[-Ly, -Ly, Ly, Ly, -Ly]+shift);
file = sprintf('Frame%d.eps', 1000+frame_iter);
disp(file);
saveas(gcf, file, 'psc2')
end
end
% Converted to gif with the UNIX command
% convert -density 100 -antialias Frame10* Snell_animation.gif
% then cropped in Gimp
function [a, b]=wfr(s, theta, n1, n2, O, sign, plot_line, color1, lw);
X=O+s*[sin(theta), sign*cos(theta)];
if( sign*X(2) > 0 )
t=-sign*O(2)/cos(theta);
X0=O+t*[sin(theta), sign*cos(theta)];
if (plot_line == 1)
plot([O(1), X0(1)], [O(2), X0(2)], 'color', color1, 'linewidth', lw, 'linestyle', '--');
end
d = norm(O-X0);
r = (s-d)*(n2/n1)^(sign);
theta2=asin(n1*sin(theta)/n2);
XE=X0+r*[sin(theta2), sign*cos(theta2)];
else
XE = X;
end
a = XE(1);
b = XE(2);
if (plot_line==1)
plot([X0(1), XE(1)], [X0(2), XE(2)], 'color', color1, 'linewidth', lw, 'linestyle', '--');
end
end