%!

/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




24 24 scale
0.1 setlinewidth
5 3 translate

newpath
-2.9 0 moveto
-2.9 -1 lineto
6.6 -1 lineto
6.6 0 lineto
-2.9 0 lineto
0 0.3 0.8 setrgbcolor
fill

/P0 [-2.5 2.5] def
/v0 [1 -1] def
/P1 P0 v0 [0 1 0] line-intersection def
/v1 [1 0] v0 1.5 refraction def
/P2 P1 v1 [0 1 1] line-intersection def
/v2 [1 0] v1 reflection def
/P3 P2 v2 [0 1 0] line-intersection def
/v3 [1 0] v2 1 1.5 div refraction def
/P4 P3 v3 [0 1 -2.5] line-intersection def

newpath
P0 0 get P0 1 get moveto
P1 0 get P1 1 get lineto
P2 0 get P2 1 get lineto
P3 0 get P3 1 get lineto
P4 0 get P4 1 get lineto
1 0.2 0.2 setrgbcolor
0.08 setlinewidth
stroke

showpage