%!
%%BoundingBox: 0 0 500 500
(lenses.inc) run

% - Inserting ps3d.inc ----------------------

% - This is ps3d.inc ----------------------------------

% - Copyright Bill Casselman
% - Original version 1.0 November, 1998
% - Version 1.1 December, 1999
% -      Took out aliases for moveto etc.
% -      Made display-matrix a 3 x 4 homogeneous matrix, added it to the 3d gstack 
% -      Allowed arbitrary eye location, 3d -> 2d projects onto z = 0 plane
% -      Although fancy stuff not yet implemented
% -      Made ght a variable
% -      Version 1.1 is *not* backwards compatible!
% - Version 1.2 August, 2002
% -      Thorough interpretation of matrices as arrays of columns and point vectors as rows
% -      And some speed up
% -      Similar rewriting of matrix.inc
% -      Again, backwards incompatible!
% - Thanks to Jim Blinn's book `A trip down the graphics pipeline'
% - for several suggestions that (I hope) made this code cleaner 
% - by suggesting how useful homogeneous coordinates were.

% ------------------------------------------------


% - Inserting matrix.inc ----------------------

% - Vector calculations (usually good in any number of dimensions) ----------

% - matrices in this package are laid out in a single array by columns ---------

% a double array: cannot be empty - just lays out all items
/to-single-array {
% [ [. . . ][ . . . ] ] => [ . . . . . . ]
[
	exch { % successive rows
		aload pop
	} forall
]
} def

% ----------------------------------------------

% [ ... ] a square matrix made into an array of columns
/to-double-array { 4 dict begin
/A exch def
/N A length sqrt round cvi def
/i 0 def
[
N {
	[
	N {
		A i get
		/i i 1 add def
	} repeat
	]
} repeat
]
end } def

% ----------------------------------------

% returns the nxn identity matrix as single array
/identity { 1 dict begin
/n exch def
[
	n 1 sub {
		1 
		n {
			0
		} repeat 
	} repeat
	1
]
end } def

% --- vector algebra --------------------------------

% u v -> u.v
/dot-product { 1 dict begin
/v exch def
0 0	     		% u s i
3 2 roll { 		% s i u[i]
  v 			% s i u[i] v
  2 index get mul 	% s i u[i]*v[i]
  3 2 roll		% i u[i]*v[i] s
  add exch 1 add	% s i
} forall
pop
end } def

% v c -> c.v
/vector-scale { 1 dict begin
/c exch def
[ exch
{ 		% s i u[i]
  c mul			% s i u[i] v 
} forall
]
end } def

% u v -> u+v
/vector-add { 1 dict begin
/v exch def
[ exch
0 	     	% u i
exch { 		% i u[i]
  v 		% i u[i] v
  2 index get add 	% i u[i]+v[i]
  exch 1 add	% i
} forall
pop
]
end } def

% u v -> u-v
/vector-sub { 1 dict begin
/v exch def
[ exch
0 	     	% u i
exch { 		% i u[i]
  v 		% i u[i] v
  2 index get sub 	% i u[i]+v[i]
  exch 1 add	% i
} forall
pop
]
end } def

% [x y z ... ] -> r
% watch out for overflow

/vector-length { 1 dict begin
dup
% find maximum entry
/max 0 def
{ % max 
  abs dup max gt {
    % if abs gt max
    /max exch def
  } {
    pop
  } ifelse
} forall

max 0 ne {
  0 exch 
  {  % 0 v[i]
    max div dup mul add
  } forall
  sqrt
  max mul
} {
  pop 0
} ifelse
end } def

% v -> v/|v|
/normalized { 1 dict begin
dup 			% v v
vector-length /r exch def
[ exch
  {
    r div
  } forall
]
end } def

% u v
% u0 u1 u2
% v0 v1 v2
% -> u x v
/cross-product {
2 dict begin
/v exch def
/u exch def
[
  u 1 get v 2 get mul
  v 1 get u 2 get mul sub

  v 0 get u 2 get mul
  u 0 get v 2 get mul sub

  u 0 get v 1 get mul
  v 0 get u 1 get mul sub
]
end
} def

% --------------------------------------------------------------

% axis A -> a matrix
/rotation-matrix3d {
8 dict begin
dup 
cos /c exch def
sin /s exch def
/a exch def
/r a vector-length def
/a0 a 0 get r div def
/a1	a 1 get r div def
/a2 a 2 get r div def

[

% e = [1 0 0] etc.
% e0 = (e.a)a, e# = e - e0, e* = a x e = a x e0 + a x e# = a x e#
/x a0 def
/e0 [a0 x mul
     a1 x mul 
     a2 x mul] def
/e# [1 e0 0 get sub 
     e0 1 get neg 
     e0 2 get neg] def
% [a0 a1 a2]
% [ 1  0  0]
/e* [0 
     a2 
     a1 neg
    ] def
e# 0 get c mul e* 0 get s mul add e0 0 get add 
e# 1 get c mul e* 1 get s mul add e0 1 get add 
e# 2 get c mul e* 2 get s mul add e0 2 get add 

/x a1 def
/e0 [a0 x mul a1 x mul a2 x mul] def
/e# [e0 0 get neg 
     1 e0 1 get sub 
     e0 2 get neg] def
% [a0 a1 a2]
% [ 0  1  0]
/e* [a2 neg 
     0
     a0
    ] def
e# 0 get c mul e* 0 get s mul add e0 0 get add 
e# 1 get c mul e* 1 get s mul add e0 1 get add 
e# 2 get c mul e* 2 get s mul add e0 2 get add 

/x a2 def
/e0 [a0 x mul a1 x mul a2 x mul] def
/e# [e0 0 get neg
     e0 1 get neg
     1 e0 2 get sub] def
% [a0 a1 a2]
% [ 0  0  1]
/e* [a1 
     a0 neg
     0
    ] def
e# 0 get c mul e* 0 get s mul add e0 0 get add 
e# 1 get c mul e* 1 get s mul add e0 1 get add 
e# 2 get c mul e* 2 get s mul add e0 2 get add 

]
% [ r0 r1 r2 r3 r4 r5 r6 r7 r8 ] -> [r0 r3 r6 r1 r4 r7 r2 r5 r8 ]
/r exch def
[
	r 0 get 
	r 3 get
	r 6 get
	r 1 get
	r 4 get
	r 7 get
	r 2 get
	r 5 get
	r 8 get
] 
end
} def

