pstricks, psSolid and pstODEsolve

Question

When using pstODEsolve with the option algebraicOutputFormat I would like to use a function \z which itself is defined by the "final value" of a solution to another ode and the argument #1 being the time to be integrated over

\documentclass[border=0pt,10pt,dvipsnames]{standalone}
\usepackage{amsmath}
\usepackage{pstricks}
\usepackage{dashrule}
\usepackage{pst-plot}%, pst-math}
\usepackage{pst-func}
\usepackage{pst-ode}
\usepackage{pst-3dplot}
\usepackage{pst-solides3d}
\usepackage[off]{auto-pst-pdf}
\usepackage[nomessages]{fp}
\psset{algebraic,unit=2}
\usepackage[active,tightpage]{preview}
\PreviewBorder=12pt
\PreviewEnvironment{pspicture}


\makeatletter
\define@key[psset]{}{valuerange}[1-]{\expandafter\pst@getrange#1\@nil}
\def\pst@getrange#1-#2\@nil{%
  \ifx\relax#1\relax \def\pst@startvalue{1 }\else\def\pst@startvalue{#1 }\fi%
  \ifx\relax#2\relax \def\pst@endvalue{1e32 }\else\def\pst@endvalue{#2 }\fi%
}
\psset{valuerange=1-}

