%!
%(lenses.inc) run









%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










72 72 scale
0.01 setlinewidth
4 6 translate

newpath
0 0 1 0 360 arc
stroke

newpath 
2 2 moveto
2 -2 lineto
stroke 


0.8 0.9 1 setrgbcolor
0 0 1 0 360 arc
fill 
stroke


/c -0.8 def


34{
/P1 [-4 c] [1 0] [0 0] 1 hit def
/x0 P1 0 get def
/y0 P1 1 get def
0 1 0 setrgbcolor
 
newpath
-4 c moveto
x0 y0 lineto
stroke

/n1 [0 0] P1 gradient normalized def
/v1 n1 [1 0] red-index refraction def
/P2 [x0 y0]v1[0 0] 1 hit def
/x1 P2 0 get def
/y1 P2 1 get def

newpath
x0 y0 moveto
x1 y1 lineto
stroke


/n2 [0 0] P2 gradient normalized def
/v2 n2 v1 1 red-index div refraction def
/P3 [x1 y1]v2[1 0 -4.5] line-intersection def
/x2 P3 0 get def
/y2 P3 1 get def

newpath
x1 y1 moveto
x2 y2 lineto
stroke


/c c 0.05 add def
}repeat








showpage