% a v -> v - 2(a.v)/(a.a) a
/euclidean-reflect { 16 dict begin
/v exch def
/a exch def
/N a length def
/d a v dot-product a dup dot-product div 2 mul def
[
  0
  v {
    % i v[i]
    exch dup 1 add					% v[i] i i+1
	3 1 roll							  % i+1 v[i] i
    a exch get d mul 				% i+1 v[i] a[i]*d
	sub										% i+1 v[i]-d*a[i]
	exch									% rv[i] i+1
  } forall
  pop
]
end } def

% f = [A B C] => linear 3d transformation: f not necessarily normalized
% Rv = v - 2 <f', v> f'

/reflection-matrix3d {
3 dict begin
aload pop
/C exch def
/B exch def
/A exch def
/r [ A B C ] vector-length def
/A A r div def
/B B r div def
/C C r div def
[
  1 A A mul dup add sub   B A mul dup add neg   C A mul dup add neg
    A B mul dup add neg 1 B B mul dup add sub   C B mul dup add neg
    A C mul dup add neg   B C mul dup add neg 1 C C mul dup add sub
]
end
} def

/matrix-mul {
8 dict begin
/B exch def
/A exch def
/n A length sqrt round cvi def
%  0    1  ...
%  n  n+1 ...
% 2n 2n+1 ...
[
  	0
  	n {
    	% i = initial index of the column on the stack = 0, n, 2n ... after all previous entries
    	dup n add exch 0	
    	n {
			% i+n i j on stack
      		% j = initial index of the row
			2 copy 1 add % i+n i j i j+1
			4 2 roll	% i+n i j+1 i j
			/ell exch def
      		/k exch def
			% i+n i j+1 
      		0
      		n {
				% i+n i j+1 s
        		A ell get B k get mul add
        		/k k 1 add def
        		/ell ell n add def
      		} repeat
			4 1 roll
			% s i+n i j+1
   	 	} repeat
    	pop pop	% s i+n
  	} repeat
	pop
]
end
} def

% A v: A = [ column 1, column 2, ... , column n ]
/matrix-vector {
8 dict begin
/v exch def
/r v length def
/A exch def
/c A length r idiv def
[
  0 1 c 1 sub {
	/i exch def
    % i = initial index of the row
    0 0
    r {
		% s j on stack
		dup 1 add	% s j j+1
		3 1 roll		% j+1 s j
      	v exch get A i get mul add	% j+1 s
		exch
   	  	/i i r add def
    } repeat
	% s r
	pop
  } for
]
end
} def

% v A: A = [ column1 column2 ... ]
/vector-matrix {
8 dict begin
/A exch def
/v exch def
/c v length def
/r A length c idiv def
[
  /i 0 def 
  r {
    % i = initial index of the row
    /j 0 def
    0
    c {
      A i get v j get mul add
      /j j 1 add def
   	  /i i 1 add def
    } repeat
  } repeat
]
end
} def

% a square matrix m x m
% [i, j] = n*i + j

/transpose {
4 dict begin
/M exch def
/n M length sqrt round cvi def
[
/k 0 def
n {
  /i k def
  n { 
    M i get
    /i i n add def
  } repeat
  /k k 1 add def
} repeat
]

end
} def

/3x3-det {
1 dict begin
/m exch def

m 0 get 
m 4 get mul
m 8 get mul

m 1 get 
m 5 get mul
m 6 get mul
add

m 2 get 
m 3 get mul
m 7 get mul
add

m 2 get 
m 4 get mul
m 6 get mul
sub

m 1 get 
m 3 get mul
m 8 get mul
sub

m 0 get 
m 5 get mul
m 7 get mul
sub

end
} def

/3x3-inverse {
2 dict begin
/m exch def
/d m 3x3-det def
[
   m 4 get m 8 get mul 
   m 5 get m 7 get mul sub
   d div
   
   m 2 get m 7 get mul 
   m 1 get m 8 get mul sub
   d div
   
   m 1 get m 5 get mul 
   m 4 get m 2 get mul sub
   d div

   m 5 get m 6 get mul 
   m 3 get m 8 get mul sub
   d div
   
   m 0 get m 8 get mul 
   m 2 get m 6 get mul sub
   d div
   
   m 2 get m 3 get mul 
   m 0 get m 5 get mul sub
   d div

   m 3 get m 7 get mul 
   m 6 get m 4 get mul sub
   d div
   
   m 1 get m 6 get mul 
   m 0 get m 7 get mul sub
   d div
   
   m 0 get m 4 get mul 
   m 1 get m 3 get mul sub
   d div
]
end
} def

/acos {
dup dup % x x x
mul 1 sub neg % x 1-x^2
sqrt exch atan
} def

% u v

/angle-between {
dup vector-length
% u v |v|
3 1 roll
% |v| u v
1 index 
% |v| u v u
dot-product
% |v| u u.v
exch vector-length
% |v| u.v |u|
div exch div
acos
} def

% x a av -> x - <a, x> av

/skew-reflect { 4 dict begin
/av exch def
/a exch def
/x exch def
/d x a dot-product def
[
	0 1 x length 1 sub {
		/i exch def
		x i get av i get d mul sub
	} for
]
end } def

% a av -> matrix
/skew-reflection-matrix { 8 dict begin
/av exch def
/a exch def
/n a length def
[
	0 1 n 1 sub {
		/i exch def
		[ n {0} repeat ] dup i 1 put
		% e[i] 
		a av skew-reflect
	} for
]
to-single-array
transpose
end } def

% - closing matrix.inc ------------------------

% - Defining PostScript commands' equivalents ----------

% - Coordinates in three dimensions -------------------

% There are three steps to drawing something in 3D:
% 1.  Converting from user 3d coords to default 3d coords
% 2.  Projecting onto (x, y)-plane
% 3.  Drawing the image on that plane
% These are more or less independent.  
% The last step is completely controlled by the usual PS stuff.
% The first and second are handled by 3d stuff here.

% - Initialize and manipulate the 3d gstack -----------

/gmax 64 def
/ght 0 def

% gstack = [[t0 d0 dm0] [t1 d2 dm1] [t2 d2 dm2] ... [t(gmax-1) d(gmax-1) dm(gmax-1)] ]

% the ctm is t[ght]
% the dual ctm is at d[ght]
% they are 4 x 4
% display-matrix = 3 x 4

/gstack3d gmax array def

% start with orthogonal projection to positive z-axis
gstack3d 0 [
  [ 1 0 0 0 
    0 1 0 0
    0 0 1 0
    0 0 0 1] 
  [ 1 0 0 0 
    0 1 0 0
    0 0 1 0
    0 0 0 1]
  [ -1 0 0 0
    0 -1 0 0
    0 0 0 -1]
] put

/gsave3d {
  /ctm gstack3d ght get def
  /ght ght 1 add def 
  ght gmax eq {
    (3d graphics stack overflow!) ==
    quit
  } if
  gstack3d 
    ght
    [ ctm 0 get ctm 1 get ctm 2 get ] 
  put
} def

