% $Id: tutFromScratch.tex 1710 2008-01-29 04:59:25Z gmcmanus $

% Copyright (C) 2000-2007
%
% Code contributed by Greg Collecutt, Joseph Hope and the xmds-devel team
%
% This file is part of xmds.
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software 
% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.

\chapter{Starting from scratch}
\label{chap:tutFromScratch}

It turns out that one of the most difficult things to do in \xmds is
to write a script from scratch, which is why most users take either
one of their own old scripts, or borrow someone else's (or one of the
examples) as a template.  However, this doesn't mean to say that one
is never going to have to write a script from scratch (just that it's
usually quite unlikely), therefore we explain how to do this here.  If
you're not interested, and want to get coding more quickly, then it's
alright to skip to \Chap{chap:usingATemplate} and learn how to
modify an existing script to your needs.

In this chapter we will outline the tags necessary for \xmds to
actually parse a document, and will implement a simple simulation
involving these tags.  More complete documentation of the tags and
what they do can be found in \Chap{chap:languageRef}.  The tags
necessary for performing more advanced operations and for solving more
complex problems will be introduced in later tutorial chapters
(e.g.~see \Chap{chap:stochasticSimsAndMPI}).

\section{A basic simulation}
\label{sec:basicSim}

\subsection{The simple beginnings of a simulation}

