%!
%%%%%%%%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

%%%%%%end lenses.inc%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%555
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% distortion %%%%%%%%%%%%%%%%%%%%%

68 dup scale 
3.2 5 translate
1 72 div setlinewidth


/x 2 def
/radius 1 def
newpath
0 x add 0 radius 119 240 arc
-1 x add 0 radius -50 61 arc
gsave
0.8 0.8 1 setrgbcolor
fill
grestore
stroke

newpath
-3 5 moveto
-3 -5 lineto
stroke
newpath
-3 0 moveto
16 0 lineto
stroke
newpath
16 3 div 5 moveto
16 3 div -5 lineto
stroke

newpath
-0.5 0 0.1 0 360 arc
%stroke

newpath
-.5 0.5 moveto
-.5 3 lineto
0 setgray
stroke

newpath
-.5 0.07 moveto
-.5 -3 lineto
stroke

/N 3 def
/vx -0.5 def
/P1 [-3 1] def

N{
1 0 0 setrgbcolor
newpath
/v1 [1 vx 2.4 div] def
/P2 P1 v1 normalized [0 x add 0] radius hit def
/n2 [0 x add 0] P2 gradient normalized def
/v2 n2 v1 1.5 refraction normalized def 
/P3 P2 v2 normalized [-1 x add 0] radius hit def
/n3 [-1 x add 0] P3 gradient normalized def
/v3 n3 v2 1 1.5 div refraction normalized def 
/P4 [16 3 div -0.5] def

P1 aload pop moveto
P2 aload pop lineto
P3 aload pop lineto
P4 aload pop lineto
stroke
/vx vx -0.2 add def
}repeat

showpage