/grestore3d {
  /ght ght 1 sub def
  ght 0 lt {
     (3d graphics stack underflow!) ==
     quit
  } if
} def

% n - restores to depth n

/gpop3d {
/ght exch def
} def

% [T T* dm]: sets ctm3d = [T T* dm]

/gset3d {
gstack3d  		% [T T* d] g 
ght			    % [T T* d] g ght
3 2 roll 		    % g ght [T T* d]
put 
} def

% => [T T* dm]

/ctm3d {
  gstack3d ght get
} def

% cpt3d isthe last 3d point drawn to

/currentpoint3d {
cpt3d
ctm3d 1 get transform3d
aload pop     % x y z w 
pop			
} def

% - Sets up projection into 2D -------------------
		
% O = [x y z w] 
% sets display-matrix perspective onto z = 0 plane

/display-matrix {
  ctm3d 2 get
} def

% [-z0 0 x0 0
%  0 -z0 y0 0
%  0 0 w0 -z0]
% (transposed)
% gives perspective onto point in z=0 plane
% from [x0 y0 z0 w0]

/set-eye {
4 dict begin
aload pop
/w0 exch def
neg /z0 exch def
/y0 exch def
/x0 exch def
gstack3d ght get
2 
[ z0 0 x0 0 
  0 z0 y0 0
  0 0 w0 z0]
put
end
} def

/get-eye { 1 dict begin
/d display-matrix def
[d 2 get d 6 get d 0 get neg d 10 get]
end } def

% - backwards compatibility -----------------------------

/origin { get-eye } def
/eye { get-eye } def
/set-display { set-eye } def

% - Manipulate the current transformation matrix -----

% x y z

/translate3d {
8 dict begin
/z exch def
/y exch def
/x exch def
/ctm ctm3d def
/T ctm 0 get def
[
[
  T 0 get 
  T 1 get 
  T 2 get 
  T 0 get x  mul
  T 1 get y  mul
  add
  T 2 get z  mul
  add
  T 3 get 
  add
  
  T 4 get 
  T 5 get 
  T 6 get
  T 4 get x  mul
  T 5 get y  mul
  add
  T 6 get z  mul
  add
  T 7 get 
  add
   
  T 8 get 
  T 9 get 
  T 10 get
  T 8 get x  mul
  T 9 get y  mul
  add
  T 10 get z  mul
  add
  T 11 get 
  add
   
  T 12 get 
  T 13 get 
  T 14 get
  T 12 get x  mul
  T 13 get y  mul
  add
  T 14 get z  mul
  add
  T 15 get 
  add
]
/T ctm 1 get def
[
  T 0 get T 12 get x mul sub 
  T 1 get T 13 get x mul sub 
  T 2 get T 14 get x mul sub 
  T 3 get T 15 get x mul sub 

  T 4 get T 12 get y mul sub 
  T 5 get T 13 get y mul sub 
  T 6 get T 14 get y mul sub 
  T 7 get T 15 get y mul sub 

  T 8 get T 12 get z mul sub 
  T 9 get T 13 get z mul sub 
  T 10 get T 14 get z mul sub 
  T 11 get T 15 get z mul sub 

  T 12 get
  T 13 get
  T 14 get
  T 15 get
]
  ctm 2 get
]
end
gset3d
} def

% ------------------------------------------------------

% axis A

/rotate3d { 4 dict begin
rotation-matrix3d
/R exch def
/C ctm3d def
/T C 0 get def
[
[
  % first row
  T 0 get R 0 get mul
  T 1 get R 3 get mul
  add
  T 2 get R 6 get mul
  add
  
  T 0 get R 1 get mul
  T 1 get R 4 get mul
  add
  T 2 get R 7 get mul
  add
  
  T 0 get R 2 get mul
  T 1 get R 5 get mul
  add
  T 2 get R 8 get mul
  add
  
  T 3 get
  
  % second row
  T 4 get R 0 get mul
  T 5 get R 3 get mul
  add
  T 6 get R 6 get mul
  add
  
  T 4 get R 1 get mul
  T 5 get R 4 get mul
  add
  T 6 get R 7 get mul
  add
  
  T 4 get R 2 get mul
  T 5 get R 5 get mul
  add
  T 6 get R 8 get mul
  add
  
  T 7 get
  
  % third row
  T 8 get R 0 get mul
  T 9 get R 3 get mul
  add
  T 10 get R 6 get mul
  add
  
  T 8 get R 1 get mul
  T 9 get R 4 get mul
  add
  T 10 get R 7 get mul
  add
  
  T 8 get R 2 get mul
  T 9 get R 5 get mul
  add
  T 10 get R 8 get mul
  add
  
  T 11 get
  
  % fourth row
  T 12 get R 0 get mul
  T 13 get R 3 get mul
  add
  T 14 get R 6 get mul
  add
  
  T 12 get R 1 get mul
  T 13 get R 4 get mul
  add
  T 14 get R 7 get mul
  add
  
  T 12 get R 2 get mul
  T 13 get R 5 get mul
  add
  T 14 get R 8 get mul
  add
  
  T 15 get
]
/T C 1 get def
% T = T^-1
% => R^-1 T^-1 
[
  R 0 get T 0 get mul
  R 3 get T 4 get mul add
  R 6 get T 8 get mul add
  
  R 0 get T 1 get mul
  R 3 get T 5 get mul add
  R 6 get T 9 get mul add
  
  R 0 get T 2 get mul
  R 3 get T 6 get mul add
  R 6 get T 10 get mul add
  
  R 0 get T 3 get mul
  R 3 get T 7 get mul add
  R 6 get T 11 get mul add
  
  % ------------------------
  
  R 1 get T 0 get mul
  R 4 get T 4 get mul add
  R 7 get T 8 get mul add
  
  R 1 get T 1 get mul
  R 4 get T 5 get mul add
  R 7 get T 9 get mul add
  
  R 1 get T 2 get mul
  R 4 get T 6 get mul add
  R 7 get T 10 get mul add
  
  R 1 get T 3 get mul
  R 4 get T 7 get mul add
  R 7 get T 11 get mul add
  
  % ------------------------
  
  R 2 get T 0 get mul
  R 5 get T 4 get mul add
  R 8 get T 8 get mul add
  
  R 2 get T 1 get mul
  R 5 get T 5 get mul add
  R 8 get T 9 get mul add
  
  R 2 get T 2 get mul
  R 5 get T 6 get mul add
  R 8 get T 10 get mul add
  
  R 2 get T 3 get mul
  R 5 get T 7 get mul add
  R 8 get T 11 get mul add

  T 12 get
  T 13 get
  T 14 get
  T 15 get  
]
  C 2 get
]
end
gset3d
} def