\def\listplotIIID{\def\pst@par{}\pst@object{listplotIIID}}
\def\listplotIIID@i#1{%
\@nameuse{beginplot@\psplotstyle}%
\addto@pscode{%
   /viewpointXYZ {\pst@solides@viewpoint} def
   /Decran \pst@solides@Decran\space def % distance de l'ecran
    viewpointXYZ /ZpointVue ED /YpointVue ED /XpointVue ED
  /THETA {YpointVue XpointVue atan} bind def
  /PHI   {ZpointVue XpointVue dup mul YpointVue dup mul add sqrt atan} bind def
  /Dobs  {XpointVue dup mul YpointVue dup mul add ZpointVue dup mul add sqrt} bind def
  XpointVue YpointVue ZpointVue /viewpoint defpoint3d
    /XYZ [#1] def
  /@tabXYZ [
    0 3 XYZ length 3 sub {
     /i exch def
     XYZ i get
     XYZ i 1 add get
     XYZ i 2 add get
     i 3 div dup /ii ED  
       \pst@startvalue ge {
         ii \pst@endvalue le { 
         3dto2d
         \pst@number\psunit mul exch
         \pst@number\psunit mul exch
       } { pop pop pop } ifelse
     }{ pop pop pop } ifelse
    }  for
    ] bind def
 [ @tabXYZ aload pop
    }%
\@nameuse{endplot@\psplotstyle}}%
\makeatother



\def\r(#1){sin(#1)}
\def\rs(#1){cos(#1)}
\def\rss(#1){-sin(#1)}
\def\rssrs(#1){-sin(#1)*cos(#1)}
\def\zsq(#1){sin(#1)*sin(#1)*(1+cos(#1)+sin(2*#1)*sin(2*#1)/16)}
\def\zsszs(#1){(1/2)*sin(2*#1)*(1+cos(#1)+sin(2*#1)*sin(2*#1)/16)+(1/2)*sin(#1)*sin(#1)*(-sin(#1)+sin(4*#1)/8)}
%%%%%
%%%%%   How to define a function as the "last" value of a numerical solution
%%%%%   to an ODE where the variable #1 is the time to be integrated over
%%%%%   starting from time = 0 in this example
%%%%%
%%%%%   the following is of course not defined!
%%%%%

\def\z(#1){ "last value of" \pstODEsolve[algebraic,algebraicT]{test}{1}{0}{#1}{50}{0 1}{ 1 | 
 -sin(x[0])*sqrt(1+cos(x[0])+sin(2*x[0])*sin(2*x[0])/16)} }


\begin{document}
\begin{pspicture}(-1.2,-1)(3.1,1.1)
\psset[pst-solides3d]{viewpoint=20 5 20,Decran=20,lightsrc=viewpoint}
\pstODEsolve[algebraic,algebraicOutputFormat]{test}{\r(x[2])*sin(x[0]) | \z(x[2]) 
| \r(x[2])*cos(x[0])}{0}{22}{200}{1.57 1 0.3 0}{ x[1] | -2*\rs(x[2])*x[1]*x[3]/(\r(x[2]))| x[3] | -(\rssrs(x[2])+\zsszs(x[2]))*x[3]*x[3]/(\rs(x[2])*\rs(x[2])+\zsq(x[2]))
+\rs(x[2])*\r(x[2])*x[1]*x[1]/(\rs(x[2])*\rs(x[2])+\zsq(x[2]))}
\listplotIIID[linecolor=black,linewidth=1pt]{test}
\axesIIID(1,4,1)(2,4.5,2)
\end{pspicture}
\end{document}

Answer 1

\pstODEsolve cannot be nested.

Fortunately, in the present case, we can work-around the problem, since the requested function z(x) is to be applied as a post-processing step on the second column of the solution of your main ODE.

This can be achieved in a two-step procedure.

1) solution of the main ODE without applying z(x) on the second element of the output vector:

\pstODEsolve[algebraic,algebraicOutputFormat,saveData]{intermediateResult}{%
  \r(x[2])*sin(x[0]) |
  x[2] | % we process this column later
  \r(x[2])*cos(x[0])
}{0}{22}{200}{1.57 1 0.3 0}{
 x[1] | -2*\rs(x[2])*x[1]*x[3]/(\r(x[2]))| x[3] |
 -(\rssrs(x[2])+\zsszs(x[2]))*x[3]*x[3]/(\rs(x[2])*\rs(x[2])+\zsq(x[2]))
 +\rs(x[2])*\r(x[2])*x[1]*x[1]/(\rs(x[2])*\rs(x[2])+\zsq(x[2]))
}

and write the solution table into a file intermediateResult.dat

2) Read the table from the intermediate file in a line-wise fashion, apply the other ODE

                                        % two output points--v
\pstODEsolve[algebraic,algebraicT]{processedValue}{1}{0}{#1}{2}{0 1}{
  1 | -sin(x[0])*sqrt(1+cos(x[0])+sin(2*x[0])*sin(2*x[0])/16)
}%

on the second element of the input line, take the last value from the result stored in processedValue and reassemble the table.

---

Since the PS --> PDF conversion step involves writing and reading text files, ps2pdf must be given the option -dNOSAFER.

Run: latex <...> dvips <...> ps2pdf -dNOSAFER <...>.ps

twice on the following input file:

\documentclass[border=10pt]{standalone}
\usepackage{pstricks}
\usepackage{pst-plot}%, pst-math}
\usepackage{pst-ode}
\usepackage{pst-solides3d}
\psset{algebraic,unit=2}

\makeatletter
\define@key[psset]{}{valuerange}[1-]{\expandafter\pst@getrange#1\@nil}
\def\pst@getrange#1-#2\@nil{%
  \ifx\relax#1\relax \def\pst@startvalue{1 }\else\def\pst@startvalue{#1 }\fi%
  \ifx\relax#2\relax \def\pst@endvalue{1e32 }\else\def\pst@endvalue{#2 }\fi%
}
\psset{valuerange=1-}

\def\listplotIIID{\def\pst@par{}\pst@object{listplotIIID}}
\def\listplotIIID@i#1{%
\@nameuse{beginplot@\psplotstyle}%
\addto@pscode{%
   /viewpointXYZ {\pst@solides@viewpoint} def
   /Decran \pst@solides@Decran\space def % distance de l'ecran
    viewpointXYZ /ZpointVue ED /YpointVue ED /XpointVue ED
  /THETA {YpointVue XpointVue atan} bind def
  /PHI   {ZpointVue XpointVue dup mul YpointVue dup mul add sqrt atan} bind def
  /Dobs  {XpointVue dup mul YpointVue dup mul add ZpointVue dup mul add sqrt} bind def
  XpointVue YpointVue ZpointVue /viewpoint defpoint3d
    /XYZ [#1] def
  /@tabXYZ [
    0 3 XYZ length 3 sub {
     /i exch def
     XYZ i get
     XYZ i 1 add get
     XYZ i 2 add get
     i 3 div dup /ii ED  
       \pst@startvalue ge {
         ii \pst@endvalue le { 
         3dto2d
         \pst@number\psunit mul exch
         \pst@number\psunit mul exch
       } { pop pop pop } ifelse
     }{ pop pop pop } ifelse
    }  for
    ] bind def
 [ @tabXYZ aload pop
    }%
\@nameuse{endplot@\psplotstyle}}%
\makeatother

\def\r(#1){sin(#1)}
\def\rs(#1){cos(#1)}
\def\rss(#1){-sin(#1)}
\def\rssrs(#1){-sin(#1)*cos(#1)}
\def\zsq(#1){sin(#1)*sin(#1)*(1+cos(#1)+sin(2*#1)*sin(2*#1)/16)}
\def\zsszs(#1){(1/2)*sin(2*#1)*(1+cos(#1)+sin(2*#1)*sin(2*#1)/16)+(1/2)*sin(#1)*sin(#1)*(-sin(#1)+sin(4*#1)/8)}

%for parsing lines of intermediateResult.dat
\def\getFirstOfThree#1 #2 #3;{#1}
\def\getSecondOfThree#1 #2 #3;{#2}
\def\getThirdOfThree#1 #2 #3;{#3}

\IfFileExists{intermediateResult.dat}{%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %  second pass:
  %
  %  post-process table --> z(second column)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % initialise Postscript variable to take the processed table
  \pstVerb{%
    true setglobal
    globaldict /processedTable {} put
    false setglobal
  }%
  \endlinechar=-1% suppress trailing space at input line end
  \newread\xIIfile%
  \openin\xIIfile=intermediateResult.dat%
  %line-wise read table from `intermediateResult.dat', treat second column
  \loop\read\xIIfile to \inputline%
  \ifeof\xIIfile\else%
    \pstODEsolve[algebraic,algebraicT]{processedVal}{1}{0}{\expandafter\getSecondOfThree\inputline;}{2}{0 1}{ 1 |
     -sin(x[0])*sqrt(1+cos(x[0])+sin(2*x[0])*sin(2*x[0])/16)}%
     %process input line and append result to postscript variable `processedTable' (executable list)
     \pstVerb{
       true setglobal
       globaldict /processedTable [
         processedTable
         \expandafter\getFirstOfThree\inputline;\space
         processedVal exch pop % insert second value and throw the first at t=0 away
         \expandafter\getThirdOfThree\inputline;] cvx put
       false setglobal
     }%
  \repeat%
}{%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %  first pass
  %
  %  solve main ODE without special treatment for second column of the solution table
  %  write solution to file `intermediateResult.dat'
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    \pstODEsolve[algebraic,algebraicOutputFormat,saveData]{intermediateResult}{%
      \r(x[2])*sin(x[0]) |
      x[2] | % we process this column later
      \r(x[2])*cos(x[0])
    }{0}{22}{200}{1.57 1 0.3 0}{
     x[1] | -2*\rs(x[2])*x[1]*x[3]/(\r(x[2]))| x[3] |
     -(\rssrs(x[2])+\zsszs(x[2]))*x[3]*x[3]/(\rs(x[2])*\rs(x[2])+\zsq(x[2]))
     +\rs(x[2])*\r(x[2])*x[1]*x[1]/(\rs(x[2])*\rs(x[2])+\zsq(x[2]))}%
}%

\begin{document}

\IfFileExists{intermediateResult.dat}{%
  %plot the resulting graph
  \begin{pspicture}(-1.2,-1)(3.1,1.1)
  \psset[pst-solides3d]{viewpoint=20 5 20,Decran=20,lightsrc=viewpoint}
  \listplotIIID[linecolor=black,linewidth=1pt]{processedTable}
  \axesIIID(1,4,1)(2,4.5,2)
  \end{pspicture}
}{
  dummy text
}  

\end{document}