\chapter{User-Defined Dynamical Systems}\label{define_ds}
In most instances, we expect that DsTool will be installed on a network of computers and
used by a variety of researchers.  It is not required, however, that each user
on a network utilize the same executable version of DsTool.  We have provided ways in
which each user may customize certain features of DsTool to suit his or her needs.  
An example of  DsTool customization is the installation of a user-defined dynamical system.
In this chapter we will detail this process.  More advanced customization, such 
as adding computational code to DsTool or adding animation to models, is
discussed in \unix{\$DSTOOL/my\_dstool/README}\index{customization}.

\section{Preliminaries}\label{prelims}
We presume DsTool is installed according to the instructions provided in
the DsTool installation information.  This procedure includes:
\begin{enumerate}
\item Defining UNIX environmental variables \unix{DSTOOL}, \unix{MY\_DSTOOL}, and \unix{ARCH}.
%\item \unix{DSTOOL\_DATA\_DIR}, 
\item Creating a local directory for DsTool.
\item Copying the contents of \unix{\$DSTOOL/my\_dstool} from the DsTool source code into the above subdirectory. 
%\item Defining \unix{UNIX} environmental variables \unix{DSTOOL}, \unix{DSTOOL\_COLOR\_DIR},
%and \unix{DSTOOL\_DATA\_DIR}.
%\item Creating a local directory for DsTool.
%\item Creating a subdirectory for user-defined models.
%\item Copying certain files from the DsTool source code into the above subdirectory.
\end{enumerate}
%The last two items are completed automatically by running the shellscript 
%\unix{install\_dsuser}.
Be certain that these steps are complete before proceeding.
%In addition, the following UNIX environmental variables are optional:
%\unix{DSTOOL\_TCL\_DIR} -- for changing the Tcl library,
%\unix{MY\_DSTOOL\_TCL\_DIR} -- for alternate Tcl interfaces,
%\unix{MY\_DSTOOL\_OOGL\_DIR} -- to hold Oogl files for custom versions,
%\unix{DSTOOL\_OOGL\_DIR} -- to hold Oogl files,
%\unix{DSTOOL\_HELP\_DIR} -- help files,
%\unix{DSTOOL\_DATA\_DIR} -- for keeping data, 
%\unix{DSTOOL\_PS\_PROLOG} -- for PostScript Prolog files.
%If any of these aren't defined, defaults will be provided.

\section{Installing a New Dynamical System in DsTool}\label{def_ds}\index{defining dynamical system}

The addition of a new dynamical system in DsTool is a two-step process.  The first
step entails writing a few procedures which define the set of governing equations for
the dynamical system (be it a vector field or a mapping) and the initial settings of
variables and parameters.  If desired, additional procedures may be written which
define derivatives (with respect to space, time, and parameters) and define an
arbitrary number of auxiliary scalar-valued functions.  The second step in the process
is to install the procedures into the libraries used to construct the executable
version of DsTool.  To help you complete the necessary steps, we provide the following checklist:
\begin{enumerate}
\item Define the dynamical system
	\begin{enumerate}
	\item Equations of motion
	\item Derivatives
	\item Inverse
	\item Auxillary equations
	\item Names and default ranges for variables
	\item Names and default ranges for parameters
	\item Names and default ranges for auxiliary functions
	\item Periodic variables
	\item Names of user-defined functions
	\end{enumerate}
\item Install the dynamical system
	\begin{enumerate}
	\item Initialization routine
	\item Title
	\item Compile source code
	\end{enumerate}
\end{enumerate}

The discussion which follows contains details of this process and a description of the
files and variables involved.  If you are not already there, change directories to
your local DsTool directory.


\subsection{Defining the Equations of Motion}\label{def_eqns}\index{defining dynamical system, equations} 
We shall use as an example a two-dimensional mapping $f=f_{\alpha,\gamma}$
defined by
\be\label{bball}\begin{array}{rcl}
  \phi_{j+1} &=& \phi_j + v_j, \\
  v_{j+1} &=& \alpha v_j - \gamma \cos(\phi_j + v_j).
\end{array}
\ee
These equations are a model for the dynamics of a ball bouncing\index{bouncing ball} on a sinusoidally
vibrating table.  The variable $\phi$ is proportional to the time, but because of the sinusoidal motion
of the table, we may think of $\phi$ as belonging to the circle parametrized by $[0,2 \pi)$.
The variable $v$ is proportional to the ball's velocity immediately after contact with the
table.  The parameter $\alpha$, $0 < \alpha \leq 1$ is the coefficient
of restitution for the ball; the parameter $\gamma > 0$ is the amplitude of the
table oscillations.  Essentially, the equations say that given the time of the $j$th impact 
($\phi_j$) and the ball's velocity at that time ($v_j$),  we know the time that the ball will next strike the
table ($\phi_{j+1}$), as well as the ball's velocity immediately after its next bounce ($v_{j+1}$).  
See~\cite{guckenheimerholmes} for information about the
dynamics of this system.

Suppose that a user wishes to study Equation~\ref{bball} numerically using DsTool.
The first task is to define the equations.  
First, copy the file \unix{GENERIC.c} to the file that you want
to edit.  For this example, call the target file \unix{bball\_def.c}.
You can use any text editor to edit \unix{bball\_def.c}.
The beginning of the file looks like:
{\csize \begin{verbatim}
#include <model_headers.h>

/* ------------------------------------------------------------------------
   function used to define the vector field or map
   ------------------------------------------------------------------------ */
int user_ds_func(f,x,p)
double *f,*x,*p;
{
}
\end{verbatim} }