% f = [A B C D] P
% f = 0 is the *affine* reflection plane
% v = v* + v0 with v* on f = 0 and v0 in P-direction => v* - v0
% The map is Q => f(P)*Q - 2*f(Q)P
% It is of order two.
%
% f(P) I -
%
% 2A*P[0] 2B*P[0] 2C*P[0] 2D*P[0]
% 2A*P[1] 2B*P[1] 2C*P[1] 2D*P[1]
% 2A*P[2] 2B*P[2] 2C*P[2] 2D*P[2]
% 2A*P[3] 2B*P[3] 2C*P[3] 2D*P[3]
%
% Matrix = f(P) I - P f
% set s0 = (T row 0)*P
% T x this =
% f(P)T[0,0]-A*s -B*s -C*s -D*s

/affine-reflect3d { 4 dict begin
aload pop
/P3 exch def
/P2 exch def
/P1 exch def
/P0 exch def
aload pop
/D exch def
/C exch def
/B exch def
/A exch def

/fP 
  A P0 mul 
  B P1 mul add
  C P2 mul add
  D P3 mul add
def
/A A dup add def
/B B dup add def
/C C dup add def
/D D dup add def
[
/T ctm3d 0 get def
[
  /s % = (T row 1)*P
  T 0 get P0 mul
  T 1 get P1 mul add
  T 2 get P2 mul add
  T 3 get P3 mul add
  def
  
  fP T 0 get mul
    A s mul sub
  fP T 1 get mul
    B s mul sub
  fP T 2 get mul
    C s mul sub
  fP T 3 get mul
    D s mul sub

  /s % = (T row 1)*P
  T 4 get P0 mul
  T 5 get P1 mul add
  T 6 get P2 mul add
  T 7 get P3 mul add
  def
  
  fP T 4 get mul
    A s mul sub
  fP T 5 get mul
    B s mul sub
  fP T 6 get mul
    C s mul sub
  fP T 7 get mul
    D s mul sub

  /s % = (T row 2)*P
  T 8 get P0 mul
  T 9 get P1 mul add
  T 10 get P2 mul add
  T 11 get P3 mul add
  def
  
  fP T 8 get mul
    A s mul sub
  fP T 9 get mul
    B s mul sub
  fP T 10 get mul
    C s mul sub
  fP T 11 get mul
    D s mul sub

  /s % = (T row 3)*P
  T 12 get P0 mul
  T 13 get P1 mul add
  T 14 get P2 mul add
  T 15 get P3 mul add
  def
  
  fP T 12 get mul
    A s mul sub
  fP T 13 get mul
    B s mul sub
  fP T 14 get mul
    C s mul sub
  fP T 15 get mul
    D s mul sub
]
/T ctm3d 1 get def
/f0 		
  % f paired with columns of T
  T 0 get A mul
  T 4 get B mul add
  T 8 get C mul add
  T 12 get D mul add
def

/f1 		
  % f paired with columns of T
  T 1 get A mul
  T 5 get B mul add
  T 9 get C mul add
  T 13 get D mul add
def

/f2 		
  % f paired with columns of T
  T 2 get A mul
  T 6 get B mul add
  T 10 get C mul add
  T 14 get D mul add
def

/f3 		
  % f paired with columns of T
  T 3 get A mul
  T 7 get B mul add
  T 11 get C mul add
  T 15 get D mul add
def

[
  fP T 0 get mul
    f0 P0 get mul sub
  fP T 1 get mul
    f1 P0 get mul sub
  fP T 2 get mul
    f2 P0 get mul sub
  fP T 3 get mul
    f3 P0 get mul sub
  fP T 4 get mul
    f0 P1 get mul sub
  fP T 5 get mul
    f1 P1 get mul sub
  fP T 6 get mul
    f2 P1 get mul sub
  fP T 7 get mul
    f3 P1 get mul sub
  fP T 8 get mul
    f0 P2 get mul sub
  fP T 9 get mul
    f1 P2 get mul sub
  fP T 10 get mul
    f2 P2 get mul sub
  fP T 11 get mul
    f3 P2 get mul sub
  fP T 12 get mul
    f0 P3 get mul sub
  fP T 13 get mul
    f1 P3 get mul sub
  fP T 14 get mul
    f2 P3 get mul sub
  fP T 15 get mul
    f3 P3 get mul sub
]
  ctm3d 2 get
] 
end
gset3d
} def

% 3x3 M 

/concat3d {
4 dict begin
/M exch def
[
/T ctm3d 0 get def
  [
    T 0 get M 0 get mul
    T 1 get M 3 get mul add
    T 2 get M 6 get mul add

    T 0 get M 1 get mul
    T 1 get M 4 get mul add
    T 2 get M 7 get mul add
     
    T 0 get M 2 get mul
    T 1 get M 5 get mul add
    T 2 get M 8 get mul add
    
    T 3 get
     
    T 4 get M 0 get mul
    T 5 get M 3 get mul add
    T 6 get M 6 get mul add

    T 4 get M 1 get mul
    T 5 get M 4 get mul add
    T 6 get M 7 get mul add
     
    T 4 get M 2 get mul
    T 5 get M 5 get mul add
    T 6 get M 8 get mul add
    
    T 7 get
     
    T 8 get M 0 get mul
    T 9 get M 3 get mul add
    T 10 get M 6 get mul add

    T 8 get M 1 get mul
    T 9 get M 4 get mul add
    T 10 get M 7 get mul add
     
    T 8 get M 2 get mul
    T 9 get M 5 get mul add
    T 10 get M 8 get mul add
    
    T 11 get
     
    T 12 get M 0 get mul
    T 13 get M 3 get mul add
    T 14 get M 6 get mul add

    T 12 get M 1 get mul
    T 13 get M 4 get mul add
    T 14 get M 7 get mul add
     
    T 12 get M 2 get mul
    T 13 get M 5 get mul add
    T 14 get M 8 get mul add
    
    T 15 get
     
  ]
/T ctm3d 1 get def
/M M 3x3-inverse def
  [
    M 0 get T 0 get mul
    M 1 get T 4 get mul add
    M 2 get T 8 get mul add
    
    M 0 get T 1 get mul
    M 1 get T 5 get mul add
    M 2 get T 9 get mul add
    
    M 0 get T 2 get mul
    M 1 get T 6 get mul add
    M 2 get T 10 get mul add
    
    M 0 get T 3 get mul
    M 1 get T 7 get mul add
    M 2 get T 11 get mul add
    
    % -----------------------------
    
    M 3 get T 0 get mul
    M 4 get T 4 get mul add
    M 5 get T 8 get mul add
    
    M 3 get T 1 get mul
    M 4 get T 5 get mul add
    M 5 get T 9 get mul add
    
    M 3 get T 2 get mul
    M 4 get T 6 get mul add
    M 5 get T 10 get mul add
    
    M 3 get T 3 get mul
    M 4 get T 7 get mul add
    M 5 get T 11 get mul add
    
    % -----------------------------
    
    M 6 get T 0 get mul
    M 7 get T 4 get mul add
    M 8 get T 8 get mul add
    
    M 6 get T 1 get mul
    M 7 get T 5 get mul add
    M 8 get T 9 get mul add
    
    M 6 get T 2 get mul
    M 7 get T 6 get mul add
    M 8 get T 10 get mul add
    
    M 6 get T 3 get mul
    M 7 get T 7 get mul add
    M 8 get T 11 get mul add
    
    % -----------------------------
    
    T 12 get
    T 13 get
    T 14 get
    T 15 get
  ]
  ctm3d 2 get
]
end
gset3d
} def

