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

% - Copyright Bill Casselman, November 13, 1998 -------
% - Thanks to Jim Blinn's book `A trip down the graphics pipeline'
% - for several suggestions that (I hope) made this code cleaner 


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

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

/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

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

/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

/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

/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 => [ x1 y1 z1 x2 y2 z2 x3 y3 z3 ] = columns of the matrix

/rotation-matrix3d {
8 dict begin
dup 
cos /c exch def
sin /s exch def
/a exch def
/r a vector-length def
/a [ 
  a 0 get r div 
  a 1 get r div 
  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 a 0 get def
/e0 [a 0 get x mul
     a 1 get x mul 
     a 2 get 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 
     a 2 get 
     a 1 get 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 a 1 get def
/e0 [a 0 get x mul a 1 get x mul a 2 get 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* [a 2 get neg 
     0
     a 0 get
    ] 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 a 2 get def
/e0 [a 0 get x mul a 1 get x mul a 2 get 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* [a 1 get 
     a 0 get 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 

]
end
} def

% square matrices

/matrix-mul {
8 dict begin
/B exch def
/A exch def
/N A length sqrt round cvi def
/n N 1 sub def
%  0    1  
%  n  n+1
% 2n 2n+1
[
  /i 0 def 
  N {
    % i = initial index of the row
    /j 0 def
    N {
      % j = initial index of the column
      /k i def
      /ell j def
      0
      N {
        A k get B ell get mul add
        /k k 1 add def
        /ell ell N add def
      } repeat
      /j j 1 add def
    } repeat
    /i i N add def
  } 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




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

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

/ps3ddict 16 dict def

ps3ddict begin

/LINETO {
lineto
} bind def

/MOVETO {
moveto
} bind def

/CLOSEPATH {
closepath
} bind def

/CURVETO {
curveto
} bind def

end

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

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

/gmax 64 def

% gstack = [n [T1 D1] [T2 D2] ... [T(gmax) D(gmax)] ]

% n = gstack3d[0] 
% the ctm is at gstack[n]
% the dual ctm is at gstack[n+1]
/gstack3d gmax 1 add array def
gstack3d 0 1 put
% start with orthogonal projection
gstack3d 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 1 0
	     0 0 0 1]
	    ] put

/gsave3d {
1 dict begin
  /n gstack3d 0 get def
  n gmax eq {
    (3d graphics stack overflow!) ==
    quit
  } if
  gstack3d 
    n 1 add 
    gstack3d n get 
  put
  gstack3d 
    0 
    n 1 add 
  put
end  
} def

/grestore3d {
  gstack3d
    0
    gstack3d 0 get 
    dup 2 lt {
      (3d graphics stack underflow!) ==
      quit
    } if
    1 sub 
  put
} def

% n - restores to depth n

/gpop3d {
gstack3d 0   % n gstack 0
3 2 roll put
} def

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

/gset3d {
gstack3d dup 		% [T T*] g g
0 get 			% [T T*] g n
3 2 roll 		% g n [T T*]
put 
} def

% => [T T*]

/ctm3d {
  gstack3d dup 0 get get
} def

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

/gdepth {
gstack3d 0 get
} def

% - Sets up projection into 2D -------------------

/display-matrix [1 0 0 0 
		0 1 0 0
		0 0 0 0
		0 0 0 1] def
		
% O = [0 0 x 0] (x > 0) or [0 0 a 1] (a > 0)

/set-display {
3 dict begin
/O exch def
/a O 2 get def
/x O 3 get def
[a 0 0 0
  0 a 0 0
  0 0 0 0
  0 0 x neg a]
end
/display-matrix exch def
} def

/origin {
[0 0 display-matrix 15 get display-matrix 14 get neg]
} def

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

% x y z

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

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

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

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

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

% axis A

/rotate3d {
4 dict begin
rotation-matrix3d
/R exch def
/T ctm3d 0 get def
[
[
  % first row
  T 0 get R 0 get mul
  T 1 get R 1 get mul
  add
  T 2 get R 2 get mul
  add
  
  T 0 get R 3 get mul
  T 1 get R 4 get mul
  add
  T 2 get R 5 get mul
  add
  
  T 0 get R 6 get mul
  T 1 get R 7 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 1 get mul
  add
  T 6 get R 2 get mul
  add
  
  T 4 get R 3 get mul
  T 5 get R 4 get mul
  add
  T 6 get R 5 get mul
  add
  
  T 4 get R 6 get mul
  T 5 get R 7 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 1 get mul
  add
  T 10 get R 2 get mul
  add
  
  T 8 get R 3 get mul
  T 9 get R 4 get mul
  add
  T 10 get R 5 get mul
  add
  
  T 8 get R 6 get mul
  T 9 get R 7 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 1 get mul
  add
  T 14 get R 2 get mul
  add
  
  T 12 get R 3 get mul
  T 13 get R 4 get mul
  add
  T 14 get R 5 get mul
  add
  
  T 12 get R 6 get mul
  T 13 get R 7 get mul
  add
  T 14 get R 8 get mul
  add
  
  T 15 get
]
/T ctm3d 1 get def
% T = T^-1
% => R^-1 T^-1 
[
  R 0 get T 0 get mul
  R 1 get T 4 get mul add
  R 2 get T 8 get mul add
  
  R 0 get T 1 get mul
  R 1 get T 5 get mul add
  R 2 get T 9 get mul add
  
  R 0 get T 2 get mul
  R 1 get T 6 get mul add
  R 2 get T 10 get mul add
  
  R 0 get T 3 get mul
  R 1 get T 7 get mul add
  R 2 get T 11 get mul add
  
  % ------------------------
  
  R 3 get T 0 get mul
  R 4 get T 4 get mul add
  R 5 get T 8 get mul add
  
  R 3 get T 1 get mul
  R 4 get T 5 get mul add
  R 5 get T 9 get mul add
  
  R 3 get T 2 get mul
  R 4 get T 6 get mul add
  R 5 get T 10 get mul add
  
  R 3 get T 3 get mul
  R 4 get T 7 get mul add
  R 5 get T 11 get mul add
  
  % ------------------------
  
  R 6 get T 0 get mul
  R 7 get T 4 get mul add
  R 8 get T 8 get mul add
  
  R 6 get T 1 get mul
  R 7 get T 5 get mul add
  R 8 get T 9 get mul add
  
  R 6 get T 2 get mul
  R 7 get T 6 get mul add
  R 8 get T 10 get mul add
  
  R 6 get T 3 get mul
  R 7 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  
]
]
end
gset3d
} def

% [A B C] v

/reflection-transform3d {
5 dict begin
/v exch def
aload pop
/C exch def
/B exch def
/A exch def
/fv 
  v 0 get A mul 
  v 1 get B mul add
  v 2 get C mul add
  dup add
def
[
  v 0 get fv A mul sub
  v 1 get fv B mul sub
  v 2 get fv C mul sub
]
end
} def

% f = [A B C] => linear 3d transformation
% Rv = v - 2 <f, v> f

/reflection-matrix3d {
3 dict begin
aload pop
/C exch def
/B exch def
/A exch 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

% f = [A B C D] u
% f = 0 is the *affine* reflection plane
% v = v* + v0 with v0 in u-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
]
] 
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 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
    
  ]
]
end
gset3d
} def

