PHYSICAL OPTICS
Maple Documentation
Maple Source Code

This problem has light rays coming parallel to the axis and refracting inside a sphere then coming back out. The idea is to show that for a range of index of refraction one will get rays coming back out which are mostly parallel to the axis, so making a spherical bead a good approximation to a retro-reflector.

Perry shows n=1.85 below, but one would want to have a way to compare n=1.85 to other indices of refraction for effectiveness.

restart;
n:=1.85;

Position of horizontal lines.

for i from 0.05 to 1 by 0.05 do


Equation of incoming ray.
L1[i]:=i;

Find where incoming ray hits the circle.

sol[i]:=fsolve({y1[i]=L1[i],x1[i]=cos(theta1[i]),y1[i]=sin(theta1[i])},{theta1[i],x1[i],y1[i]},{theta1[i]=-Pi/2..Pi/2}); assign(sol[i]);

Find the angle of the refracted ray.

snell[i]:=1*sin(theta1[i])=n*sin(theta2[i]);
theta2[i]:=solve(snell[i],theta2[i]);

Find the slope angle of the refracted ray.
thetaL2[i]:=theta1[i]-theta2[i];

Equation of the line for the refracted ray.
L2[i]:=tan(thetaL2[i])*(x-x1[i])+y1[i];

Find where the refracted ray intersects the circle. sol[i]:=fsolve({y2[i]=subs(x=x2[i],L2[i]),x2[i]=cos(theta3[i]),y2[i]=sin(theta3[i])},{theta3[i],x2[i],y2[i]},{theta3[i]=Pi/2..3*Pi/2}); > assign(sol[i]);

Find the slope angle of the reflected ray.
thetaL3[i]:=thetaL2[i]-2*theta2[i];

Equation of the line for the reflected ray.
L3[i]:=tan(thetaL3[i])*(x-x2[i])+y2[i];

Find where the reflected ray intersects the circle.
sol[i]:=fsolve({y3[i]=subs(x=x3[i],L3[i]),x3[i]=cos(theta4[i]),y3[i]=sin(theta4[i])},{theta4[i],x3[i],y3[i]},{theta4[i]=Pi..2*Pi}); > assign(sol[i]);

Find the slope angle of the reflected ray leaving the circle.
thetaL4[i]:=theta4[i]+theta1[i];

Equation of the line for the ray leaving the circle.
L4[i]:=tan(thetaL4[i])*(x-x3[i])+y3[i]; > od:

Make plots for the incomming and outgoing rays.
for i from 0.05 to 1 by 0.05 do
f.(trunc(i*20)):=plot({[x,L1[i],x=x1[i]..2],[x,L2[i],x=x2[i]..x1[i]]},color=blue,scaling=constrained): h.(trunc(i*20)):=plot({[x,L3[i],x=x2[i]..x3[i]],[x,L4[i],x=x3[i]..2]},color=green,scaling=constrained):
od:

Make a plot of the circle.
g:=plot([cos(theta),sin(theta),theta=0..2*Pi],color=black):

Display all the plots on the same graph.
with(plots):
display({g,f.(1..20),h.(1..20)});

restart;
with(plots):

Equation for the parabola
y1:=x^2/(2*R);

Find the slope of the parabola
dy_dx:=diff(y1,x);

Find the angle of this slope makes with the x axis
theta:=arctan(dy_dx);

theta := arctan(x/R)

The angle the reflected ray makes with the x axis is:
theta2:=2*theta-Pi/2;

Equation of reflected ray
y2:=tan(theta2)*(x2-x)+y1;

Find where the ray hits the parabola
x3:=min(solve(y2=subs(x=x2,y1),x2));
y3:=subs(x=x3,y1);

Radius of curvature of the mirror (try values between 5 and 100)
R:=20;
x:='x';

Draw the mirror
g17:=plot([x,y1,x=-35..35],scaling=constrained,color=black):

Draw the rays
xmax:=subs(x=32,y1);
for x from 1 to 31 by 2 do g.(x/2+1/2):=plot([[x,max(2*R,xmax)],[x,y1],[0,y1-tan(theta2)*x],[x3,y3],[x3,max(2*R,xmax)]],X=-35..35,Y=0..max(2*R,xmax),color=COLOR(RGB, evalf(sin(x/8)),evalf(sin(x/8-2*Pi/6)),evalf(sin(x/8-4*Pi/6))));
od:

Draw the graph and axis
R1:=R:
R:='R': g0:=PLOT(CURVES([[2,R1],[-2,R1]]),CURVES([[2,R1/2],[-2,R1/2]]),CURVES([[0,0],[0,max(2*R1,xmax)]]),TEXT([2,R1+2],R),
TEXT([2,R1/2+2],f),AXESSTYLE(FRAME),TITLE(Parabolic.` `.Mirror));

display([g.(0..17)]);

restart;
with(plots):

Radius of curvature of the mirror (best results obtained when greater than 50)
R:=50;

Height of the incoming ray from the optical axis
y:='y';

Compute DeltaX (distance along the optical axis from where the ray hits the mirror to the center of the mirror)
DeltaX:=R-sqrt(R^2-y^2);

Compute R' (radius of the mirror minus DeltaX)
R_p:= R - DeltaX;
theta:=arcsin(y/R);
phi:=arccos(y/R);
theta_p:=theta-phi;
xp:= y*tan(theta_p);

Draw the graph and axis
g0:=PLOT(CURVES([[0,2],[0,-2]]),CURVES([[R/2,2],[R/2,-2]]),CURVES([[-R,0],[R,0]]),TEXT([-2,2],r),TEXT([R/2-2,2],f),
AXESSTYLE(FRAME),TITLE(Spherical.` `.Mirror));

Draw the mirror
g1:=plot([R*cos(th),R*sin(th),th=-Pi/2..Pi/2],scaling=constrained,color=black):

Draw the rays
for y from 4 to 32 by 2 do
g.(y/2):=plot([[-100,y],[R_p,y],[R_p+xp,0],[R_p+6*xp,-5*y]],
x=-R..R,y=-34..34,color=COLOR(RGB, evalf(sin(y/8)),evalf(sin(y/8-2*Pi/6)),evalf(sin(y/8-4*Pi/6))));
od:

display([g.(0..16)]);

 Back to Physical Optics