%
% v => v - 2 <v, a> a
%
% Matrix = I - 2 a a
% a

/reflect3d { 4 dict begin
reflection-matrix3d 
concat3d
end } def

% [x y z w] [a00 a01 a02 a03 ... ]
% but the vector is a linear function

/dual-transform3d { 4 dict begin
/v exch def
/T exch def
[
  v 0 get T 0 get mul
  v 1 get T 4 get mul add
  v 2 get T 8 get mul add
  v 3 get T 12 get mul add
  
  v 0 get T 1 get mul
  v 1 get T 5 get mul add
  v 2 get T 9 get mul add
  v 3 get T 13 get mul add
  
  v 0 get T 2 get mul
  v 1 get T 6 get mul add
  v 2 get T 10 get mul add
  v 3 get T 14 get mul add
  
  v 0 get T 3 get mul
  v 1 get T 7 get mul add
  v 2 get T 11 get mul add
  v 3 get T 15 get mul add
]
end } def

% 4d to 3d homogeneous
/project3d { 4 dict begin
/T exch def
/v exch def
[
  T 0 get v 0 get mul
  T 1 get v 1 get mul
  add 
  T 2 get v 2 get mul 
  add
  T 3 get v 3 get mul 
  add

  T 4 get v 0 get mul
  T 5 get v 1 get mul
  add 
  T 6 get v 2 get mul 
  add
  T 7 get v 3 get mul 
  add

  T 8 get v 0 get mul
  T 9 get v 1 get mul
  add 
  T 10 get v 2 get mul 
  add
  T 11 get v 3 get mul 
  add
] 
end } def

% [x y z w] [a00 a01 a02 a03 ... ]

/transform3d {
4 dict begin
/T exch def
/v exch def
[
T 0 get v 0 get mul
T 1 get v 1 get mul
add 
T 2 get v 2 get mul 
add
T 3 get v 3 get mul 
add

T 4 get v 0 get mul
T 5 get v 1 get mul
add 
T 6 get v 2 get mul 
add
T 7 get v 3 get mul 
add

T 8 get v 0 get mul
T 9 get v 1 get mul
add 
T 10 get v 2 get mul 
add
T 11 get v 3 get mul 
add

T 12 get v 0 get mul
T 13 get v 1 get mul
add 
T 14 get v 2 get mul 
add
T 15 get v 3 get mul 
add

] 
end
} def

% sx sy sz

/scale3d {
8 dict begin
/sz exch def
/sy exch def
/sx exch def
/T ctm3d 0 get def
[
[
  T 0 get sx mul
  T 1 get sy mul
  T 2 get sz mul
  T 3 get 
  
  T 4 get sx mul
  T 5 get sy mul
  T 6 get sz mul
  T 7 get 
  
  T 8 get sx mul
  T 9 get sy mul
  T 10 get sz mul
  T 11 get 
 
  T 12 get sx mul
  T 13 get sy mul
  T 14 get sz mul
  T 15 get 
]
/T ctm3d 1 get def
[
  T 0 get sx div
  T 1 get sx div
  T 2 get sx div
  T 3 get sx div
  
  T 4 get sy div
  T 5 get sy div
  T 6 get sy div
  T 7 get sy div
  
  T 8 get sz div
  T 9 get sz div
  T 10 get sz div
  T 11 get sz div
 
  T 12 get 
  T 13 get 
  T 14 get 
  T 15 get 
]
  ctm3d 2 get
]
end
gset3d
} def

% [ <9> ] i

/row {
4 dict begin
/i exch def
/a exch def
a length 9 eq {
  /i i 3 mul def
  /n i 2 add def
} {
  /i i 4 mul def
  /n i 3 add def
} ifelse
[
  i 1 n {
    a exch get
  } for
]
end
} def

% projects from P onto f=Ax+By+Cz+D=0
% two parameters: f = [A B C D] and P
% The map is Q => f(P)*Q - f(Q)P
% It is idempotent.
%
% f(P) I -
%
% A*P[0] A*P[1] A*P[2] A*P[3]
% B*P[0] B*P[1] B*P[2] B*P[3]
% C*P[0] C*P[1] C*P[2] C*P[3]
% D*P[0] D*P[1] D*P[2] D*P[3]
%
% Matrix = f(P) I - P f
% set s0 = (T row 0)*P
% T x this =
% f(P)T[0,0]-A*s -B*s -C*s -D*s

/plane-project {
12 dict begin
aload pop
/P3 exch def
/P2 exch def
/P1 exch def
/P0 exch def
aload pop
/D exch def
/C exch def
/B exch def
/A exch def

/fP 
  A P0 mul 
  B P1 mul add
  C P2 mul add
  D P3 mul add
def
[
/T ctm3d 0 get def
[
  /s % = (T row 0)*P
  	T 0 get P0 mul
  	T 1 get P1 mul add
  	T 2 get P2 mul add
  	T 3 get P3 mul add
  def
  
  fP T 0 get mul
    A s mul sub
  fP T 1 get mul
    B s mul sub
  fP T 2 get mul
    C s mul sub
  fP T 3 get mul
    D s mul sub


  /s % = (T row 1)*P
  T 4 get P0 mul
  T 5 get P1 mul add
  T 6 get P2 mul add
  T 7 get P3 mul add
  def
  
  fP T 4 get mul
    A s mul sub
  fP T 5 get mul
    B s mul sub
  fP T 6 get mul
    C s mul sub
  fP T 7 get mul
    D s mul sub


  /s % = (T row 2)*P
  T 8 get P0 mul
  T 9 get P1 mul add
  T 10 get P2 mul add
  T 11 get P3 mul add
  def
  
  fP T 8 get mul
    A s mul sub
  fP T 9 get mul
    B s mul sub
  fP T 10 get mul
    C s mul sub
  fP T 11 get mul
    D s mul sub

  /s % = (T row 3)*P
  T 12 get P0 mul
  T 13 get P1 mul add
  T 14 get P2 mul add
  T 15 get P3 mul add
  def
  
  fP T 12 get mul
    A s mul sub
  fP T 13 get mul
    B s mul sub
  fP T 14 get mul
    C s mul sub
  fP T 15 get mul
    D s mul sub
]
/T ctm3d 1 get def
[
  /s T 0 get T 1 get add T 2 get add T 3 get add def
  s A mul
  s B mul
  s C mul
  s D mul
  /s T 4 get T 5 get add T 6 get add T 7 get add def
  s A mul
  s B mul
  s C mul
  s D mul
  /s T 8 get T 9 get add T 10 get add T 11 get add def
  s A mul
  s B mul
  s C mul
  s D mul
  /s T 12 get T 13 get add T 14 get add T 15 get add def
  s A mul
  s B mul
  s C mul
  s D mul
]
  ctm3d 2 get
] 
end
gset3d
} def