First change the name of this function from user\_ds\_func() to
bball().  Then edit the function so that it properly defines the
mapping.  When you are finished, the function should look like:
{\csize \begin{verbatim}
#include <model_headers.h>

/* ------------------------------------------------------------------------
   function used to define the map
   ------------------------------------------------------------------------ */
int bball(f,x,p)
double *f,*x,*p;
{        
   f[0] = x[0] + x[1];
   f[1] = p[0] * x[1] - p[1] * cos(x[0] + x[1]);
}
\end{verbatim} }

The mapping is now defined.  We remind  novice C-programmers that arrays in C are indexed
from 0.  When bball() is called, the calling
routine is expected to pass in arrays \unix{f}, \unix{x}, and \unix{p}.  The input variables are
the current state of variables (\unix{x}) and the current parameters (\unix{p}):
\bea
\unix{x} &=& \{\unix{x[0]}, \unix{x[1]}\} = \{\phi_j, v_j\},\\
\unix{p} &=& \{\unix{p[0]}, \unix{p[1]}\} = \{\alpha, \gamma\}.
\eea
The function bball() places the new state in the array \unix{f}:
\bea \unix{f} &=& \{ \unix{f[0]}, \unix{f[1]} \} = \{\phi_{j+1}, v_{j+1}\}.
\eea

This  function\index{defining dynamical system, required}  {\em must} be defined in order 
to begin exploration of dynamics.  If the user does not wish to
input information about the system's derivative or inverse, and if
no auxiliary functions are desired, the user could proceed to Section~\ref{ic}.
For the purpose of illustration, however, we encourage the reader to 
continue reading. 

We remark at this point that although this system is a two-dimensional mapping,
the value $x[2]$ contains the current ``time.''\index{time}  In general, if there are 
$k$ dependent variables, then \unix{x[0]}, $\ldots$, \unix{x[k-1]} contain these variables and
\unix{x[k]} contains the independent variable.\index{independent variable}  
This is true both for maps and for vector fields. 

\subsection{Defining Derivative Information}\label{dfdx}\index{defining dynamical system, derivatives} 
If available, DsTool can use analytic information in
certain computational routines.  
Jacobian information is used in root-finding algorithms 
and for implicitly iterating diffeomorphisms backwards (when an explicit inverse
does not exist).  Detailed information about these algorithms is found in the DsTool Reference
Manual\index{reference manual}.

If derivative information is not provided by the
user, then DsTool will use finite difference-approximations.  This section will describe
how to write a function which provides an explicit Jacobian for Equation~\ref{bball}.
Because time is discrete for maps, it does not make sense to define a derivative with
respect to time.  
Currently, DsTool does not make use of derivatives with respect to time and parameters,
but we include them for the convenience of the user who wishes to extend the capabilities
of DsTool.

We now continue with the bouncing ball example by defining the Jacobian of $f$.
At the $j$th instant of time, the Jacobian is
\[ \left( \begin{array}{cc}
	\partial{f_1}/\partial{\phi_j} & \partial{f_1}/\partial{v_j} \\
	\partial{f_2}/\partial{\phi_j} & \partial{f_2}/\partial{v_j} \\
\end{array} \right) = \left( \begin{array}{cc}
	1 & 1 \\
\gamma \sin( \phi_j + v_j) & \alpha + \gamma \sin(\phi_j + v_j)
\end{array} \right)
\]
Find the section of the file \unix{bball\_def.c} which reads
{\csize \begin{verbatim}
/* ------------------------------------------------------------------------
   function used to define the Jacobian
   ------------------------------------------------------------------------ */
/*
int user_jac(m,x,p)
double  **m, *x, *p;
{
}
*/
\end{verbatim}}
Using a text editor, modify this code segment to read
{\csize \begin{verbatim}
/* ------------------------------------------------------------------------
   function used to define the Jacobian
   ------------------------------------------------------------------------ */
int bball_jac(m,x,p)
double  **m, *x, *p;
{
   double    temp;

   temp    = p[1] * sin( x[0] + x[1] );
   m[0][0] = 1.0;
   m[0][1] = 1.0;
   m[1][0] = temp;
   m[1][1] = p[0] + temp;
}
\end{verbatim}}

The routine which calls bball\_jac() sends in a matrix and two arrays which contain the 
current state and the current  parameters for $f_{\alpha, \gamma}$.  When bball\_jac() returns, it has filled
the matrix with the numerical Jacobian for $Df_{\alpha, \gamma}$ evaluated at the current state.
As remarked previously, writing a Jacobian routine is optional.

\subsection{Defining Information About an Inverse}\label{inverse}\index{defining dynamical system, inverse} 
In this section we show how to define an inverse or an approximate inverse for a
mapping.  We recall that the ``inverse function'' for a vector field is trivial, since
integrating the equations of motion backwards in time is the same as integrating
forward in time---except we pass in a negative time step to a numerical integrator.
Thus if we were installing a vector field, we would skip this section.  Moreover, it
frequently happens that a mapping does not have an inverse, nor is there any sort of
an approximate inverse.  In that case, do not edit the inverse routine at all.