\xmds is coded using XML (the extensible markup language), and as such
each \xmds script must be a ``properly formed'' XML document.  One of the
main stipulations for an XML document to be properly formed is
that it have the following line at the top of each document, and so,
each \xmds script \tbf{must} have this as its first line:
\begin{xmdsCode}
<?xml version="1.0"?>
\end{xmdsCode}
There are other stipulations, but they don't really concern us too
much at the moment.  If you're interested, you can check out the World
Wide Web Consortium (W3C) web site
(\htmladdnormallink{http://www.w3c.org}{http://www.w3c.org}) for more
information.

Each \xmds simulation is enclosed within a set of
\xmdsTag{simulation} tags.  Therefore, for each simulation that you
write from scratch, the first few lines of code are going to be:
\begin{xmdsCode}
<?xml version="1.0"?>
<simulation>

  <!-- yet more xmds code to come -->

</simulation>
\end{xmdsCode}
It might be a good idea at this stage to
save this code to file.  So, for the purposes of the tutorial, we'll
refer to this script as \ttt{lorenz.xmds}.

Next, the \xmds script needs a name.  Well, to be honest, it doesn't
really need a name, but it's a good idea to add the \xmdsTag{name} tag
if only for completeness.  If the \xmdsTag{name} tag is not included,
\xmds uses the name of the \xmds script file (but with the \ttt{.xmds}
extension removed) as the name of the simulation, so in this
simulation if we were to not specify the \xmdsTag{name} then the
simulation name used by \xmds would be \ttt{lorenz}.  However, it is
still a good idea to specify the name inside the simulation text and
we do this here, setting \xmdsTag{name} to \ttt{lorenz}.  The script
now becomes:
\begin{xmdsCode}
<?xml version="1.0"?>
<simulation>

  <name>lorenz</name>

  <!-- yet more xmds code to come -->

</simulation>
\end{xmdsCode}
Notice that we have indented the \xmdsTag{name} tag from the rest of
the tags.  It is a good idea to indent tags that are nested within
other tags so that one gets an idea of the document structure just
from looking at the source code.

The next two tags that I recommend you add to your scripts are the
\xmdsTag{author} and \xmdsTag{description} tags.  These tags aren't
necessary for running a simulation or for actually calculating
anything, however, they are good for documenting the simulation, and
providing extra information that may be helpful to others if you show
your scripts to other people.  Also, this information may actually be
used in future versions of \xmds to provide extra functionality in the
output from \xmds simulations.  The code now is:
\begin{xmdsCode}
<?xml version="1.0"?>
<simulation>

  <name>lorenz</name>
  <author>Paul Cochrane</author>
  <description>
    Lorenz attractor example simulation.  Adapted from the example
    in "Numerical methods for physics" by Alejandro L. Garcia, 
    page 78 (1st ed.).
  </description>

  <!-- yet more xmds code to come -->
  
</simulation>
\end{xmdsCode}
Alternatively, you can add such information to the script by putting
comments in the source code.  The above code may then look something
like this:
\begin{xmdsCode}
<?xml version="1.0"?>
<!-- Example simulation: lorenz -->

<!-- Adapted from the example in "Numerical methods for physics"-->
<!-- by Alejandro L. Garcia, pg 78-->

<!-- Xmds script by: Paul Cochrane -->

<simulation>

  <name>lorenz</name>
  
  <!-- yet more xmds code to come -->
  
</simulation>
\end{xmdsCode}

\subsection{General simulation options}

Now that the very basic preliminaries are out of the way, we need to
focus on the problem we are trying to solve.  So, for the sake of
argument, let's try to solve the Lorenz equations,
\begin{align}
  \frac{dx}{dt} &= \sigma (y - x)\\
  \frac{dy}{dt} &= rx - y - xz\\
  \frac{dz}{dt} &= xy - bz,
  \label{eq:LorenzEquations}
\end{align}
where $\sigma$, $r$ and $b$ are positive constants, and the variables
have the initial conditions: $x(t=0) = x_0$, $y(t=0) = y_0$, $z(t=0) =
z_0$.

Even though we are trying to solve something as complex as a chaotic
system, this is very easy to write down in \xmds.  This is especially
true since the derivatives are with respect to time, for if we had
spatial derivatives we would have to use Fourier transform techniques
and mappings to simplify the calculation (we'll cover this stuff in
\Sec{sec:moreComplexSimulation}, so don't worry that we're not
discussing it here) which would complicate our script a bit more, and
we're trying to keep things simple here.

The Lorenz model can be used in the study of many interesting
phenomena, however it is possibly best known as a model of global
weather~\cite{Garcia:1994:1}.  For a system to be chaotic it must be
extremely sensitive to initial conditions such that any small
perterbation of the initial conditions will cause wildly divergent
evolution; and this is something we will hopefully see here, so let's continue.

There can be many global parameters in an \xmds simulation, although,
one in particular is special.  This is the propagation direction,
specified by the \xmdsTag{prop\_dim} tag, which is specified as part of
the global functionality of the \xmds simulation and appears within
the \xmdsTag{simulation} tags, normally after the
\xmdsTag{description}.  In the problem we are trying to solve, our
field is evolving in \emph{time}, given by the variable $t$ in the
above equations, and so the field is said to be propagating in $t$ so
we use time as the propagation dimension.  Therefore, we add the line
\begin{xmdsCode}
<prop_dim>t</prop_dim>
\end{xmdsCode}
to our \xmds script, just after the \xmdsTag{description} or
conversely just before the \xmdsTag{globals} section (in the situation
considered here, these locations are one in the same, however, this is
not the case in general).  The script is now
\begin{xmdsCode}
<?xml version="1.0"?>
<simulation>

  <name>lorenz</name>
  <author>Paul Cochrane</author>
  <description>
    Lorenz attractor example simulation.  Adapted from the example
    in "Numerical methods for physics" by Alejandro L. Garcia, 
    page 78 (1st ed.).
  </description>

  <prop_dim>t</prop_dim>
  
  <!-- yet more xmds code to come -->
  
</simulation>
\end{xmdsCode}

This problem we are trying to solve has several constants, namely
$\sigma$, $r$, $b$, $x_0$, $y_0$, and $z_0$.  These variables are
going to be used again and again in the simulation therefore it makes
sense to put them into the next element necessary to describe a
simulation in \xmds, the \xmdsTag{globals} tag.  We specify these
constants using C/C++ syntax in an XML \ttt{CDATA} block, which is
enclosed within the \xmdsTag{globals} element.  The \xmds code for
this is
\begin{xmdsCode}
<globals>
<![CDATA[
    const double sigma = 10.0;
    const double b = 8.0/3.0;
    const double r = 28.0;
    const double xo = 1.0;     // initial conditions
    const double yo = 1.0;     // initial conditions
    const double zo = 20.0;    // initial conditions
]]>
</globals>
\end{xmdsCode}
There are a couple of points to note here.  Firstly, and most
importantly, the syntax of the code within the \ttt{CDATA} block
\tbf{must} conform to C/C++ syntax rules, otherwise the simulation
won't be able to be compiled.  This is because the code within the
\ttt{CDATA} block is inserted directly into the code for the output
simulation.  Secondly, the constants that we are specifying are
declared to be double precision to \xmds (this is via the \ttt{double}
keyword in the code above) since they are continuous variables in the
problem being solved.  If for instance, we had some discrete quantity
such as the number of particles in a system, then we would specify the
variable as being integer and would therefore use the \ttt{int}
keyword to declare the variable as such.  Lastly, the odd-looking tags
for the \ttt{CDATA} block must be written correctly (i.e.~opened with
\verb|<![CDATA[| and closed with \verb|]]>|) otherwise the file won't
parse, and \xmds will give an error.  Our script is now,
\begin{xmdsCode}
<?xml version="1.0"?>
<simulation>

  <name>lorenz</name>
  <author>Paul Cochrane</author>
  <description>
    Lorenz attractor example simulation.  Adapted from the example
    in "Numerical methods for physics" by Alejandro L. Garcia, 
    page 78 (1st ed.).
  </description>

  <prop_dim>t</prop_dim>
  
  <globals>
  <![CDATA[
      const double sigma = 10.0;
      const double b = 8.0/3.0;
      const double r = 28.0;
      const double xo = 1.0;     // initial conditions
      const double yo = 1.0;     // initial conditions
      const double zo = 20.0;    // initial conditions
  ]]>
  </globals>
  
  <!-- yet more xmds code to come -->
  
</simulation>
  \end{xmdsCode}

\subsection{The field element}
\label{sec:theFieldElement}

Now that we have set up the physical constants of our problem, we need
to describe the field that we are going to be integrating over, in
other words we need to specify a discretised version of $x(t)$, $y(t)$
and $z(t)$ at the start of the problem.  To specify the field for the
simple example we are studying here, we only need to tell \xmds the
initial value of the field (i.e.~that when $t=0$).  Note that in more
complex situations (to be studied later) there are more tags to worry
about.

%% To properly specify the field we must
%% describe the dimensions that it is integrated across, the number of grid
%% points we are going to be using in that dimension, the domain over
%% which the field is to be solved, and a description of the initial
%% value of the field over the relevant dimensions.  To this end we now
%% introduce several new tags with which to specify these quantities.

The first tag is the \xmdsTag{field} element.  This is a container for
the other information that we are using to describe the field, and is
used very simply as follows
\begin{xmdsCode}
<field>

  <!-- More xmds tags in here -->
  
</field>
\end{xmdsCode}

We next give the name of the field, this is supplied with the
\xmdsTag{name} tag (note that this is a sub-element of the
\xmdsTag{field} tag, and so is different from the \xmdsTag{name} tag
of the \xmdsTag{simulation} element).  If no name is given for the
field, it defaults to ``\ttt{main}'', however, as mentioned before, it
is a good idea to specify a name for the field anyway.  We'll call it
\ttt{main} to be consistent with other scripts you are likely to see,
and because this is the main field.

%% The \xmdsTag{dimensions} tag usually comes next.  This tells \xmds the
%% name(s) of the dimension(s) that the field is to be integrated across
%% (as opposed to the propagation direction, which in this case is time
%% (\ttt{t}).

%% Computer simulations cannot model continuous fields exactly, hence one
%% must make an approximation to the field on a finite-sized discrete
%% lattice.  To tell \xmds the number of points we want to use to
%% approximate our field we specify this number within the
%% \xmdsTag{lattice} element.  For speed related reasons, especially when
%% one is using Fourier transforms, it is a good idea to specify the
%% \xmdsTag{lattice} element as a factor of two.  Since we only have a
%% finite-sized domain on which to simulate the evolution of the field,
%% we must specify this as well.  One does this by giving the end points
%% of the domain, comma-separated within paretheses (i.e.~an ordered
%% pair) within the \xmdsTag{domain} element, like so
%% \begin{xmdsCode}
%%     <domains>(-0.5,0.5)</domains>
%%   \end{xmdsCode}

We next must tell \xmds which of the field moments to sample
\emph{directly after} the field is initialised.  This is not obvious
to those new to \xmds, however this gives one the opportunity to
choose whether or not to sample the initial point of the field, and
generate output moments such as mean and standard error, before
it is propagated.  To do this we put a sequence of \ttt{1}'s or
\ttt{0}'s within the \xmdsTag{samples} element.  We must put as many
\ttt{1}'s and \ttt{0}'s as there are moment groups defined in the
\xmdsTag{output} tag (to be discussed later).  In our case, we will
only be using one output moment group here, and we do want to sample
the initial field for the moments, hence we add the code
\begin{xmdsCode}
<samples>1</samples>
\end{xmdsCode}
to our script.  The simulation at this stage looks like
\begin{xmdsCode}
<?xml version="1.0"?>
<simulation>

  <name>lorenz</name>
  <author>Paul Cochrane</author>
  <description>
    Lorenz attractor example simulation.  Adapted from the example
    in "Numerical methods for physics" by Alejandro L. Garcia, 
    page 78 (1st ed.).
  </description>

  <prop_dim>t</prop_dim>
  
  <globals>
  <![CDATA[
      const double sigma = 10.0;
      const double b = 8.0/3.0;
      const double r = 28.0;
      const double xo = 1.0;     // initial conditions
      const double yo = 1.0;     // initial conditions
      const double zo = 20.0;    // initial conditions
  ]]>
  </globals>
  
  <field>
    <name>main</name>
    <samples>1</samples>
    
    <!-- yet more xmds code to come -->
    
  </field>

</simulation>
\end{xmdsCode}

Now we have to initialise the field.  In other words, we have to
define $x(t=0)$, $y(t=0)$, and $z(t=0)$.  \xmds makes this easy by
allowing you to define the field as a vector written in terms of the
dimensions of the field.  It is possible to define other vectors that
are part of the field, but the vector of the field that we are
integrating is the \emph{main} vector, and this we name (funnily
enough) \ttt{main}.  This naming is compulsory since it is possible to
have more than one vector named within a field; however, there must be
exactly one vector named \ttt{main}.  We also need to
specify the data type of our vector, the names of the components of
the vector and we need to define how the vector should be calculated
using C code (in a \ttt{CDATA} section).  For the case we are considering
here, the \xmds code would be:
\begin{xmdsCode}
<vector>
  <name> main </name>
  <type> double </type>
  <components> x y z </components>
  <![CDATA[
      x = xo;
      y = yo;
      z = zo;
  ]]>
</vector>
  \end{xmdsCode}
So, as we can see, our vector is called \ttt{main}, it is of type
\ttt{double} and the names of its components are \ttt{x}, \ttt{y}, and
\ttt{z}.  The \ttt{CDATA} section gives the C code version of what
\eqn{eq:LorenzEquations} describes.

This code completes the \xmdsTag{field} element and we are left with the
following code listing:
\begin{xmdsCode}
<?xml version="1.0"?>
<simulation>

  <name>lorenz</name>
  <author>Paul Cochrane</author>
  <description>
    Lorenz attractor example simulation.  Adapted from the example
    in "Numerical methods for physics" by Alejandro L. Garcia, 
    page 78 (1st ed.).
  </description>

  <prop_dim>t</prop_dim>
  
  <globals>
  <![CDATA[
      const double sigma = 10.0;
      const double b = 8.0/3.0;
      const double r = 28.0;
      const double xo = 1.0;     // initial conditions
      const double yo = 1.0;     // initial conditions
      const double zo = 20.0;    // initial conditions
  ]]>
  </globals>
  
  <field>
    <name>main</name>
    <samples>1</samples>
    
    <vector>
      <name> main </name>
      <type> double </type>
      <components> x y z </components>
      <![CDATA[
	  x = xo;
	  y = yo;
	  z = zo;
      ]]>
    </vector>
      
  </field>

  <!-- yet more xmds code to come -->
    
</simulation>
  \end{xmdsCode}

\subsection{The sequence element}

We have arrived at the stage where we can tell \xmds how to actually
perform the integration of the field.  To do this we use the
\xmdsTag{sequence} element.  The \xmdsTag{sequence} element is usually
used as a container for other elements, specifically the
\xmdsTag{integrate}, \xmdsTag{filter} and other \xmdsTag{sequence}
elements.  The outermost \xmdsTag{sequence} element is referred to as
the ``parent'' \xmdsTag{sequence} and the \xmdsTag{sequence}s nested
within that as the ``child'' \xmdsTag{sequence}s.  This may sound a
bit confusing, but it is just a generalisation and a lot of the time
you will be writing scripts with just the one \xmdsTag{sequence}.  The
\xmdsTag{sequence} element may have as many of the other sub-elements
as desired to perform the calculation, and the ``child''
\xmdsTag{sequence}s can contain another element---the
\xmdsTag{cycles} element---which controls how many times a given
\xmdsTag{sequence} is repeated.  The \xmdsTag{cycles} element is
optional and defaults to one.  It is important to note that the order
of segments specified within a \xmdsTag{sequence} are significant, and
operations given will be performed in that order.

So, to summarise, most of the time you will just use one
\xmdsTag{sequence} and it will usually only contain just the one
\xmdsTag{integrate} section, hence the code will look like
\begin{xmdsCode}
<sequence>
  <integrate>
  
    <!-- More xmds tags to come -->
    
  </integrate>
</sequence>
  \end{xmdsCode}
We shall try to discuss the other features and tags in more depth later
on in more advanced tutorials.

\subsection{The integrate element}

Since the \xmdsTag{integrate} element is quite complex, and it does
all of the hard work, we'll spend some time discussing it.

We need to tell \xmds the algorithm to use to integrate the field
specified earlier.  To do this we use the \xmdsTag{algorithm} tag.
This tag is optional and will default to \ttt{SIEX} for stocastic
simulations and to \ttt{RK4EX} for non-stochastic simulations,
however, it is a very good idea to explicitly specify what algorithm
your simulation is using, if only to help yourself in six months time,
or a colleague who may end up reading your code.  At the moment (\xmds
version 1.5-1) there are six algorithms to choose from:
\ttt{RK4EX}, \ttt{RK4IP}, \ttt{ARK45EX}, \ttt{ARK45IP}, \ttt{SIEX}, \ttt{SIIP}.  \ttt{RK4EX} is a
fourth order Runge-Kutta in the explicit picture, \ttt{RK4IP} is a
fourth order Runge-Kutta in the interaction picture, the \ttt{ARK45EX} and \ttt{ARK45IP} are the corresponding adaptive time step Runge-Kutta Fehlberg methods, \ttt{SIEX} is the semi-implicit method in the explicit picture, and \ttt{SIIP} is the semi-implicit method in the interaction picture.  
For more information
about the specifics of these algorithms and techniques, see
\Sec{sec:numericalMethods}.  In solving our problem we'll use the
fourth order Runge-Kutta in the explicit picture, because we aren't
using any Fourier transforms, and the explicit picture is fine for our
purposes here.  Therefore, we specify within the
\xmdsTag{integrate} tags the line:
\begin{xmdsCode}
<algorithm>RK4EX</algorithm>
  \end{xmdsCode}
telling \xmds what algorithm to use.  

The next things \xmds needs to know are the length of the integration
interval, the total number of steps to take, and the number of samples
for each output moment to take within these steps.  These items are
denoted by the \xmdsTag{interval}, \xmdsTag{lattice} and
\xmdsTag{samples} tags respectively.  The integration interval
combined with the number of steps gives the step size internally used
by \xmds.  We'll choose some fairly arbitrary numbers here: an
interval length of 10, a large number of lattice points, namely 10000,
(which should give us a nice small step size), and we'll sample 200
points, and hence set \xmdsTag{samples} to 200.  The code for this looks like:
\begin{xmdsCode}
<interval> 10 </interval>
<lattice> 10000 </lattice>
<samples> 200 </samples>
  \end{xmdsCode}

Having told \xmds some of the parameters it must use to perform the
integration, we haven't yet told it \emph{how} to actually carry out
the integration.  By this I mean that we have to describe in terms of
C language code the differential equation that \xmds is to use to
evolve the solution forward (see \eqn{eq:LorenzEquations}).  We do
this by using a \ttt{CDATA} block, and writing the equation in a form
understandable to both us and \xmds which isn't technically speaking C
code.
\begin{xmdsCode}
<![CDATA[
    dx_dt = sigma*(y - x);
    dy_dt = r*x - y - x*z;
    dz_dt = x*y - b*z;
]]>
  \end{xmdsCode}
Notice that this looks similar to the analytical form of
\eqn{eq:LorenzEquations} in that we have described the derivatives
with respect to time, $t$ of the fields $x(t)$, $y(t)$ and $z(t)$.

With that, we have completed the \xmdsTag{integrate} element, and the
\xmdsTag{sequence} section.  The simulation script is now:
\begin{xmdsCode}
<?xml version="1.0"?>
<simulation>

  <name>lorenz</name>
  <author>Paul Cochrane</author>
  <description>
    Lorenz attractor example simulation.  Adapted from the example
    in "Numerical methods for physics" by Alejandro L. Garcia, 
    page 78 (1st ed.).
  </description>

  <prop_dim>t</prop_dim>
  
  <globals>
  <![CDATA[
      const double sigma = 10.0;
      const double b = 8.0/3.0;
      const double r = 28.0;
      const double xo = 1.0;     // initial conditions
      const double yo = 1.0;     // initial conditions
      const double zo = 20.0;    // initial conditions
  ]]>
  </globals>
  
  <field>
    <name>main</name>
    <samples>1</samples>
    
    <vector>
      <name> main </name>
      <type> double </type>
      <components> x y z </components>
      <![CDATA[
	  x = xo;
	  y = yo;
	  z = zo;
      ]]>
    </vector>
  </field>

  <sequence>
    <integrate>
      <algorithm>RK4EX</algorithm>
      <interval> 10 </interval>
      <lattice> 10000 </lattice>
      <samples> 200 </samples>
      <![CDATA[
          dx_dt = sigma*(y - x);
          dy_dt = r*x - y - x*z;
	  dz_dt = x*y - b*z;
      ]]>
    </integrate>
  </sequence>

  <!-- yet more xmds code to come -->

</simulation>
  \end{xmdsCode}

\subsection{The output element}

Ok, we're almost there!  This is the last section that we need to
worry about.  So, just to recap, we've told \xmds the general features
and variables it should use to construct the simulation, the field to
integrate over, and the way in which the integration should take place
and most importantly the differential equation that \xmds should use
to evolve the solution.  That sounds like about it doesn't it?  Well,
no.  We haven't told \xmds to output anything yet, and it's a bit
silly to spend several hours of computer time for the simulation to
come back to you and say: ``I'm done!'' and you haven't got any
results.  Fortunately, \xmds doesn't let you write a simulation
without actually specifying any output.  Therefore, we need to tell
\xmds what output we want from the simulation, and to do this we use
the \xmdsTag{output} element.

The \xmdsTag{output} element is just a container for the other tags
that specify what is to be output.  The two tags that are contained
within the \xmdsTag{output} element are: \xmdsTag{filename} and
\xmdsTag{group}.  

The \xmdsTag{filename} tag (fairly obviously) specifies the filename of
the output data file.  This tag is optional and defaults to the
simulation name (i.e. the value of \xmdsTag{name} directly within the
\xmdsTag{simulation} tag) with the string \ttt{.xsil} appended.  For
example, if we didn't specify a filename for the simulation we've
created here, then the filename \xmds would use would be
\ttt{basicSim.xsil}.  The output data file is in the XSIL
format~\cite{web:caltech_xsil} which is a handy interchange format also
using XML.

The \xmdsTag{group} tag contains a description (and to a degree the
definition) of the moments of the output data, which can be such
things as the power density of the field(s), or just means and standard
errors of the field(s).  More than one output group can be
specified, but at least one must be given.  Post-processing may be
performed before an output group is written to file, allowing one to
do some complex tasks on the data sampled in the running of the
simulation (this is put after the \xmdsTag{sampling} tag which is
coming up), however this is more involved than the discussion here, so
we won't be using it.  A good place to look if you're at all
interested is the examples directory in the distribution or see the
script repository on the \xmds web page:
\htmladdnormallink{http://www.xmds.org}{http://www.xmds.org}.

After looking for the \xmdsTag{group} tag, \xmds then expects to see a
\xmdsTag{sampling} tag within that, which defines how the group is to
be sampled.  This is also just a container for more specific tags.
Just as an aside: although it may seem a pain at this stage to have so
many containers for other containers and tags and so on, this gives
the entire document a nice structure where one builds up a simulation
from nice bite-sized (pun not intended) chunks.  Also, one really
doesn't go to the pain of writing a simulation from scratch very often
so this shouldn't be a big issue when you finally get to writing other
and more complex (and more interesting!) simulations.  Within the
\xmdsTag{sampling} tag \xmds expects to see a \xmdsTag{fourier\_space}
tag for each transverse dimension in the simulation, which tells \xmds
whether or not to Fourier transform the dimension before being sampled
and written out to disk.  For our example here, we don't have any
transverse dimensions so we won't use this tag.

%% We next tell \xmds how finely we want it to sample the output field
%% when it writes to disk.  It is possible that we could dump all of the
%% data used in the simulation to disk, but this could be \emph{very}
%% large and we probably won't get much extra information out of it.  So,
%% it is much better to take a coarser-grained subset of the data
%% actually used to do the simulation when generating output, thereby
%% producing smaller data files and making your life easier when it comes
%% to graphically interpreting the data.  To tell \xmds how finely (or
%% coarsely depending upon your point of view) to sample the data we use
%% the \xmdsTag{lattice} tag.  Exactly what this number should be depends
%% a lot upon the system being investigated, and hence it's quite a dark
%% art (or just trial and error) to get the right kind of output.  So,
%% for the purposes of our example we'll just take 50 points across the
%% $x$ dimension.  The code to do this is simply
%% {xmdsCode}
%%     <lattice> 50 </lattice>
%% {xmdsCode}
%% In general, this code is an integer specifying the sampling for each
%% of the transverse dimensions, but since there's only one here, this
%% makes the code a lot simpler.  Note that if there aren't any
%% transverse dimensions, \xmds won't search for either of
%% \xmdsTag{fourier\_space} or \xmdsTag{lattice}.

Now we need to tell \xmds the names of output moments we want to
calculate, and the C code it should use to do so.  We do this using
the \xmdsTag{moments} tag and a \ttt{CDATA} block.  For this simulation
things are quite simple since we just want the amplitude of the variables
as they evolve.  Just to make this really obvious we'll define the
output moments to be new variables called \ttt{xOut}, \ttt{yOut} and \ttt{zOut} which are just
equal to the variables \ttt{x}, \ttt{y} and \ttt{z}, but it makes it more obvious in the code
what information we are grabbing out of \xmds.  It is possible to
define as many moments as you wish, but you must define at least one.
The XML code we use is
\begin{xmdsCode}
<sampling>
  <moments> xOut yOut zOut </moments>
  <![CDATA[
      xOut = x;
      yOut = y;
      zOut = z;
  ]]>
</sampling>
\end{xmdsCode}
giving an output section which looks like
\begin{xmdsCode}
<output>
  <sampling>
    <moments> xOut yOut zOut </moments>
    <![CDATA[
	xOut = x;
	yOut = y;
	zOut = z;
    ]]>
  </sampling> 
</output> 
\end{xmdsCode}
and an overall simulation which is (finally):
\begin{xmdsCode}
<?xml version="1.0"?>
<simulation>

  <name>lorenz</name>
  <author>Paul Cochrane</author>
  <description>
    Lorenz attractor example simulation.  Adapted from the example
    in "Numerical methods for physics" by Alejandro L. Garcia, 
    page 78 (1st ed.).
  </description>

  <prop_dim>t</prop_dim>
  
  <globals>
  <![CDATA[
      const double sigma = 10.0;
      const double b = 8.0/3.0;
      const double r = 28.0;
      const double xo = 1.0;     // initial conditions
      const double yo = 1.0;     // initial conditions
      const double zo = 20.0;    // initial conditions
  ]]>
  </globals>

  <field>
    <name>main</name>
    <samples>1</samples>
    
    <vector>
      <name> main </name>
      <type> double </type>
      <components> x y z </components>
      <![CDATA[
	  x = xo;
	  y = yo;
	  z = zo;
      ]]>
    </vector>
  </field>

  <sequence>
    <integrate>
      <algorithm>RK4EX</algorithm>
      <interval> 10 </interval>
      <lattice> 10000 </lattice>
      <samples> 200 </samples>
      <![CDATA[
          dx_dt = sigma*(y - x);
          dy_dt = r*x - y - x*z;
	  dz_dt = x*y - b*z;
      ]]>
    </integrate>
  </sequence>

  <output>
    <sampling>
      <moments> xOut yOut zOut </moments>
      <![CDATA[
	  xOut = x;
	  yOut = y;
	  zOut = z;
      ]]>
    </sampling> 
  </output>
</simulation>
\end{xmdsCode}

Here is a link to the finished (gzipped) script file
\htmladdnormallink{lorenz.xmds.gz}{http://www.xmds.org/examples/lorenz.xmds.gz}
on the \xmds web site (\htmladdnormallink{http://www.xmds.org}{http://www.xmds.org}).

\subsection{Making the simulation and getting results}
\label{sec:runningXmds}

Now that the simulation script is ready, it is just a matter of
getting \xmds to generate the C++ source code and compiling that with
your system's C++ compiler.  This is a very simple process, and in the
vast majority of cases, all one has to do is enter the following
at the command prompt:
\begin{shellCode}
% xmds lorenz.xmds
\end{shellCode}

A file called \ttt{lorenz} should now appear in the same directory
as your simulation script; this file is the simulation binary
executable file.  To run the simulation merely execute the binary file
by entering its name at the command line.  You should see something
like this
\begin{shellCode}
% lorenz
Beginning full step integration ...
Sampled field (for moment group #1) at t        = 0.000000e+00
Sampled field (for moment group #1) at t        = 5.000000e-02
<snip>
Sampled field (for moment group #1) at t        = 9.950000e+00
Sampled field (for moment group #1) at t        = 1.000000e+01
maximum step error in moment group 1 was 1.248578e-05
\end{shellCode}

Once the program has finished running, you should find in the same
directory as the binary executable a file called \ttt{lorenz.xsil}.
This is the file containing your output data in a handy XML based
format that can be used to interchange data between various other
formats.  We'll look at the output here in two programs, namely Matlab (or Octave)
and Scilab.  Matlab is a commercial numerical programming language and
environment which has very powerful graphics capabilities and is used
in the scientific community extensively.  Octave is a free program which is
highly compatible with Matlab.
Scilab is very similar to
Matlab, however it is free to download and install, but doesn't have
quite the same quality as that made by Matlab.  Nevertheless, Scilab
is free, and is a handy alternative if your budget can't stretch to
Matlab.  There are subtle differences between Matlab (or Octave) and Scilab and
this is why we discuss the two here.  XSIL files can also easily be translated into 
scripts suited for input into Mathematica, gnuplot, or R using the bundled software.

Before we can start using Matlab or Scilab, we must convert the data
contained in the \ttt{.xsil} file into something that Matlab, Octave or Scilab
can understand.  To do this we use the utility program bundled with
\xmds called \ttt{xsil2graphics}.  

To generate an input file for Matlab or Octave use either
\begin{shellCode}
% xsil2graphics lorenz.xsil
\end{shellCode}
or
\begin{shellCode}
% xsil2graphics -matlab lorenz.xsil
\end{shellCode}
however the second example is redundant as a Matlab or Octave\ttt{.m} file is
the default output from \ttt{xsil2graphics}.  You should see in the
current directory a file called \ttt{lorenz.m} and a data file for
the one moment group that we sampled for \ttt{lorenz1.dat}.  For Scilab use
\begin{shellCode}
% xsil2graphics -scilab lorenz.xsil
\end{shellCode}
giving the files \ttt{lorenz.sci} and \ttt{lorenz1.dat}.
Ok, now we're ready to fire up our relevant numerical processing and
graphical environment and visualise the results.

\subsubsection{Matlab and Octave}

Start Matlab or Octave, and once at the command prompt load the
information contained in the data file by using the command
\begin{matlabCode}
>> lorenz
\end{matlabCode}
doing a \ttt{whos} should give you something similar to this
\begin{matlabCode}
>> whos
  Name               Size                   Bytes  Class

  error_xOut_1       1x201                   1608  double array
  error_yOut_1       1x201                   1608  double array
  error_zOut_1       1x201                   1608  double array
  t_1                1x201                   1608  double array
  xOut_1             1x201                   1608  double array
  yOut_1             1x201                   1608  double array
  zOut_1             1x201                   1608  double array

Grand total is 1407 elements using 11256 bytes

\end{matlabCode}
We can see that we've loaded our output data from the simulation into
the variables \ttt{xOut\_1}, \ttt{yOut\_1} and \ttt{zOut\_1} in Matlab or Octave.
The reason why the \ttt{\_1} is appended to the variable names is so
that if one defines a variable in two different moment groups, but of
the same name, then data isn't lost.  The number refers to the label
of the moment group in the simulation.  At present you don't need to
worry about these details, but just realise that they are there for
when you write more complex scripts in the future.  The error
variables seen in the \ttt{whos} listing are the differences between
the full-step integration and the half-step integration.  The
half-step integration is used for error checking purposes, so that you
can check if your simulation is likely to be giving you reliable
answers.  Now plot the data by going
\begin{matlabCode}
>> plot3(xOut_1, yOut_1, zOut_1)
>> xlabel('x')
>> ylabel('y')
>> zlabel('z')
\end{matlabCode}
and you should see a figure similar to \fig{fig:lorenzMatlabPlot}.
\begin{figure}[!h]
  \centerline{\includegraphics[width=\figwidth]{figures/lorenzMatlabPlot}}
  \caption{Three dimensional plot in Matlab of the trajectories of a
  Lorenz attractor.  Parameters used were: $\sigma = 10$, $b = 8/3$,
  $r=28$, with initial conditions of $x_0 = 1.0$, $y_0 = 1.0$, and
  $z_0 = 20.0$.}
  \label{fig:lorenzMatlabPlot}
\end{figure}

\subsubsection{Scilab}

A very similar process is necessary for viewing the results in
Scilab.  Start up scilab, and at its command prompt run the command
\begin{scilabCode}
-->exec('lorenz.sci')
 
-->temp_d1 = zeros(1,201);
 
-->t_1 = zeros(1,201);
 
-->xOut_1 = zeros(1,201);
 
-->yOut_1 = zeros(1,201);
 
-->zOut_1 = zeros(1,201);
 
-->error_xOut_1 = zeros(1,201);
 
-->error_yOut_1 = zeros(1,201);
 
-->error_zOut_1 = zeros(1,201);
 
 
-->lorenz1 = fscanfMat('lorenz1.dat');
Error Info buffer is too small (too many columns in your file ?)
 
-->temp_d1(:) = lorenz1(:,1);
 
-->xOut_1(:) = lorenz1(:,2);
 
-->yOut_1(:) = lorenz1(:,3);
 
-->zOut_1(:) = lorenz1(:,4);
 
-->error_xOut_1(:) = lorenz1(:,5);
 
-->error_yOut_1(:) = lorenz1(:,6);
 
-->error_zOut_1(:) = lorenz1(:,7);
 
-->t_1(:) = temp_d1(:);
 
 
-->clear lorenz1 temp_d1
\end{scilabCode}
to load the data into Scilab (you can safely ignore the warning), and
then to obtain a graphical output of the data, run the command
\begin{scilabCode}
-->param3d(xOut_1, yOut_1, zOut_1)    
\end{scilabCode}
which should give something along the lines of that in
\fig{fig:lorenzScilabPlot}
\begin{figure}[!h]
  \centerline{\includegraphics[width=\figwidth]{figures/lorenzScilabPlot}}
  \caption{Three dimensional plot in Scilab of the trajectories of a
  Lorenz attractor.  Parameters used were: $\sigma = 10$, $b = 8/3$,
  $r=28$, with initial conditions of $x_0 = 1.0$, $y_0 = 1.0$, and
  $z_0 = 20.0$.}
  \label{fig:lorenzScilabPlot}
\end{figure}

As we can see from both of these figures that we get the usual strange
attractor ``butterfly'' shape.  Now that we've spent a lot of time
going over the very basics of writing a simulation from scratch, we
now speed up a bit, and introduce some new \xmds tags, but still with
the theme that we are writing this all from a clean slate.  Hopefully
you will be able to see the other more powerful abilities of \xmds and
be able to start writing your own simulations.

\section{A more complex simulation}
\label{sec:moreComplexSimulation}

In this section I'll introduce how to use Fourier space in your
simulations (and the extra tags required), and explain why it is
sometimes easier to perform part of the calculation in Fourier space
and then transform back to position space.  To illustrate these
extensions to what we already know from \Sec{sec:basicSim} we'll look
at solving the one-dimensional diffusion equation
\begin{equation}
  \frac{\del a(x,t)}{\del t} = \kappa \frac{\del^2 a(x,t)}{\del x^2}
  \label{eq:diffusionEquation}
\end{equation}
where $a(x,t)$ is the field to be evolved by the differential equation
and is a function of time, $t$, and space $x$, and $\kappa$ is a
constant describing how quickly the solution diffuses.  An example of
an application of the diffusion equation is for modelling the
diffusion of (some initial distribution of) temperature in a metal
rod; over time the temperature distribution will flow from areas of
higher temperature to areas of lower temperature, eventually achieving
a uniform distribution over the entire rod.  

So why do we use Fourier space when solving this differential
equation?  The main reason is that it's a lot easier to calculate some
of the differentials in Fourier space than it is in position space.
It turns out that if one transforms position space into its respective
Fourier domain, that a partial derivative with respect to position,
just becomes $i$ (i.e.~$\sqrt{-1}$) times the coordinate in Fourier
space.  For our example, such a mapping would be:
\begin{equation}
\frac{\del}{\del x} \mapsto i k_x.
\end{equation}
Hence, in Fourier space, the second derivative on the right hand side
of \eqn{eq:diffusionEquation} is just $-k_x^2$.  Given that (discrete)
Fourier transforms aren't that hard to do on a computer (and in \xmds
we use the Fastest Fourier Transforms in the
West~\htmladdnormallink{http://www.fftw.org}{http://www.fftw.org}, so
they're pretty fast) using such a transformation improves the
calculation somewhat.

\subsection{Specifying the problem}

We now need to specify the problem properly.  To do this we must
specify an initial condition for the solution we wish to evolve, and
we must specify the boundary conditions of the domain over which we
wish to solve this particular problem.  Boundary conditions are
necessary here since we have a transverse dimension (i.e.~$x$) in this
system (recall in \Sec{sec:basicSim} we had no transverse dimensions,
only the propagation dimension of time).  In following
Garcia~\cite{Garcia:1994:1}, as we do here again, we are trying to
model the temperature diffusion of an initial temperature distribution
in a one-dimensional rod, the ends of which are kept at a constant
temperature of $T=0$.  Unfortunately, this implies Dirichlet boundary
conditions, and \xmds only implements periodic boundary conditions.
This isn't strictly true, as one \emph{can} implement absorbing
boundary conditions in \xmds (and people do in practice), however, one
has to jump through some hoops that we don't really want to bother
with here.  So, to imitate Dirichlet boundary conditions, we'll not
let the solution evolve outside the domain of the transverse
dimension, $x$.  All this means is that we make that particular domain
rather larger than was necessary, and we make sure we don't evolve it
in time for too long.  We also start with a very narrow initial
condition, and not have the diffusion coefficient too high, so as to
inhibit the diffusion a bit.  Please note that this is an example
simulation to get people used to using the syntax of \xmds and not
necessarily to pedantically solve certain physical problems, we
merely use the physical situations here to illustrate \xmds, not the
other way around.

An important solution of the diffusion equation is a Gaussian of the
form~\cite{Garcia:1994:1}
\begin{equation}
  T(x,t) = \frac{1}{\sigma(t) \sqrt{2\pi}} \exp\left[\frac{-(x-x_0)^2}{2
      \sigma^2(t)}\right],
  \label{eq:diffusionSolution}
\end{equation}
where $x_0$ is the location of the maximum and the standard deviation,
$\sigma(t)$, increases with time as
\begin{equation}
\sigma(t) = \sqrt{2 \kappa t}.
\end{equation}
This solution is also a handy initial condition, and this is the
analytical form of the initial condition we will be giving to \xmds.

\subsection{Starting off the simulation code}

Now with the boundary conditions, and the initial condition specified
analytically we are now in a position to start writing some \xmds
code.  As per usual, we start with the \ttt{<?xml version="1.0"?>}
and \xmdsTag{simulation} tags.  Then add the simulation name, which
we'll call \ttt{diffusion}, and we'll put the author name and a brief
description of what the simulation is supposed to do.  The global
variables we have in the problem we are solving are the diffusion
coefficient $\kappa$, the standard deviation of the initial Gaussian
distribution $\sigma$, and the positon of the mean of the Gaussian
distribution $x_0$.  These variables we'll call respectively,
\ttt{kappa}, \ttt{sigma}, and \ttt{x0}.  The code for this looks like:
\begin{xmdsCode}
<?xml version="1.0"?>
<simulation>

  <!-- Global system parameters and functionality -->
  <name> diffusion </name>
  <author> Paul Cochrane </author>
  <description>
    Solves the one-dimensional diffusion equation for an initial
    Gaussian pulse.  Adapted from A. L. Garcia, "Numerical Methods
    in Physics" (1994).
  </description>

  <!-- Global variables for the simulation -->
  <globals>
  <![CDATA[
      const double kappa = 0.1;
      const double sigma = 0.1;
      const double x0 = 0.0;
  ]]>
  </globals>

  <!-- more xmds code to come -->

</simulation>
\end{xmdsCode}
which we put into a file called \ttt{diffusion.xmds}.

\subsection{Describing the field}

Again, the next thing to tell \xmds about is the \xmdsTag{field}
element.  This time we have a transverse dimension which is in the $x$
direction, and we mention this using the \xmdsTag{dimensions} tag like
so:
\begin{xmdsCode}
<dimensions> x </dimensions>
\end{xmdsCode}
In the general case, \xmds expects a space-separated list of
transverse dimensions here, but in our example things are a bit
simpler.  

We now have to tell \xmds the number of grid points of the lattice in
this dimension, and over what domain in this dimension the grid is
defined.  To do these two things we use the \xmdsTag{lattice} and
\xmdsTag{domains} tags.  In our simulation here, we want to sample
from $x = -1$ to $x = 1$, and use 100 points.  Therefore, we set
\xmdsTag{lattice} to \ttt{100}, and \xmdsTag{domains} to \ttt{(-1,1)}.
Notice that the \xmdsTag{domains} tag is defined by using an ordered
pair syntax.  \xmds expects to see the domain for each transverse
dimension defined as an ordered pair i.e.~two comma-separated values
enclosed in parentheses; if there is more than one transverse
dimension, then the domains for each are defined by a list of
space-separated ordered pairs.  The last thing we need to mention
before discussing the \xmdsTag{vector} tag is that the number of
samples is set to the same value as that in \Sec{sec:basicSim},
i.e.~1, and that the name of the field, as per usual, is \ttt{main}.
The \xmdsTag{field} element now looks like
\begin{xmdsCode}
<field>
  <name> main </name>
  <dimensions> x </dimensions>
  <lattice> 100 </lattice>
  <domains> (-1,1) </domains>
  <samples> 1 </samples>

  <!-- more xmds code to come -->

</field>
\end{xmdsCode}

We now need to specify the \xmdsTag{vector} element.  This is much the
same as previously discussed in \Sec{sec:basicSim}, but with some
changes and one addition: the \xmdsTag{fourier\_space} tag.
The vector \xmdsTag{name} is again \ttt{main}, the \xmdsTag{type} this
time is \ttt{complex} though.  The reason this might be confusing is
because the temperature is a real quantity, and therefore, those who
have been reading carefully may question why \ttt{float} wasn't chosen
as the type instead.  A type of \ttt{complex} is chosen because we are
using a Fourier transform technique whereby the solution is
transformed into Fourier space to calculate the second partial
derivative in $x$ and then transformed back again, and to be able to
transform into Fourier space, we need our variables to be of complex
type.  The \xmdsTag{components} tag is set to \ttt{T}, since this is the name
of the variable about to be defined in the \ttt{CDATA} block to come,
and we tell \xmds that this component is \emph{not} defined in Fourier
space by setting the \xmdsTag{fourier\_space} tag to \ttt{no}.

We next define the \ttt{CDATA} block.  This is our C++ language
representation of \eqn{eq:diffusionSolution}, which we are using as
our initial condition.  The \ttt{CDATA} block is
\begin{xmdsCode}
<![CDATA[
    T = rcomplex(
          exp(-(x - x0)*(x - x0)/(2.0*sigma*sigma))
	     /(sigma*sqrt(2.0*M_PI))
	         ,0.0);
]]>
\end{xmdsCode}
This may look quite complicated, and possibly because we've tried to
split the code over several lines in an attempt to break up the
various parts and because this is intended for a fixed width page.
The code does nevertheless introduce some important concepts.  These
are: the \ttt{rcomplex()} function, the ability to split lines of code
over multiple lines if necessary, and the \ttt{M\_PI} variable

The \ttt{rcomplex()} function
is one of a set of utility functions added as part of \xmds to allow
users to define complex variables.  The syntax of \ttt{rcomplex()} is
\ttt{rcomplex(x,y)} where \ttt{x} and \ttt{y} are real variables
representing the real and imaginary parts of the complex number
respectively.  This explains one line of the \ttt{CDATA}
block reads merely: \ttt{,0.0);}.  The lonely comma is just separating
the two arguments to \ttt{rcomplex()}, the \ttt{0.0} is the
imaginary part of the variable we are defining, which is real, but of
\ttt{complex} type, and the closing parenthesis and semicolon just
finish off the function call syntax.  

Notice that the variable \ttt{T} is assigned to a quantity which is
defined over multiple lines.  Although this may seem strange to some,
for those familiar with C language rules will know that all the C
compiler is looking for is the semicolon as the character denoting the
end of the expression.  Therefore, it is possible to split equations
over many lines to break up complicated expressions, or to highlight
certain parts of the expression that may be important.

The \ttt{M\_PI} variable is an automatically set variable by \xmds
that is the value of $\pi$, i.e.~3.14159\ldots

Note also that one can use any of the standard mathematical functions
defined in C/C++, such as \ttt{sqrt()}, \ttt{log()} etc.

This completes the discussion of the \xmdsTag{field} element, which is
now:
\begin{xmdsCode}
<!-- Field to be integrated -->
<field>
  <name> main </name>
  <dimensions> x </dimensions>
  <lattice> 100 </lattice>
  <domains> (-1,1) </domains>
  <samples> 1 </samples>

  <vector>
    <name> main </name>
    <type> complex </type>
    <components> T </components>
    <fourier_space> no </fourier_space>
    <![CDATA[
	T = rcomplex(
              exp(-(x - x0)*(x - x0)/(2.0*sigma*sigma))
	        /(sigma*sqrt(2.0*M_PI))
	             ,0.0);
    ]]>
  </vector>
</field>
\end{xmdsCode}

\subsection{The sequence and integrate elements}

Using what we know from \Sec{sec:basicSim}, we now use
\xmdsTag{sequence} and \xmdsTag{integrate} elements to describe the
guts of what \xmds has to do.  We'll use here the \ttt{RK4EX}
algorithm, an interval of length 1, a lattice of 1000, and take 50
samples along the propagation direction.  

Now, because we're evolving the solution partially in Fourier space,
we need to define the operators that are going to be performing the
evolution.  This is done with the \xmdsTag{k\_operators} element.  The
reason why we call this the \xmdsTag{k\_operators} element is because
Fourier space is often referred to as $k$-space, and position space as
$x$-space, hence these operators are operating in $k$-space, and so
they are $k$-operators.  We next tell \xmds that the $k$-operators
(there is actually only one here) are constant over the course of the
simulation by setting the \xmdsTag{constant} tag to \ttt{yes}.  We do
this because if \xmds has to assume that the $k$-operators
\emph{aren't} constant then it has to use much slower code to evolve
the solution, and we are fortunate that for our simulation here the
$k$-operators are constant since they can be calculated via the (also
constant) \ttt{x} variable.  \xmds needs to know the name of the
operator we are going to use for our $k$-operator, and this we set
with the \xmdsTag{operator\_names} tag to be \ttt{L}.  In general, the
\xmdsTag{operator\_names} tag expects a space-separated list of the
operator names you wish to define.  We then give the C++ code
necessary to define the operator in a \ttt{CDATA} block, and for our
simulation this is
\begin{xmdsCode}
<![CDATA[
    L = -kappa*kx*kx;
]]>
\end{xmdsCode}
Note that we have one variable here that we haven't defined before:
\ttt{kx}.  This is a variable automatically defined by \xmds when we
define $k$-operators.  If we had another variable called \ttt{y} and
defined a $k$-operator for it, then it would be called \ttt{ky}.

The last thing we need to do within this section of the script is tell
\xmds how to evolve the solution---in other words, the differential
equation!  As in \Sec{sec:basicSim} we use a \ttt{CDATA} block with a
modified C++ syntax that \xmds understands to write down the
differential equation.
\begin{xmdsCode}
<![CDATA[
    dT_dt = L[T];
]]>
\end{xmdsCode}
This completes the work necessary to integrate the solution forward,
and completes the \xmdsTag{integrate} and \xmdsTag{sequence}
elements, which are:
\begin{xmdsCode}
<!-- The sequence of integrations to perform -->
<sequence>
  <integrate>
    <algorithm>RK4EX</algorithm>
    <interval>1</interval>
    <lattice>1000</lattice>
    <samples>50</samples>
    <k_operators>
      <constant>yes</constant>
      <operator_names>L</operator_names>
      <![CDATA[
          L = -kappa*kx*kx;
      ]]>
    </k_operators>
    <![CDATA[
        dT_dt =  L[T];
    ]]>
  </integrate>
</sequence>
\end{xmdsCode}

\subsection{The output element}

The \xmdsTag{output} element is much like that discussed in
\Sec{sec:basicSim}.  We need to specify a \xmdsTag{filename} tag, just
so that we are documenting everything nicely, and then \xmdsTag{group}
and \xmdsTag{sampling} elements to describe the rest of the
information \xmds needs to properly sample the data and save it out to
file.  There are two new tags that are needed because we have a
transverse dimension to worry about; these are
\xmdsTag{fourier\_space} and \xmdsTag{lattice}.  The
\xmdsTag{fourier\_space} tag is necessary to tell \xmds if the
sampling of the output is to be performed in Fourier space.  \xmds
expects to see a space-separated list of \ttt{yes} or \ttt{no} values
for each of the transverse dimensions, and for the situation here we
don't want the output sampled in Fourier space and we only have one
transverse dimension, so we set the \xmdsTag{fourier\_space} tag to
\ttt{no}.  The \xmdsTag{lattice} tag tells \xmds how finely the
transverse dimensions should be sampled, and in general expects to see
a space separated list of values telling it how many samples in the
particular transverse dimension to take.  Here we just set
\xmdsTag{lattice} to \ttt{50}.  The moments we are interested in
sampling is of the temperature, so we specify a \xmdsTag{moments} tag
value of \ttt{temperature} and give the code to calculate this moment
in a \ttt{CDATA} block as
\begin{xmdsCode}
<![CDATA[
    temperature = T;
]]>
\end{xmdsCode}

All of this information gives the \xmdsTag{output} element to be
\begin{xmdsCode}
<!-- The output to generate -->
<output>
  <filename>diffusion.xsil</filename>
  <group>
    <sampling>
      <fourier_space> no </fourier_space>
      <lattice>       50 </lattice>
      <moments>temperature</moments>
      <![CDATA[
	  temperature = T;
      ]]>
    </sampling>
  </group>
</output>
\end{xmdsCode}

\subsection{The final script}

We now have a complete script!  And this is, in its entirety:
\begin{xmdsCode}
<?xml version="1.0"?>
<!-- Example simulation: Diffusion Equation -->

<simulation>

  <!-- Global system parameters and functionality -->
  <name>diffusion</name>
  <author> Paul Cochrane </author>
  <description>
    Solves the one-dimensional diffusion equation for an initial
    Gaussian pulse.  Adapted from A. L. Garcia, "Numerical Methods
    in Physics" (1994).
  </description>
  <prop_dim>t</prop_dim>

  <!-- Global variables for the simulation -->
  <globals>
  <![CDATA[
    const double kappa = 0.1;  // diffusion coefficient
    const double sigma = 0.1;  // std dev of initial Gaussian
    const double x0 = 0.0;     // mean position of initial Gaussian
  ]]>
  </globals>

  <!-- Field to be integrated over -->
  <field>
    <name>main</name>
    <dimensions>  x   </dimensions>
    <lattice>    100  </lattice>
    <domains>  (-1,1) </domains>

    <samples>1</samples>
    <vector>
      <name>main</name>
      <type>complex</type>
      <components>T</components>
      <fourier_space>no</fourier_space>
      <![CDATA[
        T = rcomplex(
	      exp(-(x - x0)*(x - x0)/(2.0*sigma*sigma))/
	        (sigma*sqrt(2.0*M_PI))
	             ,0.0);
      ]]>
    </vector>
  </field>

  <!-- The sequence of integrations to perform -->
  <sequence>
    <integrate>

      <algorithm>RK4EX</algorithm>
      <interval>1</interval>
      <lattice>1000</lattice>
      <samples>50</samples>
      <k_operators>
        <constant>yes</constant>
        <operator_names>L</operator_names>
        <![CDATA[
          L = -kappa*kx*kx;
        ]]>
      </k_operators>
      <![CDATA[
          dT_dt =  L[T];
      ]]>
    </integrate>
  </sequence>

  <!-- The output to generate -->
  <output>
    <filename>diffusion.xsil</filename>
    <group>
      <sampling>
        <fourier_space> no </fourier_space>
        <lattice>       50 </lattice>
        <moments>temperature</moments>
	<![CDATA[
	  temperature = T;
	]]>
      </sampling>
    </group>
  </output>

</simulation>
\end{xmdsCode}

Here is a link to the finished (gzipped) script file
\htmladdnormallink{diffusion.xmds.gz}{http://www.xmds.org/examples/diffusion.xmds.gz}
on the \xmds web site (\htmladdnormallink{http://www.xmds.org}{http://www.xmds.org}).

\subsection{Making the simulation and getting results}

Running through the sequence of events necessary to generate a
simulation binary executable file, running it and producing the
results for Matlab, Octave or Scilab, you should see a sequence of
events (and output) something like following:

Generating the simulation binary:
\begin{shellCode}
% xmds diffusion.xmds
\end{shellCode}

Running the simulation:
\begin{shellCode}
% diffusion
Making forward plan
Making backward plan
Beginning full step integration ...
Sampled field (for moment group #1) at t        = 0.000000e+00
Sampled field (for moment group #1) at t        = 2.000000e-02
<snip>
Sampled field (for moment group #1) at t        = 9.800000e-01
Sampled field (for moment group #1) at t        = 1.000000e-00
maximum step error in moment group 1 was 9.909189e-09
\end{shellCode}

The forward and backward plans are the fftw routines calculating the
necessary Fourier transforms.

Generating the Matlab or Octave output:
\begin{shellCode}
% xsil2graphics lorenz.xsil
Output file format defaulting to matlab.
Output file name defaulting to 'diffusion.m'
Proccessing xsil data container 1 ...
Writing data container 1 to file ...
\end{shellCode}

Generating the Scilab output:
\begin{shellCode}
% xsil2graphics -scilab lorenz.xsil
Output file name defaulting to 'diffusion.sci'
Proccessing xsil data container 1 ...
Writing data container 1 to file ...
\end{shellCode}

\subsubsection{Matlab}

Loading the data into Matlab or Octave:
\begin{matlabCode}
>> diffusion
\end{matlabCode}

Doing a \ttt{whos}:
\begin{matlabCode}
  Name                      Size                   Bytes  Class

  error_temperature_1      50x51              20400  double array
  t_1                       1x51                408  double array
  temperature_1            50x51              20400  double array
  x_1                       1x50                400  double array

Grand total is 5201 elements using 41608 bytes

\end{matlabCode}

Plotting the data:
\begin{matlabCode}
>> mesh(t_1, x_1, temperature_1)
>> xlabel('t')
>> ylabel('x')
>> zlabel('T')
\end{matlabCode}
you should see a figure similar to \fig{fig:lorenzMatlabPlot}.
\begin{figure}[!h]
  \centerline{\includegraphics[width=\figwidth]{figures/diffusionMatlabPlot}}
  \caption{Three dimensional plot in Matlab of the diffusion of
  Gaussian pulse according to the diffusion equation.  Parameters used
  were: $\kappa = 0.1$, $\sigma = 0.1$, $x_0=0$}
  \label{fig:diffusionMatlabPlot}
\end{figure}

\subsubsection{Scilab}

Loading the data into Scilab:
\begin{scilabCode}
-->exec('diffusion.sci')
 
-->temp_d1 = zeros(50,51);
 
-->t_1 = zeros(1,51);
 
-->temp_d2 = zeros(50,51);
 
-->x_1 = zeros(1,50);
 
-->temperature_1 = zeros(50,51);
 
-->error_temperature_1 = zeros(50,51);
 
 
-->diffusion1 = fscanfMat('diffusion1.dat');
Error Info buffer is too small (too many columns in your file ?)
 
-->temp_d1(:) = diffusion1(:,1);
 
-->temp_d2(:) = diffusion1(:,2);
 
-->temperature_1(:) = diffusion1(:,3);
 
-->error_temperature_1(:) = diffusion1(:,4);
 
-->t_1(:) = temp_d1(1,:);
 
-->x_1(:) = temp_d2(:,1);
 
 
-->clear diffusion1 temp_d1 temp_d2
\end{scilabCode}

Plotting the data:
\begin{scilabCode}
-->plot3d(x_1, t_1, temperature_1)    
\end{scilabCode}
should give something along the lines of that in
\fig{fig:lorenzScilabPlot}
\begin{figure}[!h]
  \centerline{\includegraphics[width=\figwidth]{figures/diffusionScilabPlot}}
  \caption{Three dimensional plot in Scilab of the diffusion of
  Gaussian pulse according to the diffusion equation.  Parameters used
  were: $\kappa = 0.1$, $\sigma = 0.1$, $x_0=0$}
  \label{fig:diffusionScilabPlot}
\end{figure}

Note that in both \fig{fig:diffusionMatlabPlot} and
\fig{fig:diffusionScilabPlot} we have an initial Gaussian pulse at
$t=0$, which then spreads out and loses amplitude as $t$ increases.
This is the expected evolution of the solution according to the
diffusion equation.  The data points near the back of the graph, close
to $t=1$ are unlikely to be an accurate representation of the solution
at this point, and indeed, are unlikely to be correct.  Nevertheless,
we have the expected behaviour, and have now demonstrated sufficient
\xmds tags for you, the user, to be confident to go off and write your
own simulations.  We wish you the very best of luck!