%
% v => v - 2 <v, a> a
%
% Matrix = I - 2 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
/T exch def
/v 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

% [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 
]
]
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

% f = [A B C D] P
% The map is Q => f(P)*Q - f(Q)P
% It is idempotent.
%
% f(P) I -
%
% A*P[0] B*P[0] C*P[0] D*P[0]
% A*P[1] B*P[1] C*P[1] D*P[1]
% A*P[2] B*P[2] C*P[2] D*P[2]
% A*P[3] B*P[3] C*P[3] 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 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
[
  /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
]
] gset3d
end
} def

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

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

/project {
  render
} def

/render {
8 dict begin
/v exch def
/T display-matrix def
/x
  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
def
/y 
  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
def
/w 
  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
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 {
ps3ddict begin
1 [ 5 1 roll ] 
ctm3d 0 get 
transform3d 
aload pop cpt3d astore 
project MOVETO
end
cpt3d aload pop lm3d astore pop
} def

% x y z

/lineto3d {
ps3ddict begin
1 [ 5 1 roll ] 
ctm3d 0 get 
transform3d 
aload pop cpt3d astore 
project LINETO
end
} def

% x y z

/rmoveto3d {
ps3ddict begin
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
project MOVETO
end
cpt3d aload pop lm3d astore pop
} def

% x y z

/rlineto3d {
ps3ddict begin
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  
project LINETO
end
} def

% x1 y1 z1 x2 y2 z2 x3 y3 z3 

/curveto3d {
ps3ddict begin
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 transform3d def
/P1 [x1 y1 z1 1] ctm3d 0 get transform3d display-matrix transform3d 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 transform3d /P3 exch def
/P2 [x2 y2 z2 1] ctm3d 0 get transform3d display-matrix transform3d 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
end
} def

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

% dP dQ Q

/dcurveto3d {
ps3ddict begin
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 transform3d def
/dP0 [dx0 dy0 dz0 0] ctm3d 0 get transform3d display-matrix transform3d 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 transform3d /P3 exch def
/dP3 [dx3 dy3 dz3 0] ctm3d 0 get transform3d display-matrix transform3d 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
end
} def

% dP dQ Q

/dcurveto {
ps3ddict begin
16 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
end
} def

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

/closepath3d {
ps3ddict begin
CLOSEPATH
end
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 transform3d def
/v [ velocity aload pop 0 ] 
ctm3d 0 get transform3d display-matrix transform3d 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 transform3d def
/v [ velocity aload pop 0 ] 
ctm3d 0 get transform3d display-matrix transform3d 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: [normal polygon array]
% The normal is a linear function [A B C D] f(x) > 0 = outside

% face eye

/visible-to {
exch 0 get 
dot-product 
0 ge
} def

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

/shademin 0.4 def

% [-1 1]

/shade {
1 add 0.5 mul		% [0 1]
dup 2 sub mul neg
1 shademin sub mul
shademin add
} def

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

% [pars] /fcn s0 s1 t0 t1 ns nt

/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

/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
[
  n 0 get n 1 get n 2 get 
  a 0 get n dot-product neg
]
end
} def

1 setlinejoin
1 setlinecap