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

  % the zoom lens with big focal length
8 dup scale
1 7 div setlinewidth
5 10 translate
	%%red 
0 0 0.3 0 360 arc
1 0 0 setrgbcolor
%stroke
0 setgray	
	%lens 1 **********
	%% R1 @@@@@@@@@@
newpath
-13.5 35  14.3416 -20 20 arc
stroke
	%%% R2 @@@@@@@@@
newpath
-18.3 35 24.109 -12 12 arc
stroke
	newpath
0 30 moveto 
5.15700 30 lineto
stroke
	newpath
0 40 moveto 
5.15700 40 lineto
stroke
	%lens 2 *********
	% R3 @@@@@@@@@
newpath
32 35 26.025 167 193 arc
stroke
	% R4 @@@@@@@@@
newpath
-3 35 18.5430 -19 19 arc
stroke
0 0 1 setrgbcolor % blue
	
0 setgray
	newpath
6.5 41 moveto 
14.5 41 lineto
stroke
	newpath
6.5 29 moveto 
14.5 29 lineto
stroke
	%lens 3 ************
	% R5 @@@@@@@@@@@
/x 8 -5 add def
newpath
1 x add 35 27.285 -26 26 arc
stroke
	% R6 @@@@@@@@@
newpath
274.5 x add 35 245.74 177 183 arc
stroke

0 setgray
	newpath
25.5 x add 48 moveto 
29 x add 48 lineto
stroke
	newpath
25.5 x add 22 moveto 
29 x add 22 lineto
stroke
	newpath %vertical lines
25.5 x add 48 moveto 
25.5 x add 47 lineto
stroke
	newpath
25.5 x add 22 moveto 
25.5 x add 23 lineto
stroke

	% trace 1 *****************
	% distance between
newpath
15.5 35 moveto 
15.5 15.79 add 35 lineto
1 0 0 setrgbcolor
0.5 setlinewidth
%stroke

/N 10 def
/y 31 def
/dy 9 N div def

N {
/P0 [-10 y] def
/v0 [1 0] def

/P1 P0 v0 [-13.5 35] 14.3416 hit def
/n2 [-13.5 35] P1 gradient normalized def
/v1 n2 v0 1.5 refraction normalized def
/P2 P1 v1 [-18.3 35] 24.109 hit  def

/n3 [-18.3 35] P2 gradient normalized def
/v2 n3 v1 1 1.5 div refraction normalized def
/P3 P2 v2 [32 35] 26.025 hit def

/n4 [32 35] P3 gradient normalized def
/v3 n4 v2 1.5 refraction normalized def
/P4 P3 v3 [-3 35] 18.543 hit def

/n5 [-3 35] P4 gradient normalized def
/v4 n5 v3 1 1.5 div refraction normalized def
/P5 P4 v4 [1 x add 35] 27.285 hit def

/n6 [1 x add 35] P5 gradient normalized def
/v5 n6 v4 1.5 refraction normalized def
/P6 P5 v5 [274.5 x add 35] 245.74 hit def

/n7 [274.5 x add 35] P6 gradient normalized def
/v6 n7 v5 1 1.5 div refraction normalized def
/P7 P6 v6 [1 0 -200] line-intersection def

newpath
P0 aload pop moveto
P1 aload pop lineto
P2 aload pop lineto
P3 aload pop lineto
P4 aload pop lineto
P5 aload pop lineto
P6 aload pop lineto
P7 aload pop lineto
1 0 0 setrgbcolor
0 setlinewidth
stroke
	/y y dy add def
}repeat

showpage