diff options
Diffstat (limited to 'Master/texmf-dist/dvips/pst-3dplot/pst-3dplot.pro')
| -rw-r--r-- | Master/texmf-dist/dvips/pst-3dplot/pst-3dplot.pro | 381 |
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 + |