%!PS-Adobe-2.0
%%BoundingBox: 0 0 160 160
%lenses.inc
/blue-index 1.3440 def
/red-index 1.3309 def

% A B C -> [ ], [ r1 (= r2) ], [ r1 < r2 ] (for Ax^2 + B + C = 0)
% 0, 1, or 2 roots
/roots { 8 dict begin
aload pop
/C exch def
/B exch def
/A exch def
/D B dup mul 4 A mul C mul sub def
D 0 lt {
	[ ] 
}{
D 0 eq {
[ B neg A div 2 div ]
}{
% D > 0
  B 0 le {
	/r B neg D sqrt add A div 2 div def
  }{ 
  % B > 0
	/r B neg D sqrt sub A div 2 div def
  } ifelse
  /R C A div r div def
  R r gt {
    [r R]
  }{
    [R r]
  } ifelse
} ifelse
} ifelse
end } def

% x -> asin in degrees
/asin { 1 dict begin
/x exch def
x
1 x dup mul sub sqrt 
atan
end } def

% n v 
% v = an + b n^perp -> -an + b n^perp
% n^perp = [-ny nx]
% = [ -a nx - b ny, -a ny + b nx ]
/reflection { 8 dict begin
normalized
aload pop
/vy exch def
/vx exch def
normalized
aload pop
/ny exch def
/nx exch def
/a vx nx mul vy ny mul add def
/b vy nx mul vx ny mul sub def
[
  a neg nx mul b ny mul sub
  b nx mul a ny mul sub
]
end } def

% n v index 
% v = a n + b n^perp -> cos(r) n + sin(r) n^perp
% n^perp = [-ny nx]
% = [ -a nx - b ny, -a ny + b nx ]
/refraction { 8 dict begin
/index exch def
normalized
aload pop
/vy exch def
/vx exch def
normalized
aload pop
/ny exch def
/nx exch def
% a = v.n, b = v.n^perp
/a vx nx mul vy ny mul add def
/b vy nx mul vx ny mul sub def
/b' b index div def
b' abs 1 gt {
  v n reflection
}{
  a 0 lt {
  	/a' 1 b' dup mul sub sqrt neg def
  }{
  	/a' 1 b' dup mul sub sqrt def
  } ifelse
  [
    a' nx mul b' ny mul sub
    a' ny mul b' nx mul add
  ]
} ifelse
end } def

% u v
/vadd {
aload pop exch	% u vy vx
3 2 roll aload pop	% vy vx ux uy
3 1 roll			% vy uy vx ux
add 				% vy uy ux+vx
3 1 roll add		% ux+vx uy+vy
[ 3 2 roll ] 
} def

% v c
/vscale { 
[ 3 1 roll dup % [ v c c
  3 2 roll % [ c c v
  aload pop % [ c c x y
  4 1 roll	% [ y c c x
  mul 	% [ y c cx
  3 1 roll
  mul ]
} def

% P=[x y] [vx vy] [cx cy] R -> [x0 y0] where P + tv hits circle cx, cy, R 
/hit { 16 dict begin
/R exch def
aload pop
/cy exch def
/cx exch def
aload pop
/vy exch def
/vx exch def
aload pop
/y exch def
/x exch def
/dx x cx sub def
/dy y cy sub def
/A 1 def
/B vx dx mul vy dy mul add 2 mul def
/C dx dup mul dy dup mul add R dup mul sub def
/r [A B C] roots def
r length 0 eq {
  []
}{
r length 1 eq {
  % one hit
  /t r 0 get def
  t 0.001 gt {
    [ x vx t mul add y vy t mul add ]
  }{
    []
  } ifelse
}{
% two hits
  /t r 0 get def
  t 0.001 gt {
    [ x vx t mul add y vy t mul add ]
  }{
    /t r 1 get def
    t 0.001 gt {
      [ x vx t mul add y vy t mul add ]
    }{
      []
    } ifelse
  } ifelse
} ifelse
} ifelse
end } def

% [cx cy] [x y]
% (x-cx)^2+(y-cy)^2 -> [2(x-cx), 2(y-cy)]
/gradient { 8 dict begin
aload pop 
/y exch def
/x exch def
aload pop
/cy exch def
/cx exch def
[ 
  x cx sub 2 mul 
  y cy sub 2 mul 
]

end } def

% n
/perp {
[ exch aload pop % [ x y
neg exch ]
} def

% [x y] -> unit vector
/normalized { 8 dict begin
aload pop 
/y exch def
/x exch def
/r x dup mul y dup mul add sqrt def
[ x r div y r div ]
end } def

/vlength { 4 dict begin
aload pop
dup mul exch dup mul add sqrt

} def

% [x y] [a b c]
/evaluate { 8 dict begin
aload pop
/c exch def
/b exch def
/a exch def
aload pop
/y exch def
/x exch def
a x mul b y mul add c add
end } def

% p v [A B C]
/line-intersection { 8 dict begin
aload pop
/C exch def
/B exch def
/A exch def
aload pop 
/vy exch def
/vx exch def
aload pop
/py exch def
/px exch def
/t 
  [px py][A B C] evaluate
  neg
  A vx mul B vy mul add 
  div
def
/x px t vx mul add def
/y py t vy mul add def
[x y] 
end } def









gsave

44 20 translate



/f { 4 dict begin
/y exch def
/i y asin def
/r y 1.33 div asin def
4 r mul 2 i mul sub
end } def

100 dup scale
1 100 div setlinewidth

gsave
0.7 setgray
newpath
-0.125 0 moveto
1.25 0 rlineto
0 -0.125 moveto
0 1.375 rlineto
1 -0.125 moveto
0 1.375 rlineto
stroke

newpath
0 10 50 { /i exch def
-0.125 i 90 div 2 mul moveto
1.25 0 rlineto
} for
stroke

grestore

/N 100 def
/y 0 def
/dy 0.99 N div def

newpath
1 0 0 setrgbcolor
currentlinewidth 2 mul setlinewidth
y y f 90 div 2 mul moveto
N {
/y y dy add def
y y f 90 div 2 mul lineto
} repeat
stroke

grestore

/Helvetica-Bold findfont
11 scalefont
setfont

gsave

14 0 translate

30 20 moveto
-8 2 rlineto
(0) show

/Y 22.1 def

0 Y translate
30 20 moveto
-14.5 2 rlineto
(10) show

0 Y translate
30 20 moveto
-14.5 2 rlineto
(20) show

0 Y translate
30 20 moveto
-14.5 2 rlineto
(30) show

0 Y translate
30 20 moveto
-14.5 2 rlineto
(40) show

0 Y translate
30 20 moveto
-14.5 2 rlineto
(50) show

grestore

64 8 moveto
(deflection) show

20 58 translate
90 rotate 
0 0 moveto
(incidence) show

    showpage