For diffeomorphisms, it is sometimes difficult to iterate backwards in time.
The Reference Manual\index{reference manual} discusses Newton's method\index{Newton's method}
and implicitly iterating a mapping backwards\index{implicit inverse iteration}; 
it also briefly indicates some of the ways that Newton's method
can fail.  One of the most important  aspects of Newton's algorithm is its dependence on a ``good''
initial guess.  In the present section, we will discuss one way to provide a good guess.

For our bouncing ball example, it turns out that we can find an explicit inverse\index{inverse, explicit} for $f$
since $\alpha > 0$.
This inverse is the map given by
\be\label{bball_inv}\begin{array}{rcl}
  \phi_{j-1} &=& \phi_j - \frac{1}{\alpha}(\gamma \cos \phi_j + v_j), \\
  v_{j-1} &=& \frac{1}{\alpha}(\gamma \cos \phi_j + v_j).
\end{array}
\ee

To enter this inverse into DsTool, locate the section of \unix{bball\_def.c} which looks like
{\csize \begin{verbatim}
/* ------------------------------------------------------------------------
   function used to define inverse or approximate inverse
   ------------------------------------------------------------------------ */
/*
int user_inv(y,x,p)
double    *y,*x,*p;
{
}
*/
\end{verbatim}}
\noindent and modify this section to produce 
{\csize \begin{verbatim}
/* ------------------------------------------------------------------------
   function used to define inverse
   ------------------------------------------------------------------------ */
int bball_inv(y,x,p)
double    *y,*x,*p;
{
   double     temp;

   temp = ( p[1] * cos( x[0] ) + x[1] ) / p[0];
   y[0] = x[0] - temp;
   y[1] = temp;
}
\end{verbatim}}
% WAS:  f[0] = x[0] - temp;
%       f[1] = temp;
Suppose that our map did not have an explicit inverse (or we could not calculate it).
Then we would have to resort to using Newton's method to numerically solve for the inverse iterate
at each backward step.  To do so requires that we 
provide an initial guess for the inverse image of the current state.
If we provide a good guess for the inverse image, then possibly we will require only a few
Newton steps to ``sharpen'' our guess and converge to a solution (See the Reference Manual\index{reference manual}).
Suppose, for the sake of an example, that we were only interested in the bouncing ball map $f_{\alpha, \gamma}$
in the special case that $\gamma$ is small.  Then $f$ can be thought of as a perturbation
of the linear map $g$ defined by
\be\begin{array}{rcl}
  \phi_{j+1} &=& \phi_j + v_j, \\
  v_{j+1} &=& \alpha v_j,
\end{array}
\ee
and $g^{-1}$ can be written as:
\be\begin{array}{rcl}
  \phi_{j-1} &=& \phi_j - \frac{1}{\alpha} v_j, \\
  v_{j-1} &=& \frac{1}{\alpha} v_j.
\end{array}
\ee
In this case, we could code $g^{-1}$ into DsTool in place of $f^{-1}$.\index{inverse, approximate}
If we were to do this, we would modify \unix{bball\_def.c} to read:
{\csize \begin{verbatim}
/* ------------------------------------------------------------------------
   function used to define approximate inverse
   ------------------------------------------------------------------------ */
int bball_approx_inv(y,x,p)
double    *y,*x,*p;
{
   y[0] = x[0] - x[1] / p[0];
   y[1] = x[1] / p[0];
}
\end{verbatim}}
Of course, we must inform DsTool whether to interpret the procedure as a definition of the actual inverse
or merely as an approximation; Section~\ref{ic} explains how this is done.

\subsection{Defining Auxiliary Functions}\label{auxeqn}\index{defining dynamical system, auxiliary functions} 
We continue with the above example.  Suppose that there is some scalar function
of the phase space variables, time, and parameters which we wish to monitor.
The scalar function may be a physically meaningful quantity such as the 
total energy of an (almost) conservative system,
or it could be a mathematically interesting quantity such as the
largest eigenvalue of the Jacobian matrix.  Sometimes it is useful to have functions 
which represent coordinate transformations.  For example, a particular dynamical system may 
be best described in terms of spherical coordinates.   Auxiliary functions can also
be used to monitor the distance from some set of interest.  For example, if a
system spends most of its time close to the sphere $|x|=1$ then it may be interesting
to define a function $|x|$. 

For our bouncing ball example, we will monitor the quantity $v^2$.  Since $v_j$ is proportional to 
the ball's velocity just after it strikes the table for the  $j$th time, we can think 
of $v_j^2$ as being a measure of the ball's kinetic energy at time $j$.  To define this function,
find the location in the file \unix{bball\_def.c} which reads
{\csize \begin{verbatim}
/* ------------------------------------------------------------------------
   function used to define aux functions of the varbs, time, or params
   ------------------------------------------------------------------------ */
/*
int user_aux_func(f,x,p)
double *f,*x,*p;
{
}
*/
\end{verbatim}}
\noindent and edit it to produce
{\csize \begin{verbatim}
/* ------------------------------------------------------------------------
   function used to define aux functions of the varbs, time, or params
   ------------------------------------------------------------------------ */
int bball_aux(f,x,p)
double *f,*x,*p;
{
   f[0] = x[1] * x[1];
}
	\end{verbatim}}
The number of auxiliary functions is completely arbitrary; it is not necessary to have
any auxiliary functions.  When this equation is called, the calling routine passes in
the current state of the phase space variables (in \unix{x}), the current parameters (in \unix{p}),
and an array (\unix{f}) which is filled up by this function.

