pst-ode (27mar14)

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

1 files changed, 126 insertions, 83 deletions

diff --git a/Master/texmf-dist/dvips/pst-ode/pst-ode.pro b/Master/texmf-dist/dvips/pst-ode/pst-ode.pro
index 02ea9f726e5..e1de2a52680 100644
--- a/Master/texmf-dist/dvips/pst-ode/pst-ode.pro
+++ b/Master/texmf-dist/dvips/pst-ode/pst-ode.pro
@@ -1,7 +1,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
 % PostScript prologue for pst-ode.tex.
-% Version 0.6, 2013/11/26
+% Version 0.7, 2014/03/27
 %
 % Alexander Grahn (C) 2012, 2013
 %
@@ -12,102 +12,145 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 /tx@odeDict 1 dict def
 tx@odeDict begin
-%some constants for step size calculation
-/sfty 0.9 def /pgrow -0.2 def /pshrink -0.25 def
-%helper functions
-/addvect { % [1 2 3] [4 5 6] addvect => [5 7 9]
-  aload pop xlength1 -1 roll {xlength1 -1 roll add} forall
-  xlength array astore
-} bind def
-/subvect { % [1 2 3] [4 5 6] subvect => [-3 -3 -3]
-  aload pop xlength1 -1 roll {xlength1 -1 roll sub} forall
-  xlength array astore
-} bind def
-/mulvect { % [1 2 3] 4 mulvect => [4 8 12]
-  /mul cvx 2 array astore cvx forall xlength array astore
-} bind def
-/edivvect { % [1 2 3] [4 5 6] edivvect => [0.25 0.4 0.5]
-  aload pop xlength1 -1 roll {xlength1 -1 roll div} forall
-  xlength array astore
-} bind def
-/eabsvect { % [-1 2 -3] eabsvect => [1 2 3]
-  {abs} forall xlength array astore
-} bind def
-%/revstack { % (a) (b) (c) (d) 3 revstack => (a) (d) (c) (b)
-%  -1 2 {dup 1 sub neg roll} for
-%} bind def
-/min { 2 copy gt { exch } if pop } bind def
-/max { 2 copy lt { exch } if pop } bind def
-%coefficient table (Butcher table) of RKF45
-/a41 25 216 div def /a43 1408 2565 div def
-/a44 2197 4104 div def /a45 1 5 div neg def
-/a51 16 135 div def /a53 6656 12825 div def
-/a54 28561 56430 div def /a55 9 50 div neg def
-/a56 2 55 div def
-/b21 1 4 div def /b31 3 32 div def /b32 9 32 div def
-/b41 1932 2197 div def /b42 7200 2197 div neg def /b43 7296 2197 div def
-/b51 439 216 div def /b52 8 neg def /b53 3680 513 div def
-/b54 845 4104 div neg def /b61 8 27 div neg def /b62 2 def
-/b63 3544 2565 div neg def /b64 1859 4104 div def /b65 11 40 div neg def
+/ode@@dict 1 dict def
+/ode@dict {ode@@dict begin} bind def
+ode@dict
+  %some constants for step size calculation
+  /sfty 0.9 def /pgrow -0.2 def /pshrink -0.25 def
+  %helper functions
+  /addvect { % [1 2 3] [4 5 6] addvect => [5 7 9]
+    ode@dict
+      aload pop xlength1 -1 roll {xlength1 -1 roll add} forall
+      xlength array astore
+    end
+  } bind def
+  /subvect { % [1 2 3] [4 5 6] subvect => [-3 -3 -3]
+    ode@dict
+      aload pop xlength1 -1 roll {xlength1 -1 roll sub} forall
+      xlength array astore
+    end
+  } bind def
+  /mulvect { % [1 2 3] 4 mulvect => [4 8 12]
+    ode@dict /mul cvx 2 array astore cvx forall xlength array astore end
+  } bind def
+  /edivvect { % [1 2 3] [4 5 6] edivvect => [0.25 0.4 0.5]
+    ode@dict
+      aload pop xlength1 -1 roll {xlength1 -1 roll div} forall
+      xlength array astore
+    end
+  } bind def
+  /eabsvect { % [-1 2 -3] eabsvect => [1 2 3]
+    ode@dict {abs} forall xlength array astore end
+  } bind def
+  %/revstack { % (a) (b) (c) (d) 3 revstack => (a) (d) (c) (b)
+  %  -1 2 {dup 1 sub neg roll} for
+  %} bind def
+  /min { 2 copy gt { exch } if pop } bind def
+  /max { 2 copy lt { exch } if pop } bind def
+  %coefficient table (Butcher table) of RKF45
+  /a41 25 216 div def /a43 1408 2565 div def
+  /a44 2197 4104 div def /a45 1 5 div neg def
+  /a51 16 135 div def /a53 6656 12825 div def
+  /a54 28561 56430 div def /a55 9 50 div neg def
+  /a56 2 55 div def
+  /b21 1 4 div def /b31 3 32 div def /b32 9 32 div def
+  /b41 1932 2197 div def /b42 7200 2197 div neg def /b43 7296 2197 div def
+  /b51 439 216 div def /b52 8 neg def /b53 3680 513 div def
+  /b54 845 4104 div neg def /b61 8 27 div neg def /b62 2 def
+  /b63 3544 2565 div neg def /b64 1859 4104 div def /b65 11 40 div neg def
+end
 %Runge-Kutta-Fehlberg (RKF45) method
 %performs one integration step over tentative step size ddt
 %[state vector x(t)] RKF45 => [x(t)] [x(t+ddt) by RKF4] errmax
 /RKF45 {
-  /t tcur def
-  dup ODESET ddt mulvect /k1 exch def
-  /t tcur ddt 4 div add def
-  dup k1 b21 mulvect addvect ODESET ddt mulvect /k2 exch def
-  /t tcur ddt 3 mul 8 div add def
-  dup k1 b31 mulvect addvect k2 b32 mulvect addvect ODESET ddt mulvect
-    /k3 exch def
-  /t tcur ddt 12 mul 13 div add def
-  dup k1 b41 mulvect addvect k2 b42 mulvect addvect k3 b43 mulvect addvect
-    ODESET ddt mulvect /k4 exch def
-  /t tcur ddt add def
-  dup k1 b51 mulvect addvect k2 b52 mulvect addvect k3 b53 mulvect addvect
-    k4 b54 mulvect addvect ODESET ddt mulvect /k5 exch def
-  /t tcur ddt 2 div add def
-  dup k1 b61 mulvect addvect k2 b62 mulvect addvect k3 b63 mulvect addvect
-    k4 b64 mulvect addvect k5 b65 mulvect addvect ODESET ddt mulvect
-    /k6 exch def % => [x(t)]
-  %fourth order solution (increment dx)
-  dup dup k1 a41 mulvect k3 a43 mulvect addvect k4 a44 mulvect addvect
-    k5 a45 mulvect addvect dup
-    % => [x(t)] [x(t)] [x(t)] [dx by RKF4] [dx by RKF4]
-  %fifth order solution (abs. error)
-  k1 a51 mulvect k3 a53 mulvect addvect k4 a54 mulvect addvect
-    k5 a55 mulvect addvect k6 a56 mulvect addvect subvect
-    % => [x(t)] [x(t)] [x(t)] [dx by RKF4] [err]
-  5 1 roll addvect 4 -2 roll % => [x(t)] [x(t+ddt) by RKF4] [err] [x(t)]
-  %scaling vector for step size adjustment (Numerical Recipies)
-  eabsvect k1 eabsvect addvect {1e-30 add} forall xlength array astore
-  % => [x(t)] [x(t+ddt) by RKF4] [err] [xscale]
-  %maximum rel. error
-  /errmax 0 def
-  edivvect {abs dup errmax gt {/errmax exch def}{pop} ifelse} forall errmax
-  % => [x(t)] [x(t+ddt) by RKF4] errmax
+  dup
+  ode@dict tcur end /t exch def
+  ODESET
+  ode@dict
+    ddt mulvect /k1 exch def
+    dup k1 b21 mulvect addvect
+    tcur ddt 4 div add
+  end
+  /t exch def
+  ODESET
+  ode@dict
+    ddt mulvect /k2 exch def
+    dup k1 b31 mulvect addvect k2 b32 mulvect addvect
+    tcur ddt 3 mul 8 div add
+  end
+  /t exch def
+  ODESET
+  ode@dict
+    ddt mulvect /k3 exch def
+    dup k1 b41 mulvect addvect k2 b42 mulvect addvect k3 b43 mulvect addvect
+    tcur ddt 12 mul 13 div add
+  end
+  /t exch def
+  ODESET
+  ode@dict
+    ddt mulvect /k4 exch def
+    dup k1 b51 mulvect addvect k2 b52 mulvect addvect k3 b53 mulvect addvect
+    k4 b54 mulvect addvect
+    tcur ddt add
+  end
+  /t exch def
+  ODESET
+  ode@dict
+    ddt mulvect /k5 exch def
+    dup k1 b61 mulvect addvect k2 b62 mulvect addvect k3 b63 mulvect addvect
+    k4 b64 mulvect addvect k5 b65 mulvect addvect
+    tcur ddt 2 div add
+  end
+  /t exch def
+  ODESET
+  ode@dict
+    ddt mulvect /k6 exch def % => [x(t)]
+    %fourth order solution (increment dx)
+    dup dup k1 a41 mulvect k3 a43 mulvect addvect k4 a44 mulvect addvect
+      k5 a45 mulvect addvect dup
+      % => [x(t)] [x(t)] [x(t)] [dx by RKF4] [dx by RKF4]
+    %fifth order solution (abs. error)
+    k1 a51 mulvect k3 a53 mulvect addvect k4 a54 mulvect addvect
+      k5 a55 mulvect addvect k6 a56 mulvect addvect subvect
+      % => [x(t)] [x(t)] [x(t)] [dx by RKF4] [err]
+    5 1 roll addvect 4 -2 roll % => [x(t)] [x(t+ddt) by RKF4] [err] [x(t)]
+    %scaling vector for step size adjustment (Numerical Recipies)
+    eabsvect k1 eabsvect addvect {1e-30 add} forall xlength array astore
+    % => [x(t)] [x(t+ddt) by RKF4] [err] [xscale]
+    %maximum rel. error
+    /errmax 0 def
+    edivvect {abs dup errmax gt {/errmax exch def}{pop} ifelse} forall errmax
+    % => [x(t)] [x(t+ddt) by RKF4] errmax
+  end
 } bind def
 /ODEINT { % performs integration over output step [t,t+dt]
           % [ state vector x(t) ] ODEINT => [ state vector x(t+dt) ]
   %decrease overshooting step size
-  tcur ddt add tout sub tcur ddt add tout sub dt add mul 0 gt
-  {/ddt tout tcur sub def} if
+  ode@dict
+    tcur ddt add tout sub tcur ddt add tout sub dt add mul 0 gt
+    {/ddt tout tcur sub def} if
+  end
   RKF45
   ode@tol div dup 1 gt {
     %failed step -> reduce step size
-    exch pop pshrink exp 0.1 max sfty mul ddt mul /ddt exch def
-    (-) odeprint tcur ddt add tcur eq {
-      (error: step size underflow in ODEINT) print quit
-    } if
+    ode@dict
+      exch pop pshrink exp 0.1 max sfty mul ddt mul /ddt exch def
+      (-) odeprint ode@dict tcur ddt add tcur end eq {
+        (error: step size underflow in ODEINT) print quit
+      } if
+    end
     ODEINT %repeat step with new ddt
   }{
     %success
     3 -1 roll pop % remove previous state vector
-    /tcur tcur ddt add def
-    dup 0 ne {pgrow exp 5 min sfty mul ddt mul /ddt exch def}{pop} ifelse
-    %output step completed?
-    tcur tout sub dt mul 0 ge {
-      (o) odeprint /tcur tout def /tout tout dt add def
+    ode@dict
+      /tcur tcur ddt add def
+      dup 0 ne {pgrow exp 5 min sfty mul ddt mul /ddt exch def}{pop} ifelse
+      %output step completed?
+      tcur tout sub dt mul
+    end
+    0 ge {
+      (o) odeprint ode@dict /tcur tout def /tout tout dt add def end
     }{
       (+) odeprint ODEINT %continue integration
     } ifelse
@@ -123,7 +166,7 @@ tx@odeDict begin
 /assembleresult { %assembles state vector for building table of results
   {
     dup (t) eq {
-      pop tcur
+      pop ode@dict tcur end
     }{
       ode@laststate exch get
     } ifelse