% - Drawing commands in 3D --------------------------

% P displayed = [x y w z] => x/w y/w 

/render {
8 dict begin
aload pop
/v3 exch def
/v2 exch def
/v1 exch def
/v0 exch def
/T display-matrix def
/x
  T 0 get v0 mul 
  T 1 get v1 mul add
  T 2 get v2 mul add
  T 3 get v3 mul add
def
/y 
  T 4 get v0 mul 
  T 5 get v1 mul add
  T 6 get v2 mul add
  T 7 get v3 mul add
def
/w 
  T 8 get v0 mul 
  T 9 get v1 mul add
  T 10 get v2 mul add
  T 11 get v3 mul add
def
w 0 eq {
  (Perspective: division by zero!) ==
  quit
} if
x w div
y w div
end
} def

/cpt3d 4 array def
/lm3d 4 array def

% cpt3d is a point in the "real" 3d world

% Should we build the current 3d path for reuse?

% x y z
/moveto3d {
1 [ 5 1 roll ] 
ctm3d 0 get 
transform3d 
aload pop cpt3d astore 
render moveto
cpt3d aload pop lm3d astore pop
} def

% x y z

/lineto3d {
1 [ 5 1 roll ] 
ctm3d 0 get 
transform3d 
aload pop cpt3d astore 
render lineto
} def

% x y z

/rmoveto3d {
0 [ 5 1 roll ] ctm3d 0 get
transform3d 			% [dx dy dz 0]
dup 0 get cpt3d 0 get add	% [dx dy dz 0] x+dx
exch 
dup 1 get cpt3d 1 get add	% x+dx [dx dy dz dw] y+dy
exch
dup 2 get cpt3d 2 get add	% x+dx x+dy [dx dy dz dw] z+dz
exch				% x+dx y+dy z+dz [dx dy dz dw]
3 get cpt3d 3 get add		% x+dx x+dy z+dz w+dw
cpt3d astore
render moveto
cpt3d aload pop lm3d astore pop
} def

% x y z

/rlineto3d {
0 [ 5 1 roll ] ctm3d 0 get
transform3d 		% [dx dy dz 0]
dup 0 get cpt3d 0 get add	% [dx dy dz 0] x+dx
exch 
dup 1 get cpt3d 1 get add	% x+dx [dx dy dz dw] y+dy
exch
dup 2 get cpt3d 2 get add	% x+dx x+dy [dx dy dz dw] z+dz
exch				% x+dx y+dy z+dz [dx dy dz dw]
3 get cpt3d 3 get add		% x+dx x+dy z+dz w+dw
cpt3d astore  
render lineto
} def

% x1 y1 z1 x2 y2 z2 x3 y3 z3 

/curveto3d { 16 dict begin
/z3 exch def
/y3 exch def
/x3 exch def
/z2 exch def
/y2 exch def
/x2 exch def
/z1 exch def
/y1 exch def
/x1 exch def
% F(t) 
/P0 cpt3d display-matrix project3d def
/P1 [x1 y1 z1 1] ctm3d 0 get transform3d display-matrix project3d def
/w P0 3 get def
/c P1 3 get w div 1 sub def
/x0 P0 0 get w div def
/x1 P1 0 get w div def
/y0 P0 1 get w div def
/y1 P1 1 get w div def
x1 x0 c mul sub 
y1 y0 c mul sub 
[x3 y3 z3 1] ctm3d 0 get transform3d aload pop cpt3d astore 
display-matrix project3d /P3 exch def
/P2 [x2 y2 z2 1] ctm3d 0 get transform3d display-matrix project3d def
% We are assuming the display-matrix has images on { z = 0 }
/w P3 3 get def
/c P2 3 get w div 1 sub def
/x3 P3 0 get w div def
/x2 P2 0 get w div def
/y3 P3 1 get w div def
/y2 P2 1 get w div def
x2 x3 c mul sub 
y2 y3 c mul sub 
x3 y3
curveto
end } def

% - are the next two used? --------------------------------

% dP dQ Q

/dcurveto3d { 16 dict begin
/z3 exch def
/y3 exch def
/x3 exch def
/dz3 exch def
/dy3 exch def
/dx3 exch def
/dz0 exch def
/dy0 exch def
/dx0 exch def
/P0 cpt3d display-matrix project3d def
/dP0 [dx0 dy0 dz0 0] ctm3d 0 get transform3d display-matrix project3d def
/w P0 3 get def 
% c = 1 - w'/w
/c 1 dP0 3 get w div sub def
  P0 0 get w div c mul dP0 0 get w div add
  P0 1 get w div c mul dP0 1 get w div add
[x3 y3 z3 1] ctm3d 0 get transform3d aload pop cpt3d astore 
display-matrix project3d /P3 exch def
/dP3 [dx3 dy3 dz3 0] ctm3d 0 get transform3d display-matrix project3d def
/w P3 3 get def
/c 1 dP3 3 get w div add def
  P3 0 get w div c mul dP3 0 get w div sub
  P3 1 get w div c mul dP3 1 get w div sub
  P3 0 get w div
  P3 1 get w div
curveto
end
} def

% dP dQ Q

/dcurveto { 8 dict begin
/y3 exch def
/x3 exch def
/dy3 exch def
/dx3 exch def
/dy0 exch def
/dx0 exch def
currentpoint dy0 add exch dx0 add exch
x3 dx3 sub y3 dy3 sub
x3 y3 
curveto
end } def

% - are the above two used? ----------------------------

/closepath3d {
closepath
lm3d aload pop cpt3d astore pop
} def

