pst-3dplot 1.78 (23dec07)

git-svn-id: svn://tug.org/texlive/trunk@5827 c570f23f-e606-0410-a88d-b1316a301751

1 files changed, 321 insertions, 60 deletions

diff --git a/Master/texmf-dist/dvips/pst-3dplot/pst-3dplot.pro b/Master/texmf-dist/dvips/pst-3dplot/pst-3dplot.pro
index 12f52ad1443..251f0612ff5 100644
--- a/Master/texmf-dist/dvips/pst-3dplot/pst-3dplot.pro
+++ b/Master/texmf-dist/dvips/pst-3dplot/pst-3dplot.pro
@@ -15,22 +15,24 @@
 %%   `pst-3dplot' is a PSTricks package to draw 3d curves and graphical objects
 %%
 %%
-%% version 0.23 / 2007-08-18  Herbert Voss <voss _at_ PSTricks.de>
+%% version 0.25 / 2007-12-22  Herbert Voss <voss _at_ PSTricks.de>
+%% with contributions of Darrell Lamm <darrell.lamm _at_ gtri.gatech.edu<
+%%            
 %
-/tx@3DPlotDict 50 dict def
+/tx@3DPlotDict 200 dict def
 tx@3DPlotDict begin
 %
 /printDot { gsave 2 copy 2 0 360 arc fill stroke grestore } def
 %
 /saveCoor { 
-  dzUnit mul /z exch def
-  dyUnit mul /y exch def
-  dxUnit mul /x exch def
+  dzUnit mul /z ED
+  dyUnit mul /y ED
+  dxUnit mul /x ED
 } def
 %
 /ConvertTo2D {
   RotatePoint
-  1 { %  dummy loop, will run only 1 time
+  1 { %  dummy loop, will run only 1 time, allows exit 
     coorType 0 le {                                               % the default |
       /x2D x leftHanded not { neg } if Alpha cos mul y Alpha sin mul add def %  /\  co system
       /y2D x leftHanded { neg } if Alpha sin mul y Alpha cos mul add neg Beta sin mul z Beta cos mul add def
@@ -43,6 +45,10 @@ tx@3DPlotDict begin
       /x2D y x 0.5 mul sub def
       /y2D z x 0.5 mul sub def 
       exit } if
+    coorType 3 le { % coorType |/_ with a 1/sqrt(2) shortend x-axis and 135 degrees 
+      /x2D y x -0.5 mul sub def
+      /y2D z x -0.5 mul sub def 
+      exit } if
   } repeat
 } def
 %
@@ -75,7 +81,7 @@ tx@3DPlotDict begin
 %
 /convertStackTo2D {
   counttomark
-  /n exch def /n3 n 3 div cvi def
+  /n ED /n3 n 3 div cvi def
   n3 {
     n -3 roll
     SphericalCoor { ConvertToCartesian } { saveCoor } ifelse
@@ -85,43 +91,181 @@ tx@3DPlotDict begin
   } repeat
 } def
 %
+/RotSet (set ) def
+%
+/eulerRotation false def
+% Matrix multiplication procedure
+/matmul {
+
+  /M@tMulDict 20 dict def
+  M@tMulDict begin
+  /m2 ED
+  /m1 ED
+  m1 dup length 2 sub 2 getinterval aload pop
+  /col1max ED
+  /row1max ED
+  m2 dup length 2 sub 2 getinterval aload pop
+  /col2max ED
+  /row2max ED
+  /m3 row1max col2max mul 2 add array def
+  m3 dup length 2 sub row1max col2max 2 array astore putinterval
+  0 1 row1max 1 sub {
+   /row ED
+   0 1 col2max 1 sub {
+    /col ED
+    /sum 0 def
+    0 1 col1max 1 sub{
+    /rowcol ED
+    sum
+    m1 row col1max mul rowcol add get
+    m2 rowcol col2max mul col add get
+    mul add 
+    /sum ED
+    } for
+    m3 row col2max mul col add sum put
+   } for
+  } for
+  m3
+  end % end of M@tMulDict
+
+} def
+%
+/SetMQuaternion {
+
+  /MnewTOold 11 array def
+
+  /Qu@ternionDict 30 dict def
+  Qu@ternionDict begin
+
+  /normRotVec  xRotVec yRotVec zRotVec 3 array astore VecNorm  def
+  normRotVec 0 gt
+  {/xRotVecNorm xRotVec normRotVec div def
+   /yRotVecNorm yRotVec normRotVec div def
+   /zRotVecNorm zRotVec normRotVec div def
+   RotAngle}
+  {/xRotVecNorm 1 def
+   /yRotVecNorm 0 def
+   /zRotVecNorm 0 def 
+   0} ifelse
+
+  2 div dup
+  /q0 exch cos def
+      sin dup dup
+  /q1 exch xRotVecNorm mul def
+  /q2 exch yRotVecNorm mul def
+  /q3 exch zRotVecNorm mul def
+
+  /q0q0 q0 q0 mul def
+  /q0q1 q0 q1 mul def
+  /q0q2 q0 q2 mul def
+  /q0q3 q0 q3 mul def
+
+  /q1q1 q1 q1 mul def
+  /q1q2 q1 q2 mul def
+  /q1q3 q1 q3 mul def
+
+  /q2q2 q2 q2 mul def
+  /q2q3 q2 q3 mul def
+
+  /q3q3 q3 q3 mul def
+
+  MnewTOold 0 q0q0 q1q1 add q2q2 sub q3q3 sub put
+  MnewTOold 1 q1q2 q0q3 sub 2 mul put
+  MnewTOold 2 q1q3 q0q2 add 2 mul put
+
+  MnewTOold 3 q1q2 q0q3 add 2 mul put
+  MnewTOold 4 q0q0 q1q1 sub q2q2 add q3q3 sub put
+  MnewTOold 5 q2q3 q0q1 sub 2 mul put
+
+  MnewTOold 6 q1q3 q0q2 sub 2 mul put
+  MnewTOold 7 q2q3 q0q1 add 2 mul put
+  MnewTOold 8 q0q0 q1q1 sub q2q2 sub q3q3 add put
+
+  MnewTOold 9 3 put
+  MnewTOold 10 3 put
+
+  end % end of Qu@ternionDict
+
+} def
+%
+/SetMxyz {
+
+  1.0 0.0 0.0  0.0 1.0 0.0  0.0 0.0 1.0  3 3  11 array astore /MnewTOold ED
+  RotSequence cvx exec % Now create a new MnewTOold using xyz, etc.
+
+} def
 %
-/CalcCoordinates {% from pst-vue3d
-    formulesTroisD
-    Xi 28.45 mul Yi 28.45 mul
+/ConcatMQuaternion {
+
+  MnewTOold % Push onto stack
+  SetMQuaternion % Uses [xyz]RotVec and RotAngle to make MnewToOld 
+  MnewTOold matmul /MnewTOold ED
+
+} def
+%
+/ConcatMxyz {
+
+  MnewTOold % Push onto stack
+  SetMxyz % Uses RotX, etc. to set MnewTOold 
+  MnewTOold matmul /MnewTOold ED
+
 } def
-% pour la 3D conventionnelle
-/formulesTroisD{%
-  /xObservateur x Sin1 mul neg y Cos1 mul add def
-  /yObservateur x Cos1Sin2 mul neg y Sin1Sin2 mul sub z Cos2 mul add def
-  /zObservateur x neg Cos1Cos2 mul y Sin1Cos2 mul sub z Sin2 mul sub def
-  /Xi xObservateur zObservateur Div def
-  /Yi yObservateur zObservateur Div def
+%
+/RotatePoint{
+  MnewTOold x y z  3 1  5 array astore matmul
+  0 3 getinterval aload pop 
+  /z ED 
+  /y ED 
+  /x ED 
 } def
 %
-/RotatePointXYZ{% Mxx are defined in the TeX file
-  /xi M11 x mul M12 y mul add M13 z mul add def
-  /yi M21 x mul M22 y mul add M23 z mul add def
-  /zi M31 x mul M32 y mul add M33 z mul add def
-  /x xi def
-  /y yi def
-  /z zi def
+/makeMoldTOnew {
+  /MoldTOnew 11 array def
+  MoldTOnew 0 MnewTOold 0 get put
+  MoldTOnew 1 MnewTOold 3 get put
+  MoldTOnew 2 MnewTOold 6 get put
+  MoldTOnew 3 MnewTOold 1 get put
+  MoldTOnew 4 MnewTOold 4 get put
+  MoldTOnew 5 MnewTOold 7 get put
+  MoldTOnew 6 MnewTOold 2 get put
+  MoldTOnew 7 MnewTOold 5 get put
+  MoldTOnew 8 MnewTOold 8 get put
+  MoldTOnew 9               3 put
+  MoldTOnew 10              3 put
 } def
 %
 /RotXaxis { 
-  /yTemp y RotX cos mul z RotX sin mul sub def
-  /z  y RotX sin mul z RotX cos mul add def
-  /y yTemp def
+  eulerRotation 
+  {1 0 0}
+  {makeMoldTOnew MoldTOnew  1 0 0  3 1  5 array astore matmul
+   0 3 getinterval aload pop} ifelse
+  /zRotVec ED
+  /yRotVec ED
+  /xRotVec ED
+  /RotAngle RotX def
+  ConcatMQuaternion
 } def
 /RotYaxis { 
-  /xTemp x RotY cos mul z RotY sin mul add def
-  /z  x RotY sin mul neg z RotY cos mul add def
-  /x xTemp def
+  eulerRotation 
+  {0 1 0}
+  {makeMoldTOnew MoldTOnew  0 1 0  3 1  5 array astore matmul
+   0 3 getinterval aload pop} ifelse
+  /zRotVec ED
+  /yRotVec ED
+  /xRotVec ED
+  /RotAngle RotY def
+  ConcatMQuaternion
 } def
 /RotZaxis { 
-  /xTemp x RotZ cos mul y RotZ sin mul sub def
-  /y  x RotZ sin mul y RotZ cos mul add def
-  /x xTemp def
+  eulerRotation 
+  {0 0 1}
+  {makeMoldTOnew MoldTOnew  0 0 1  3 1  5 array astore matmul
+   0 3 getinterval aload pop} ifelse
+  /zRotVec ED
+  /yRotVec ED
+  /xRotVec ED
+  /RotAngle RotZ def
+  ConcatMQuaternion
 } def
 /xyz { RotXaxis RotYaxis RotZaxis } def
 /yxz { RotYaxis RotXaxis RotZaxis } def
@@ -129,13 +273,13 @@ tx@3DPlotDict begin
 /xzy { RotXaxis RotZaxis RotYaxis } def
 /zxy { RotZaxis RotXaxis RotYaxis } def
 /zyx { RotZaxis RotYaxis RotXaxis } def
-%
-/RotatePoint { RotSequence cvx exec } def
+/quaternion { } def % Null
 %
 /VecNorm { 0 exch { dup mul add } forall sqrt } def
 %
 /UnitVec {			% on stack is [a]; returns a vector with [a][a]/|a|=1 
   dup VecNorm /norm ED
+  norm 0 lt {/norm 0 def} if
   { norm div } forall 3 array astore } def
 %
 /AxB {				% on the stack are the two vectors [a][b]
@@ -169,10 +313,10 @@ tx@3DPlotDict begin
 % les rayons de lumière
   xLight dup mul yLight dup mul zLight dup mul add add sqrt /NormeLight ED
 % the color values
-  /K exch def
-  /Yellow exch def
-  /Magenta exch def
-  /Cyan exch def
+  /K ED
+  /Yellow ED
+  /Magenta ED
+  /Cyan ED
 } def
 %
 /facetteSphere {
@@ -182,28 +326,28 @@ tx@3DPlotDict begin
   /Zpoint Rsphere phi sin mul CZ add def
   Xpoint Ypoint Zpoint tx@ProjThreeD moveto
   theta 1 theta increment add {%
-    /theta1 exch def
+    /theta1 ED
     /Xpoint Rsphere theta1 cos mul phi cos mul CX add def
     /Ypoint Rsphere theta1 sin mul phi cos mul CY add def
     /Zpoint Rsphere phi sin mul CZ add def
     Xpoint Ypoint Zpoint tx@ProjThreeD  lineto
   } for
   phi 1 phi increment add {
-    /phi1 exch def
+    /phi1 ED
     /Xpoint Rsphere theta increment add cos mul phi1 cos mul CX add def
     /Ypoint Rsphere theta increment add sin mul phi1 cos mul CY add def
     /Zpoint Rsphere phi1 sin mul CZ add def
     Xpoint Ypoint Zpoint tx@ProjThreeD lineto
   } for
   theta increment add -1 theta {%
-    /theta1 exch def
+    /theta1 ED
     /Xpoint Rsphere theta1 cos mul phi increment add cos mul CX add def
     /Ypoint Rsphere theta1 sin mul phi increment add cos mul CY add def
     /Zpoint Rsphere phi increment add sin mul CZ add def
     Xpoint Ypoint Zpoint tx@ProjThreeD lineto
   } for
   phi increment add -1 phi {
-    /phi1 exch def
+    /phi1 ED
     /Xpoint Rsphere theta cos mul phi1 cos mul CX add def
     /Ypoint Rsphere theta sin mul phi1 cos mul CY add def
     /Zpoint Rsphere phi1 sin mul CZ add def
@@ -215,17 +359,17 @@ tx@3DPlotDict begin
 /MaillageSphere { 
 % on stack must be x y z Radius increment C M Y K 
   setColorLight
-  /increment exch def
-  /Rsphere exch def
-  /CZ exch def
-  /CY exch def
-  /CX exch def
+  /increment ED
+  /Rsphere ED
+  /CZ ED
+  /CY ED
+  /CX ED
   /StartTheta 0 def
   /condition { PSfacetteSphere 0 ge } def
   -90 increment 90 increment sub {%
-    /phi exch def
+    /phi ED
     StartTheta increment 360 StartTheta add increment sub {%
-      /theta exch def
+      /theta ED
       % Centre de la facette
       /Xpoint Rsphere theta increment 2 div add cos mul phi increment 2 div add cos mul CX add def
       /Ypoint Rsphere theta increment 2 div add sin mul phi increment 2 div add cos mul CY add def
@@ -267,7 +411,7 @@ tx@3DPlotDict begin
 /PlanCoupeCylinder { %
   /TableauxPoints [
     0 1 359 { 
-      /phi exch def 
+      /phi ED 
       [ Radius phi Height ConvCyl2d ] % on décrit le cercle
     } for
   ] def
@@ -374,7 +518,7 @@ tx@3DPlotDict begin
   /height ED
   /radius ED
   startAngle increment endAngle {
-    /Angle exch def
+    /Angle ED
     radius Angle 0 ConvCylToCartesian MoveTo  
     radius Angle height ConvCylToCartesian LineTo  
   } for
@@ -494,13 +638,13 @@ tx@3DPlotDict begin
     /V {Z sqrt} bind def
     /TableauxPoints [
       0 1 359 { 
-        /U exch def [ U U Z V calculate2DPoint ] % on décrit le cercle
+        /U ED [ U U Z V calculate2DPoint ] % on décrit le cercle
       } for
     ] def
     newpath
     TableauxPoints 0 get aload pop moveto
     0 1 359 {
-      /compteur exch def
+      /compteur ED
       TableauxPoints compteur get aload pop
       lineto } for
     closepath
@@ -510,20 +654,20 @@ tx@3DPlotDict begin
     newpath
     U U Z V calculate2DPoint moveto
     U 1 U increment add  {%
-      /U1 exch def
+      /U1 ED
       U1 U1 Z V calculate2DPoint lineto
     } for
     Z pas10 Z pas add pas10 add{
-      /Z1 exch def
+      /Z1 ED
       /V {Z1 sqrt} bind def
       U1 U1 Z1 V calculate2DPoint lineto
     } for
     U increment add -1 U {%
-      /U2 exch def
+      /U2 ED
       U2 U2 Z pas add V calculate2DPoint lineto
     } for
     Z pas add pas10 sub pas10 neg Z pas10 sub {
-      /Z2 exch def
+      /Z2 ED
       /V Z2 abs sqrt def
       U U Z2 V calculate2DPoint lineto
     } for
@@ -539,7 +683,7 @@ tx@3DPlotDict begin
       /Z ED
       /V Z sqrt def
       0 increment 360 increment sub {%
-        /U exch def
+        /U ED
 % Centre de la facette
         /Ucentre U increment 2 div add def
         /Vcentre Z pas 2 div add sqrt def
@@ -576,5 +720,122 @@ tx@3DPlotDict begin
     } for
 } def
 %
+% ------------------------------------ math stuff ----------------------------------
+%
+% Matrix A in arrays of rows A[[row1][row2]...]
+% with [row1]=[a11 a12 ... b1]
+% returns on stack solution vector X=[x1 x2 ... xn]
+/SolveLinEqSystem { 				% on stack matrix M=[A,b] (A*x=b)
+  10 dict begin					% hold all ocal
+    /A exch def
+    /Rows A length def         			% Rows = number of rows
+    /Cols A 0 get length def   			% Cols = number of columns
+    /Index [ 0 1 Rows 1 sub { } for ] def	% Index = [0 1 2 ... Rows-1]
+    /col 0 def
+    /row  0 def
+    /PR Rows array def 				% PR[c] = pivot row for row row
+  { 						% starts the loop, find pivot entry in row r
+    col Cols ge row  Rows ge or { exit } if	% col < Cols and row < Rows else exit
+    /pRow row def  				% pRow = pivot row		
+    /max A row  get col get abs def		% get A[row[col]], first A[0,0] 
+    row 1 add 1 Rows 1 sub { 			% starts for loop 1 1 Rows-1
+      /j exch def				% index counter
+      /x A j get col get abs def		% get A[j[r]]
+      x max gt {				% x>max, then save position
+        /pRow j def
+        /max x def
+      } if
+    } for					% now we have the row with biggest A[0,1]
+						% with pRow = the pivot row
+    max 0 gt {					% swap entries pRow and row  in i 
+      /tmp Index row  get def
+      Index row  Index pRow get put
+      Index pRow tmp put			% and columns pRow and row  in A
+      /tmp A row get def
+      A row  A pRow get put
+      A pRow tmp put   				% pivot
+      /row0  A row  get def 			% the pivoting row
+      /p0 row0  col get def 			% the pivot value
+      row 1 add 1 Rows 1 sub { 			% start for loop
+        /j exch def
+        /c1 A j get def
+        /p c1 col get p0 div def
+        c1 col p put				% subtract (p1/p0)*row[i] from row[j]
+        col 1 add 1 Cols 1 sub {		% start for loop
+          /i exch def
+          c1 dup i exch 			% c1 i c1
+          i get row0 i get p mul sub put
+        } for
+      } for
+      PR row col put
+      /col col 1 add def
+      /row row 1 add def
+    }{						% all zero entries
+      /row row 1 add def			% continue loop with same row
+    } ifelse
+  } loop
+  /X A def					% solution vector
+  A Rows 1 sub get dup
+  Cols 1 sub get exch
+  Cols 2 sub get div
+  X Rows 1 sub 3 -1 roll put  			% X[n]
+  Rows 2 sub -1 0 {				% for loop to calculate X[i]
+    /xi exch def				% current index
+    A xi get 					% i-th row
+    /Axi exch def
+    /sum 0 def
+    Cols 2 sub -1 xi 1 add { 
+      /n exch def
+      /sum sum Axi n get X n get mul add def 
+    } for
+    Axi Cols 1 sub get 				% b=Axi[Cols-1]
+    sum sub 					% b-sum
+    Axi xi get div				% b-sum / Axi[xi]
+    X xi 3 -1 roll put  			% X[xi]
+  } for
+  X
+  end 
+} def
 %
-end
+% 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
+%
+% [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
+%
+%
+end % tx@3DPlotDict
+