\subsection{Defining Labels and Initial Conditions}\label{ic}\index{defining dynamical system, initial conditions} 
This section explains how to complete the definition of a dynamical
system by giving DsTool information about the structure of phase space,
the names of variables, and initial conditions.  

This section is split into subsections.  The subsections describe how to provide 
information on   variables,  parameters,  auxiliary functions,   periodic variables,
and numerical  algorithms.

Before we begin,
find the code segment of \unix{bball\_def.c} which looks like
{\csize \begin{verbatim}
/* ------------------------------------------------------------------------
   Procedure to define default data for the dynamical system. NOTE: You
   may change the entries for each variable but DO NOT change the list of
   items.  If a variable is unused, NULL or zero the entry, as appropriate.
   ------------------------------------------------------------------------ */
int user_init()
{
/* ------------ define the dynamical system in this segment --------------- */
...
\end{verbatim}}
\noindent and change the name of the function from user\_init() to bball\_init().

\subsubsection{Variables}\label{varbs}
For our example, we have two spatial phase space variables, namely $\phi$ and $v$.
Because the equations of motion (Equation~\ref{bball}) are invariant under the
coordinate transformation $\phi \to \phi + 2n \pi$, $n \in \zed$, 
it is natural to display $\phi$ on the interval $[0, 2 \pi]$.  There is no ``natural''
interval to use in displaying the $v$ coordinate, since $v$
can be any real number, but we will choose the interval $[-30, 30]$ as a default range
on which to display $v$.
We will also need to choose a default initial condition $(\phi_0, v_0)$, which we arbitrarily
select to be the origin, $(0,0)$.

After implementing these choices, the relevant code in  bball\_init() looks like:
{\csize \begin{verbatim}
int            n_varb=2;                      /* dim of phase space           */
static char    *variable_names[]={"phi","v"}; /* list of phase varb names     */
static double  variables[]={0.,0.};           /* default varb initial values  */
static double  variable_min[]={0.,-30.};      /* default varb min for display */
static double  variable_max[]={TWOPI,30};     /* default varb max for display */
	\end{verbatim}}
We remark that \unix{TWOPI} ($2 \pi$)\index{TWOPI} and \unix{PI} ($\pi$)\index{PI} 
are two constants\index{constants} which the user may use in
defining a dynamical system. 

Although we have defined labels
and initial values for the {\em spatial} variables, the independent variable
(usually thought of as time)  is also considered to be a member
of every phase space.\index{time, member of phase space}  
The code which provides this information is given by:
{\csize \begin{verbatim}
static char    *indep_varb_name="time";  /* name of indep variable             */
static double  indep_varb_min=0.;        /* default indep varb min for display */
static double  indep_varb_max=10000.;    /* default indep varb max for display */
	\end{verbatim}}
In fact, this is the way the code looked when we copied it from \unix{GENERIC.c}, so
we do not need to make any changes to the code.  If we wanted to call the independent
variable ``iteration''  instead of ``time,'' or if we wanted to change the default plotting
range, then the code segment above would have to be appropriately modified.

\subsubsection{Parameters}
When examining bifurcations of our bouncing ball map, it is useful to graphically
depict (a subset of) the parameter space.  Since $0 < \alpha \leq 1$, we
choose the interval $[0,1]$ as the default range to display $\alpha$.  For 
$\gamma$, we choose $[0,25]$ as a default range.
As initial parameter values, we choose $(\alpha, \gamma) = (0.8, 10)$, since the dynamics
for these parameter values are fairly interesting.

To implement these choices, we edit a few more lines in bball\_init()
so that the code looks like:
{\csize \begin{verbatim}
int            n_param=2;                /* dim of parameter space          */
static char    *parameter_names[]={"alpha","gamma"}; /* list of param names */
static double  parameters[]={0.8,10.};   /* initial parameter values        */
static double  parameter_min[]={0.,0.};  /* default param min for display   */
static double  parameter_max[]={1.,25.}; /* default param max for display   */
	\end{verbatim}}

\subsubsection{Auxiliary Functions}
For our example, we have only a single function which we will call ``KE'' for kinetic energy.
From the definition of the function ($v^2$) and from the default plotting range for the
variable $v$, the interval $[0, 1000]$ is probably a suitable range on which to plot
the value of $v^2$.  Accordingly, we edit a few more lines in bball\_init():
{\csize \begin{verbatim}
int            n_funct=1;                /* number of user-defined functions */
static char    *funct_names[]={"KE"};    /* list of funct names; {""} if none*/
static double  funct_min[]={0};          /* default funct min for display    */
static double  funct_max[]={1000};       /* default funct max for display    */
	\end{verbatim}}

We remark that if we did not want to monitor any auxiliary functions then we would set 
\code{n\_funct=0}.  The array of function names, however, must contain at least an empty string
or else our code will not compile properly.  In other words, if there were no auxiliary quantities
of interest, then we could write \code{*funct\_names[]={""}}, but {\em not}
\code{*funct\_names[]={}}.

\subsubsection{Periodic Variables}
We must tell DsTool whether our phase space contains any periodic variables.  The
coding of this information is accomplished through two C-language constants:
\unix{EUCLIDEAN} or \unix{PERIODIC}.  If there are no periodic variables, then
the phase space is \unix{EUCLIDEAN}; if there is at least one periodic variable, the
phase space is classified as \unix{PERIODIC}.  For \unix{PERIODIC} phase spaces, we must
also inform DsTool as to which variables are periodic and what interval to use as
a fundamental domain.

For our example, since the equations of motion are invariant under the transformation
$\phi \to \phi + 2 n \pi$ for any integer $n$, we may consider $\phi$ to be a 
periodic variable with period $2 \pi$.  We may choose {\em any} interval of length
$2 \pi$ as a fundamental domain for the variable $\phi$.  Common choices are 
the intervals $[-\pi, \pi]$ and $[0, 2 \pi]$.  We make the latter choice.
To pass this information into DsTool, we edit yet a few more lines in bball\_init():
{\csize \begin{verbatim}
int           manifold_type=PERIODIC;   /* EUCLIDEAN or PERIODIC                */
static int    periodic_varb[]={TRUE, FALSE}; /* if PERIODIC, which varbs periodic? */
static double period_start[]={0.,0.};   /*if PERIODIC, begin fundamental domain */
static double period_end[]={TWOPI, 1.}; /*if PERIODIC, end of fundamental domain*/
	\end{verbatim}}

We remark on the variables \code{period\_start} and
\code{period\_end}.  If the $j$th coordinate is not periodic (\ie, the
value of \code{periodic\_varb}$[j]$ is \unix{FALSE}) then it does not
matter what \code{period\_start}$[j]$ and \code{period\_end}$[j]$ are because
the entries are ignored by DsTool.  Similarly, if the variable
\code{manifold\_type} is \unix{EUCLIDEAN}, then it doesn't matter what
values are given for the entries of \code{periodic\_varb}.  It is
always safe, of course, to set each entry of \code{periodic\_varb} to
\unix{FALSE}.  As mentioned in Section~\ref{varbs}, \unix{TWOPI} is a global
constant.

\subsubsection{Numerical Algorithms}
Continuing with our example, there remain a few pieces of information which 
we have not yet entered into DsTool.  We modify the last lines of code in bball\_init()
to obtain:
{\csize \begin{verbatim}
int        mapping_toggle=TRUE;         /* this is a map? TRUE or FALSE         */
int        inverse_toggle=EXPLICIT_INV; /* if so, is inverse FALSE, APPROX_INV, */
                                        /* or EXPLICIT_INV? FALSE for vec field */ 

/*  In this section, input NULL or the name of the function which contains...   */
int     (*def_name)()=bball;            /* the eqns of motion                   */
int     (*jac_name)()=bball_jac;        /* the jacobian (deriv w.r.t. space)    */
int     (*aux_func_name)()=bball_aux;   /* the auxiliary functions              */
int     (*inv_name)()=bball_inv;        /* the inverse or approx inverse        */
int     (*dfdt_name)()=NULL;            /* the deriv w.r.t time                 */
int     (*dfdparam_name)()=NULL;        /* the derivs w.r.t. parameters         */
	\end{verbatim}}
The first line means that our system is a mapping; if it were a vector field, we would
set \code{mapping\_toggle} to \unix{FALSE}.  The second line means that we have an
explicit inverse for our system.  In the event that a mapping does not have an explicit
inverse, then Newton's method will be used to calculate inverse iterates.  In this
event, we need to tell DsTool if there exists an approximate inverse (see
Section~\ref{inverse} or the Reference Manual\index{reference manual}) which can be used as an initial guess for
Newton's method.  If there is, then \code{inverse\_toggle} is set to \unix{APPROX\_INV};
otherwise, it is set to \unix{FALSE}.  We remark that if \code{mapping\_toggle} is set
to \unix{FALSE} (that is, the system is a vector field), then the value of
\code{inverse\_toggle} is ignored by DsTool.


In the last section of code, we provide DsTool with the names of any functions we defined previously.
For our example, this means filling in the names of the functions which compute the
equations of motion, the Jacobian, auxiliary functions, and the inverse.  If any of
these functions were omitted, then the user should enter \unix{NULL} instead of a
name. For example, the function which defines the inverse will be \unix{NULL} if you are
installing a vector field, and the function which defines the derivative with respect to
time should always be \unix{NULL} for mappings.

At this point, the user should save and exit the file \unix{bball\_def.c}.

\subsection{Installing a Defined Dynamical System}\label{install_sys}\index{defining dynamical system, installing}\index{installing dynamical system}
It is assumed that the reader has completed the steps in the preceding sections of this chapter.
In this section, we describe the installation of the dynamical system  defined using the
above  steps.  As an  illustration, we continue with the bouncing ball example.

Our first task is to modify the file \unix{user.c}.  There are two steps to this process: 
\begin{itemize}
\item Inform DsTool about the name of the initialization routine for the dynamical system we
desire to install.  In our bouncing ball example, this function is  bball\_init().
\item Enter this name into the DsTool list of dynamical systems, along with a title.
\end{itemize}

During the execution of DsTool, the user may select one of several installed dynamical systems.
The program maintains a list of installed models, along with the procedures which initialize
and define them.  
The user-defined models are maintained in the 
file \unix{user.c} of the user's local DsTool directory.  This file contains 
three relevant blocks of C code: the first block defines an array of categories used to collect dynamical systems
into related classes, the second consists of several lines beginning with the declaration
\unix{extern~int}, the last block is a C-language structure which contains titles for the 
models along with names of the models' initialization procedures.  The relevant pieces
of the file \unix{user.c}
delivered with DsTool contain the following code:
{\csize \begin{verbatim}
/* ----------------------------------------------------------------
 * INCLUDE USER DYNAMICAL SYSTEMS CATEGORIES HERE
 *
 * ----------------------------------------------------------------
 */
char *USER_DS_Category[] = { "User models" /* Category 0 */
                           };
 
/* ----------------------------------------------------------------
 * INCLUDE USER DYNAMICAL SYSTEM NAMES HERE
 *
 * We have put in the Lorenz system as an example
 * ----------------------------------------------------------------
 */
 
/* declare the model initialization procedure here */
extern int lorenz_init();
 
 
/* list the category, a short model name, and the initialization proc here */
struct DS_DataS USER_DS_Sel[]= {
  { 0, "Lorenz system", lorenz_init }
};
\end{verbatim}}

To tell DsTool the name of the initialization routine which defines our bouncing ball model,
we need to add the line
{\csize \begin{verbatim}
extern int  bball_init();
\end{verbatim}}
\noindent to the second block of code in  \unix{user.c}. 
To install our model, we replace the line of code
{\csize \begin{verbatim}
  { 0, "Lorenz system", lorenz_init }
\end{verbatim}}
with the line of code
{\csize \begin{verbatim}
  { 0, "Bouncing Ball", bball_init }
\end{verbatim}}
\noindent in the definition of the data structure \unix{USER\_DS\_Sel} within the third block of code.  
The variable \unix{USER\_DS\_Sel} is an array, each element of which is a 
data structure which defines a dynamical system.  The first element of the structure is 
a number which specifies the user category
to which the dynamical system belongs.  The category is
used to group together systems which share similar properties.  
These categories are defined in the first block of code.
Users cannot add dynamical systems to the standard list of categories.
However, they can create as many new categories as they wish, by modifying the category list like so:
{\csize \begin{verbatim}
char *USER_DS_Category[] = { "User models", /* Category 0 */
                             "Project with Bob", /* Category 1 */
                             "My favorite dynamical systems" /* Category 2 */
                           };
	\end{verbatim}}
This array provides the title for each category.  For our example, the category index is 0, so our bouncing ball
model belongs to the category named ``User models.''  (Recall
that arrays in C are indexed from
0 so \code{DS\_Category}$[0]$ is not an invalid array index.)  Every dynamical 
system must belong to a valid category.  

The next element of the dynamical system data structure is the title of the system being defined.
We have chosen ``Bouncing Ball'' as the title of our system.  Any descriptive title will do, provided
that the length of the title is not larger than the global constant \unix{MAX\_LEN\_DS\_TITLE}.
The value of this constant is found in the  file \unix{\$DSTOOL/src/include/constants.h}.

The last element is the name of the function which contains the initialization routine
for the new dynamical system.  For our example, the name is \unix{bball\_init}.

Note that each dynamical system data structure must be enclosed by braces and followed
by a comma.  That is, all except the last.  The last element in the array of dynamical systems
(our new dynamical system ``Bouncing Ball'' in our example code) should not have a trailing comma.

We are now ready to compile the source code we have written.
To do this, edit the \unix{Makefile} in your local DsTool directory.
Add the file \unix{bball\_def.c} to the \unix{USER\_SRCS}
list of files, and the file \unix{bball\_def.o} to the \unix{USER\_OBJS} list of files.
In the example described above, the corresponding Makefile entry would be:
{\csize \begin{verbatim}
        USER_SRCS = user.c bball_def.c
        USER_OBJS = user.o bball_def.o
	\end{verbatim}}
Now save and exit the \unix{Makefile}. 

If our model compiles without errors, we can create a custom version of DsTool
called \unix{my\_dstool}.  To do this, execute the command \unix{make} in your
local DsTool directory.  

Once this is done, you can ensure that this version of DsTool is executed by settting
the UNIX environmental variable \unix{MY\_DSTOOL} to your local DsTool directory.  
When the script \unix{dstool\_tk} is run from any directory, it will execute the binary 
\unix{\$MY\_DSTOOL/bin/\$ARCH/my\_dstool}.  A message will be printed confirming which binary
is being executed.
This version of DsTool will include the bouncing ball equations among its installed dynamical
systems.

\section{A Vector Field Example: Duffing's Equations}\label{vfield}
In this section we present an example of installing a user-defined vector field into DsTool.
It is assumed that the reader is familiar with the material in Section~\ref{def_ds}.

Our task is to install the vector field defined by the equations
\be\label{duff1}\begin{array}{rcl}
  \dot{u} &=& v, \\
  \dot{v} &=& u - u^3 - \delta v - \gamma \cos \omega t.
\end{array}
\ee
This system is known as Duffing's equations and we direct the reader to 
\cite{guckenheimerholmes} and references therein for an exposition of the dynamics of this system.

To define and install this system, change to your local DsTool directory
and copy the file \unix{GENERIC.c} to the file \unix{userduffing\_def.c}.
Use any text editor to edit \unix{userduffing\_def.c}.

The segment of code which defines the equations of motion (Equation~\ref{duff1}) 
will be called userduffing().  This function is defined as follows:
{\csize \begin{verbatim}
/* ------------------------------------------------------------------------
   function used to define the vector field 
   ------------------------------------------------------------------------ */
int userduffing(f,x,p)
double *f,*x,*p;
{
   f[0] = x[1];
   f[1] = x[0] - x[0]*x[0]*x[0] - p[0]*x[1] - p[1]*cos( p[2]*x[2] );
}
	\end{verbatim}}
\noindent Here we have defined the variables and parameters as
\bea
	\unix{x} &=& \{ \unix{x[0]}, \unix{x[1]} \} = \{ u, v \}, \\
	\unix{p} &=& \{ \unix{p[0]}, \unix{p[1]}, \unix{p[2]} \} = \{ \delta, \gamma, \omega \}.
\eea
As usual, the independent variable is stored after the spatial variables so that
\unix{x[2]} is time.  The function userduffing() returns in the array \unix{f} the strength of the
vector field with parameters \unix{p}, evaluated at the point \unix{x}.

We now define the Jacobian of Equation~\ref{duff1} by writing the function userduffing\_jac():
{\csize \begin{verbatim}
/* ------------------------------------------------------------------------
   function used to define the Jacobian
   ------------------------------------------------------------------------ */
int userduffing_jac(m,x,p)
double  **m, *x, *p;
{
   m[0][0] = 0.0;
   m[0][1] = 1.0;
   m[1][0] = 1.0 - 3.0 * x[0] * x[0];
   m[1][1] = -p[0];
}
\end{verbatim}}

Since our equations are for a vector field, we do not need to define an inverse function.
Our vector field is time-dependent, so we can define a function which returns
the temporal derivatives of the vector field.   This function is not yet used by DsTool, and
there is no template for such a function 
in \unix{GENERIC.c}, but it is easy enough to write:
{\csize \begin{verbatim}
/* ------------------------------------------------------------------------
   function used to define temporal derivatives
   ------------------------------------------------------------------------ */
int userduffing_dfdt(f,x,p)
double  *f, *x, *p;
{
   f[0] = 0.0;
   f[1] = -p[1] * p[2] * sin( p[2] * x[2] ); 
}
\end{verbatim}}

Since our vector field is Hamiltonian for $\delta = 0$, we choose 
\[ H(u,v) = \frac{v^2}{2} - \frac{u^2}{2} + \frac{u^4}{4}
\]
as an auxiliary function.  We will also define a second auxiliary function.  Since our vector field is
periodic in time (with period $2 \pi / \omega$), it is a common technique to look at the time-$2 \pi / \omega$
map in order to better understand the dynamics.  This ``stroboscopic map'' can also be thought
of as a Poincar\'e map for the extended phase space.  
We are thus interested in the times $t$  for which $\sin( wt + t_0) = 0$.  Choosing $t_0 = 0$
we define the procedure userduffing\_aux() by:
{\csize \begin{verbatim}
/* ------------------------------------------------------------------------
   function used to define aux functions of the varbs, time, or params
   ------------------------------------------------------------------------ */
int userduffing_aux(f,x,p)
double *f,*x,*p;
{
   double   temp;

   temp = 0.5 * x[0] * x[0];
   f[0] = 0.5 * x[1] * x[1] - temp + temp * temp;
   f[1] = sin( p[2] * x[2] );
}
\end{verbatim}}

We are now ready to define the labels and initial conditions for Duffing's equations
by writing the function userduffing\_init().  We choose $[-2, 2]$ as default plotting ranges for both 
$u$ and $v$.  Initializing the variables is straightforward:
{\csize \begin{verbatim}
int userduffing_init()
{
/* ------------ define the dynamical system in this segment --------------- */

int            n_varb=2;                  /* dim of phase space           */
static char    *variable_names[]={"u","v"}; /* list of phase varb names   */
static double  variables[]={0.5,0.5};     /* default varb initial values  */
static double  variable_min[]={-2.,-2.};  /* default varb min for display */
static double  variable_max[]={2.,2.};    /* default varb max for display */

static char    *indep_varb_name="time";   /* name of indep variable             */
static double  indep_varb_min=0.;         /* default indep varb min for display */
static double  indep_varb_max=1000;       /* default indep varb max for display */
\end{verbatim}}

Defining the parameter ranges is somewhat arbitrary. Sometimes it is difficult to 
tell {\em a priori} what range of parameters will provide interesting bifurcations.
This is not a big problem
since it is a trivial matter to change the range upon which a function (or parameter) is plotted
once within DsTool.  
We choose $[0,1]$ as a range for each parameter.  We choose $(\delta_0, \gamma_0, \omega_0)
= (0.25, 0.4, 1.0)$:
{\csize \begin{verbatim}
int           n_param=3;                  /* dim of parameter space        */
static char   *parameter_names[]={"delta","gamma", "w"}; /* list of param names */
static double parameters[]={0.25,0.4,1.}; /* initial parameter values      */
static double parameter_min[]={0.,0.,0.}; /* default param min for display */
static double parameter_max[]={1.,1.,1.}; /* default param max for display */
\end{verbatim}}

The initialization of our two auxiliary functions is accomplished by the code segment:
{\csize \begin{verbatim}
int            n_funct=2;               /* number of user-defined functions */
static char    *funct_names[]={"H", "P_Section"}; /* list of funct names; {""} if none */
static double  funct_min[]={-4.,-1.};   /* default funct min for display    */
static double  funct_max[]={4.,1.};     /* default funct max for display    */
\end{verbatim}}
\noindent As in the case with parameters, it is sometimes difficult to choose {\em a priori}
what appropriate ranges for the functions should be.  

The manifold type for Duffing's equations is \unix{EUCLIDEAN} since we do not have any 
periodic spatial variables.  Thus we do not need to modify the following segment of code:
{\csize \begin{verbatim}
int           manifold_type=EUCLIDEAN; /* EUCLIDEAN or PERIODIC                 */
static int    periodic_varb[]={FALSE, FALSE}; /* if PERIODIC, which varbs periodic? */
static double period_start[]={0.,0.};  /* if PERIODIC, begin fundamental domain */
static double period_end[]={1., 1.};   /* if PERIODIC, end of fundamental domain*/
\end{verbatim}}

The last segment of code we need to modify is the segment which tells DsTool which numerical algorithms
may be used on the Duffing's equations.  
We complete the definition of userduffing\_init() with the code segment:
{\csize \begin{verbatim}
int           mapping_toggle=FALSE;       /* this is a map? TRUE or FALSE         */
int           inverse_toggle=FALSE;       /* if so, is inverse FALSE, APPROX_INV, */
                                        /* or EXPLICIT_INV? FALSE for vec field */ 

/*  In this section, input NULL or the name of the function which contains...   */
int     (*def_name)()=userduffing;          /* the eqns of motion                   */
int     (*jac_name)()=userduffing_jac;      /* the jacobian (deriv w.r.t. space)    */
int     (*aux_func_name)()=userduffing_aux; /* the auxiliary functions              */
int     (*inv_name)()=NULL;             /* the inverse or approx inverse        */
int     (*dfdt_name)()=userduffing_dfdt;    /* the deriv w.r.t time                 */
int     (*dfdparam_name)()=NULL;        /* the derivs w.r.t. parameters         */
\end{verbatim}}


As in Section~\ref{def_ds}, we now need to edit two other files.
We edit \unix{Makefile} and add \unix{userduffing\_def.c}
to the \unix{USER\_SRCS} list, and add \unix{userduffing\_def.o} to the
\unix{USER\_OBJS} list.  We also edit the file
\unix{user.c} and add the lines
{\csize \begin{verbatim}
extern int  userduffing_init();
	\end{verbatim}}
\noindent and 
{\csize \begin{verbatim}
{0, "Duffing's Eqns", userduffing_init},
	\end{verbatim}}
\noindent to  the second and third blocks of code, respectively.

Typing \unix{make}
will create a local executable of DsTool which includes Duffing's equations among its installed dynamical
systems.


\subsection{A Few Remarks on the Definition of Duffing's Equations}
Recall that when we installed Duffing's equations as a time-dependent vector field,
we defined $\sin \omega t$ as an auxiliary function and claimed that we could
use it to study the time-$2 \pi/\omega$ stroboscopic map.  
In theory, there is nothing wrong with
this, however in practice we will encounter numerical errors in the evaluation of 
transcendental functions such as $\sin \omega t$ for large values of $t$.  
Since we are often interested in generating
Poincar\'e maps for extremely long times, and since the function $\cos \omega t$
also appears in the definition of our vector field, the user may want to extend phase space
by introducing the variable $\theta$.  Then we can rewrite Duffing's
equations in the form
\be\begin{array}{rcl}
  \dot{u} &=& v, \\
  \dot{v} &=& u - u^3 - \delta v - \gamma \cos \omega \theta, \\
  \dot{\theta} &=& 1,
\end{array}
\ee
where $(u,v,\theta) \in \real^2 \times S^1$, and $S^1$ is the circle of length $2 \pi/ \omega$.
That is, $\theta$ takes values in $[0, 2 \pi/ \omega)$.
The  problem with this formulation is that DsTool {\em cannot handle periodic variables
whose length depends on a parameter!}  To overcome this difficulty, we change coordinates
via the transformation $\psi = \omega \theta$.  Thus we could study Duffing's  equations
on extended phase space in the form
\be\label{duff2}\begin{array}{rcl}
  \dot{u} &=& v, \\
  \dot{v} &=& u - u^3 - \delta v - \gamma \cos \psi, \\
  \dot{\psi} &=& \omega,
\end{array}
\ee
where $(u,v,\psi) \in \real^2 \times S^1$, and $S^1$ is now the circle of length $2 \pi$.

The advantage to an extended phase space such as we have for Equation~\ref{duff2}
is that it is trivial to plot Poincar\'e sections for this set of equations because we
can make $\psi$ a periodic variable.  This allows us to request that DsTool only plot
points when $\theta=\theta_0$ for some $\theta_0$. 
In contrast, DsTool {\em never treats time as a periodic variable},
so we needed to define the auxiliary function $\sin \omega t$ in order to be able to generate a
Poincar\'e map for Duffing's equations.


\section{Deleting Dynamical Systems}\index{deleting dynamical system}
Suppose that we want to remove the dynamical system ``Bouncing Ball'' from 
section \ref{def_ds}, 
which is defined in the file \unix{bball\_def.c}.  To remove this:
\begin{enumerate}
\item Edit the \unix{Makefile} in your local DsTool directory.  Delete the name of the file
which contains the dynamical system from the \unix{USER\_SRCS}
list and the name of the corresponding object file from the \unix{USER\_OBJS} list. 
In our example, we would remove the file names \unix{bball\_def.c} and \unix{bball\_def.o} 
from the \unix{Makefile}.

\item Edit the file \unix{user.c} and
delete all references to the dynamical system you wish to remove.  To remove the Bouncing ball system,
delete  the lines
{\csize \begin{verbatim}
extern int	bball_init();
	\end{verbatim}}
and
{\csize \begin{verbatim}
{0, "Bouncing Ball", bball_init},
	\end{verbatim}}
from the file.  

\item Remake your custom executable as explained in previous sections.
\end{enumerate}


We note that you cannot delete any of the standard models that come 
with DsTool.  You can only delete models that you have previously 
added to your own custom executable.