% - Conversion from 2d to 3D ---------------------------

/simple-convert {
[
  { % x y
    [ 3 1 roll 0 {moveto3d}]
  }
  { % x y
    [ 3 1 roll 0 {lineto3d}]
  }
  { % x1 y1 x2 y2 x3 y3
    [ 7 1 roll
    0 5 1 roll
    0 3 1 roll
    0 {curveto3d} ]
  }
  {
    [ {closepath3d} ]
  }
  pathforall
]
newpath
{
  aload pop exec
} forall
} def

% -----------------------------------------------

% For a simple currentpoint:

/mkpath3dDict 8 dict def

mkpath3dDict begin

/pathcount {
0 
{pop pop 
  pop 1 exit}
{pop pop 
  pop 1 exit}
{pop pop pop pop pop pop 
  pop 1 exit}
{ pop 1 exit}
pathforall
} def

/thereisacurrentpoint {
pathcount
0 gt
{true}
{false}
ifelse
} def

end

% ---------------------------------------------

% stack: [parameters] /f t0 t1 N

/mkpath3d {
ps3ddict begin
mkpath3dDict begin
12 dict begin
/N exch def
/t1 exch def
/t exch def
/f exch cvx def
/pars exch def

/h t1 t sub N div def
/h3 h 0.333333 mul def

pars t f aload pop
/velocity exch def
/position exch def
/p [ position aload pop 1 ] 
ctm3d 0 get transform3d display-matrix project3d def
/v [ velocity aload pop 0 ] 
ctm3d 0 get transform3d display-matrix project3d def

% p v = homogeneous position and velocity
% p/p[3] v/p[3] - (v[3]/p[3])(p/p[3]) = inhomogeneous ones
% = p/w v/w - c*p/w

/w p 3 get def
/c v 3 get w div def
thereisacurrentpoint {
  p 0 get w div p 1 get w div lineto
} {
  p 0 get w div p 1 get w div moveto
} ifelse

N {				% x y = currentpoint
  p 0 get 
    v 0 get 
    p 0 get c mul 
    sub 
    h3 mul 
  add 
  w div
  
  p 1 get 
    v 1 get 
    p 1 get c mul 
    sub 
    h3 mul 
  add
  w div

/t t h add def

pars t f aload pop
/velocity exch def
/position exch def
/p [ position aload pop 1 ] 
ctm3d 0 get transform3d display-matrix project3d def
/v [ velocity aload pop 0 ] 
ctm3d 0 get transform3d display-matrix project3d def
/w p 3 get def
/c v 3 get w div def

  p 0 get 
    v 0 get 
    p 0 get c mul 
    sub 
    h3 mul 
  sub 
  w div
  
  p 1 get 
    v 1 get 
    p 1 get c mul 
    sub 
    h3 mul 
  sub
  w div

p 0 get w div
p 1 get w div
curveto

} repeat
end % local dict
end % mkpath3d dict
end % ps3ddict
} def

% makes polygon out of control points

/mkcontrolpath3d {
mkpath3dDict begin
12 dict begin
/N exch def
/t1 exch def
/t exch def
/f exch cvx def
/pars exch def

/h t1 t sub N div def
/h3 h 0.333333 mul def

pars t f aload pop
/velocity exch def
/position exch def

position aload pop
thereisacurrentpoint {
	lineto3d
} {
	moveto3d
} ifelse

N {				% x y = currentpoint
% currentpoint pixel pop

position 0 get
velocity 0 get		% x dx/dt
h3 mul add
position 1 get
velocity 1 get		% y dy/dt
h3 mul add
position 2 get
velocity 2 get		% z dz/dt
h3 mul add

lineto3d

/t t h add def

pars t f aload pop
/velocity exch def
/position exch def

position 0 get
velocity 0 get 
h3 mul sub
position 1 get
velocity 1 get 
h3 mul sub
position 2 get
velocity 2 get 
h3 mul sub

lineto3d

position 0 get
position 1 get
position 2 get

lineto3d

} repeat
end % local dict
end % mkpath3d dict
} def

% -----------------------------------------------

% makes polygon from endpoints

/mkpolypath3d {
mkpath3dDict begin
12 dict begin
/N exch def
/t1 exch def
/t exch def
/f exch cvx def
/pars exch def

/h t1 t sub N div def
/h3 h 0.333333 mul def

pars t f aload pop
/velocity exch def
/position exch def

position aload pop
thereisacurrentpoint {
	lineto3d
} {
	moveto3d
} ifelse

N {				% x y = currentpoint
% currentpoint pixel pop

/t t h add def

pars t f aload pop
/velocity exch def
/position exch def

position 0 get
position 1 get
position 2 get
lineto3d

} repeat
end % local dict
end % mkpath3d dict
} def

% ---------------------------------------------

% length width

/plainarrow3d {
5 dict begin
/shaftwidth exch def
/arrowlength exch def
/headwidth shaftwidth 3 mul def
/headlength headwidth def
/shaftlength arrowlength shaftwidth 2.5 mul sub def
0 0 0 moveto3d
0 shaftwidth 0.5 mul 0 lineto3d
% shaftlength 0 0 rlineto3d
shaftlength shaftwidth 0.5 mul 0 lineto3d
arrowlength headlength sub headwidth 0.5 mul 0 lineto3d
arrowlength 0 0 lineto3d
arrowlength headlength sub headwidth -0.5 mul 0 lineto3d
shaftlength shaftwidth -0.5 mul 0 lineto3d
0 shaftwidth -0.5 mul 0 lineto3d
0 0 0 lineto3d
end
} def

% length width

/plainarrowtail3d {
5 dict begin
/shaftwidth exch def
/arrowlength exch def
0 0 0 moveto3d
0 shaftwidth 0.5 mul 0 lineto3d
arrowlength 0 0 rlineto3d
0 shaftwidth neg 0 rlineto3d
arrowlength neg 0 0 rlineto3d
0 0 0 lineto3d
end
} def

% A face is an array of two items: [polygon array + normal]
% The normal is a linear function [A B C D] f(x) > 0 = outside

% normal-function eye

/is-outside {
dot-product 
-0.0001 ge
} def

% smin + (1-smin)(1 - (1-x)^2)
% ... x(2 - x)
% -1 => smin
%  1 => 1

% [A B C D] + # in [-1 1]: A etc = control points, A = min, D = max
% 1 = facing towards, -1 = facing away from light
% x -> (x+1)/2 = 0.5(x+1) takes [-1, 1] -> [0, 1]
/shade { 8 dict begin
1 add 0.5 mul /s exch def	% s now in [0 1]
/t 1 s sub def % t = a-s
aload pop 
/D exch def
/C exch def
/B exch def
/A exch def

/A' A t mul B s mul add def
/B' B t mul C s mul add def
/C' C t mul D s mul add def
A' t mul B' s mul add t mul 
B' t mul C' s mul add s mul 
add

end } def

% ----------------------------------------------------

% input: [pars] /fcn s0 s1 t0 t1 ns nt
% the fcn: [pars] s t -> f_{pars}(s, t)
% output: a polygonal surface of faces [ normal-fcn triangle ]

/mksurface {
16 dict begin
/nt exch def
/ns exch def
/t1 exch def
/t0 exch def
/s1 exch def
/s0 exch def
/ds s1 s0 sub ns div def
/dt t1 t0 sub nt div def
/f exch cvx def
/pars exch def
/P [
  /s s0 def
  ns 1 add {
    [
      /t t0 def
      nt 1 add {
        pars s t f
        /t t dt add def
      } repeat
    ]
    /s s ds add def
  } repeat
] def
% P[i][j] = f(s0 + i.ds, t0 + j.dt)
[
  0 1 ns 1 sub {
    /i exch def
    0 1 nt 1 sub {
      /j exch def
      % an array of triangles (i, j, 0) + (i, j, 1) 
      % dividing the rectangles in two
      /P00 P i get j get def
      /P10 P i 1 add get j get def
      /P01 P i get j 1 add get def
      /P11 P i 1 add get j 1 add get def
      % normal
      /Q P10 P00 vector-sub P01 P00 vector-sub cross-product def
      /r Q vector-length def
      r 0 ne {
        [
          [ Q aload pop Q P10 dot-product neg ]
          % array of pointers to three vertices
          [ P00 P10 P01 ]
        ] 
      } if
      /Q P01 P11 vector-sub P10 P11 vector-sub cross-product def
      /r Q vector-length def
      r 0 ne {
        [
          [ Q aload pop Q P01 dot-product neg ]
          % array of pointers to three vertices
          [ P10 P11 P01 ]
        ]
      } if
    } for
  } for
]
end
} def

% an array of vertices 
% traversed according to right hand rule
% output normalized

/normal-function {
2 dict begin
/a exch def
/n 
  a 1 get a 0 get vector-sub 
  a 2 get a 1 get vector-sub 
  cross-product def
/r n 0 get dup mul n 1 get dup mul add n 2 get dup mul add sqrt def
/n [ n 0 get r div n 1 get r div n 2 get r div ] def
[
  n aload pop
  a 0 get n dot-product neg
]
end
} def








% - closing ps3d.inc ------------------------


gsave 
72 72 scale 
0.02 setlinewidth 
1 2 translate 



newpath
0.2 setgray
0 -1 moveto 
5 -1 lineto 
5 5 lineto 
0 5 lineto 
closepath
fill


grestore




% -----------------------


gsave
gsave3d
/P {
2 dict begin
/B exch def
/A exch def
% lat long
[A cos B cos mul A sin B cos mul B sin]
end
} def

/page-begin {
gsave
72 dup scale
1 72 div setlinewidth
4 5 translate
1 setlinecap
1 setlinejoin
} def

/page-end {
grestore
showpage
} def

% latitudes in rows 
/sphere-vertex { 1 dict begin
/N exch def
/dA 179.8 N div def
/dB 360 N div 2 div def
/A -89.9 def
% A = latitude
N 1 add { 
  [
  /B 0 def
  % B = longitude
  2 N mul 1 add {
    B A P 
    /B B dB add def
  } repeat
  ]
  /A A dA add def
} repeat
end } def

/N 16 def
/S [ N sphere-vertex ] def

/sphere [
0 1 N 1 sub {
  /i exch def
  0 1 2 N mul 1 sub { /j exch def
    [
      [ S i get j get 
      S i get j 1 add get 
      S i 1 add get j 1 add get ] dup normal-function ] 
    [
      [ S i get j get 
      S i 1 add get j 1 add get 
      S i 1 add get j get ] dup normal-function ] 
  } for
} for
] def

% ------------------------------------------------------------

[0 0 6 1] set-eye 
/light [-1 2 1 0 ] normalized def
/Sh [ 0 0.2 0.8 0.9 ] def

/A 0 def
/dA 10 def
/E get-eye def

 {

page-begin
gsave3d

[0 0 1] 90 rotate3d
[0 1 0] 90 rotate3d
[0 0 1] A rotate3d
/Tinv ctm3d 1 get def
/T ctm3d 0 get def

/O E Tinv transform3d def
/L light Tinv transform3d def

0 1 sphere length 1 sub {
  /i exch def
  /f sphere i get def
  % f = [ point-array visibility ]
  /v f 1 get def
  v O dot-product 0 ge {
     /p f 0 get def

/ds [
	0 [ p 0 get aload pop 1] T transform3d render p 0 get L dot-product Sh exch shade
	0 [ p 1 get aload pop 1] T transform3d render p 1 get L dot-product Sh exch shade
	0 [ p 2 get aload pop 1] T transform3d render p 2 get L dot-product Sh exch shade
] def
gsave
<<
	/ShadingType 4
	/ColorSpace [ /DeviceGray ]
	/DataSource ds
>>
shfill
grestore
 
  } if
} for

grestore3d



0.02 setlinewidth
-4 -5 translate
1 2 translate 


/z 0.35 def
/g 0.6 def
/N 30 def
/x 3.3 def
/dx 3.3 def
/y  0.008 def 


8 {
z setgray


N { 
newpath 
0 x moveto 
/P0 [0 x] def 
/V0 [1 0] def 
/C [3 3] def 
/R 1 def 

/P1 P0 V0 C R hit def 
P1 0 get P1 1 get lineto 

stroke
g setgray
P1 0 get P1 1 get moveto


/n C P1 gradient normalized def 
/V1 n V0 1.33 refraction def 
/P2 P1 V1 C R hit def 
P2 0 get P2 1 get lineto 

/n C P2 gradient normalized def 
/V2 n V1 reflection def 
/P3 P2 V2 C R hit def 
P3 0 get P3 1 get lineto 

stroke
z setgray
P3 0 get P3 1 get moveto

/n C P3 gradient normalized def 
/V3 n V2 1 1.33 div refraction def 
/P4 P3 V3 [1 0 0] line-intersection def 
P4 0 get P4 1 get lineto 
stroke 
/x x y add def 
} repeat 

/z z 0.075 add def
/g g 0.04 add def
/N N 4 sub def
/dx dx 2 y mul add def
/x dx def
} repeat


0.2 setgray
3 3 1.03 0 360 arc
stroke

page-end

/A A dA add def

} loop

grestore3d
grestore

% ------------------------