\documentclass[a4paper]{article}
\usepackage{amsmath}
\usepackage{txfonts}
\usepackage{natbib}

\usepackage{graphicx}
%\usepackage{amssymb}
\usepackage{supertabular}
\usepackage{array}
%\usepackage[colorlinks,final]{hyperref}

\usepackage{aas_macros}

\usepackage{standalone}
\usepackage{sphinx}

\newcommand{\VERNO}{0.4.0 }
\newcommand{\var}[1]{\texttt{#1}}
\newcommand{\singlecommand}[1]{\\ \newline \indent \var{#1} \\ \newline \noindent}
\renewcommand{\vec}[1]{\ensuremath{\mathbf{#1}}}
\usepackage{listings}

\textheight 9.2in
\textwidth 6.2in
\oddsidemargin 0pt
\topmargin -40pt


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}

\lstset{language=python,keywordstyle=\color[rgb]{0,0,1},
        commentstyle=\color[rgb]{0.133,0.545,0.133},
        stringstyle=\color[rgb]{0.627,0.126,0.941},
        basicstyle=\ttfamily}

\title{GANDALF v\VERNO userguide}
\author{David Hubber, Giovanni Rosotti \& Richard Booth}

\maketitle
\tableofcontents

\newpage


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Overview of code}
GANDALF is a new SPH and N-body code written in C++ and Python.  Although partly based on some of the algorithms and code structures used in SEREN (Hubber et al. 2011), it has been written from scratch and contains many new features and optimisations which significantly improve the speed, functionality and usability of the code.  It has been written for several purposes.
\begin{itemize}
\item GANDALF has been written with a fully object-oriented philosophy in C++.  This is to improve the maintainability of the code and also to allow the code to be easily extended in the future with as little disruption to the original code-base as possible.
\item GANDALF can be run in three different modes; as a standard C++ executable run from the command line, run from a python script or run inside an interactive python environment.  Both the script and interactive python modes include a graphical output that can be used to visualise previously run simulations, or to interactively visualise simulations as they are run on-the-fly.
\item GANDALF contains a variety of SPH algorithms such as conservative 'grad-h' SPH (Springel \& Hernquist 2002, Price \& Monaghan 2004) and Saitoh \& Makino (2013).
\item The N-body component of GANDALF contains a variety of common integration schemes such as 4th/6th-order Hermite with time-symmetric integration.
\item GANDALF can generate initial conditions at run-time, as opposed to SEREN where the initial conditions had to be prepared a priori to running the simulation.  The initial conditions may either be hard-coded into the C++ code, or set-up via a python script using numpy arrays.
\item GANDALF uses additional parameters to switch on various physics options, as well as requiring them to be activated in the Makefile.  In comparison, SEREN controlled physics options exclusively in the Makefile.  Therefore, GANDALF need not be re-compiled everytime a different setting is required.
\end{itemize}

GANDALF has been developed as a fully open-source project hosted on the site github at the address \var{https://github.com/gandalfcode/gandalf}.  Since GANDALF is currently in the beta-testing phase, it would be very helpful for the authors if new users could provide feedback of any problems with the code, or with any suggestions for improvements of existing features.  Ideally, users would submit bug reports to the GANDALF github page (\var{https://github.com/gandalfcode/gandalf/issues}) since other users can see that a bug has been identified before sending duplicate bug reports.


\newpage


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Download \& Installation} \label{S:INSTALL}
The easiest way to download the code is using git. Just cd to the folder you prefer and type:\\
\newline
\noindent \var{git clone https://github.com/gandalfcode/gandalf.git} \\
A new folder called ``gandalf'' will be created in the current directory and a fresh version of the code will be downloaded inside it. The use of the software version control system git makes easy to stay up-to-date with the changes in the code while keeping your own version. If you want to know more about git, it is plenty of tutorials on-line. A good starting point is \var{http://git-scm.com/book/en/v2}.

\subsection{Requirements}

If you are not interested in the Python library, you need only a C++ compiler. In order to use all features in GANDALF, the following programs and libraries must be installed :
\begin{itemize}
\item C++ compiler
\item Python 2.7
\item swig compatible with python 2.7
\item matplotlib compatible with python 2.7
\item numpy
\item scipy
\end{itemize}

To reduce the number of needed libraries, we ship together with GANDALF the following packages:
\begin{itemize}
\item pyparsing
\item cmd2
%\item ffmpeg (for generating movies automatically)
\end{itemize}
We thank the authors for writing this software and releasing them under the MIT license.

Finally, in order to generate movies through the Python library, you need to have ffmpeg installed and invokable from the command line. You do not need it to compile the code or to use it if you are not generating movies.

\subsection{Linux}
All of these programs/libraries can be found in most standard Linux installations, and if not, will be available to download from most package managers (e.g. apt-get, yum, pacman).

\subsection{Mac OS X}
For Mac users, all programs can be installed with either fink, MacPorts or homebrew.  It is preferable that all are downloaded with just the one package manager in order to ensure they are compatible and function together correctly. Other options for installation are the Anaconda or Enthought Python distributions. A good reference which compares the available possibilities is the \href{http://python4astronomers.github.io/installation/recommended_options.html}{Python4Astronomers website}.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Command-line code compilation}

To compile GANDALF only as a C++ executable to be run from the command line without visualisation via the python program, first set the chosen C++ compiler in the Makefile (see section \ref{S:MAKEFILE} for a full list of the variables you can change) by setting the variable \var{CPP} (if you leave it blank, you will get the default c++ compiler as defined by your system) and type \\
\newline
\noindent \var{make executable} \\

\noindent The code is compiled and linked with the chosen C++ compiler and all python components are ignored.  If you do not have python installed on your system, or are having trouble getting the python components to function correctly, then the C++ executable can still be compiled and run stand-alone. Differently from other codes, in GANDALF we adopt the philosophy that the code should be compiled only once. If you need to change parameters, you can do it by setting them via the parameter files; in most cases no recompilation is needed. For this reason, if your physical set up does not require modifying the code, it could be that you need to compile GANDALF only once. If that is the case, you might want to copy the executable to one of your system directories. In this case, the command\\
\newline
\noindent \var{make install} \\

will copy the executable (provided you have already compiled the code) to your \var{/usr/local/bin} folder. In this way you should be able to invoke GANDALF from the command line just typing \var{gandalf} (i.e. without the need to specify the full path). Note that you might need to be root to write to that folder (in this case, use the sudo command).


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Python code compilation}
To compile all code elements including the python components, then complete the following steps :
\begin{itemize}

\item Install all python-related programs listed in Section \ref{S:INSTALL}

\begin{itemize}
\item For Linux machines, python is usually installed by default.  The additional python libraries should be easily installed using either the main package manager,  {\var pip} or {\var easy\_install}.
\item For Mac OS X, python is installed by default.  However, it does not have directly compatible versions of all the required libraries.  Therefore, it is required to install an additional version of python 2.7 using a 3rd party package manager like fink, Macports or homebrew.
\end{itemize}

\item Set the required version of python in your Makefile (See Section \ref{S:MAKEFILE}).  Since operating systems usually have more than one version of python installed, it is important to ensure that \var{make} uses the correct version when compiling the code.  This is set using the \var{PYTHON} variable in the Makefile.

Also, \var{make} requires the location of the python and numpy libraries.  In most cases, \var{make} will be able to locate these libraries automatically (see description of Makefile options in section \ref{S:MAKEFILE}).  However, if there is a problem, or you wish to use an alternative version of these libraries installed elsewhere on your system, then these can be set by the \var{PYLIB} and \var{NUMPY} variables in the Makefile.


\item To compile all components of the code, including the python libraries with swig, type \\
\newline
\noindent \var{make} or (equivalently) \var{make all}\\

\item You now need to tell python how to find GANDALF when you want to use it. You have two options:
\begin{itemize}
\item Add the location of the folder containing the main GANDALF directory to the \var{PYTHONPATH} environment variable.  If you are using bash or related shells, then add the line \\
\newline
\noindent \var{export PYTHONPATH=XXX/YYY:\$PYTHONPATH} \\

\noindent to your \var{.bashrc} (or \var{.bash\_profile} on a Mac) script where \var{XXX/YYY/gandalf} is the absolute path of the GANDALF directory (note that you must not include ``gandalf'' in the folder name!). If you are using csh, tcsh or related shells, then add \\
\newline
\noindent \var{setenv PYTHONPATH ``XXX/YYY:\$PYTHONPATH''} \\

\noindent to your relevant shell configuration file. Remember that you have to close and reopen the shell for this change to take effect, or to \var{source} the shell configuration file. In this case, python will use the version of GANDALF present in the folder you specified. Note that this \textit{requires} the folder to be called \var{gandalf}; if you give another name, python will not be able to find GANDALF.

\item If you know you will not be making changes to GANDALF, consider installing GANDALF in one of your system folders. In this case, just type \\
\newline
\var{make installpython}
\\

\noindent and the GANDALF python library will be installed (this requires that you have already used \var{make} to compile the code). In this case, there is no need to set the variable \var{PYTHONPATH}. Also, the folder where you keep your original code does not need to have any special name. If later you do make changes to the code, remember to re-run this command - otherwise the installed version remains the old one! Note that you might need to be root (or to use sudo) to run this command.

\end{itemize}

\end{itemize}



\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Makefile options} \label{S:MAKEFILE}
The GANDALF Makefile is used to select options which are used to compile the code with.  If the user wishes to change any compile-time options, the code must be recompiled from scratch by typing {\var make clean} and then \var{make}.  For users of SEREN, the GANDALF Makefile has been simplified with many options either not present (since various specialist algorithms have not been implemented) or have been transferred to the parameters file.  This has been done in order to make it less likely that the wrong Makefile options are used in simulations, and also to stop the need to recompile the code completely so often when using slightly different options.

Note that the Makefile in the root folder is just a convenience to set the variables you need - the actual Makefile is contained in the \var{src} folder.

\begin{itemize}

\item CPP : C++ compiler. If left blank (default), uses the default defined by your system. Depending on the values set here, we try to guess the compiler flags for your system. If you want to use a compiler different from the ones we support, you need to manually set the flags in the Makefile in the \var{src} folder. Note also that when using MPI we assume that you use g++/clang. If this is not the case, you need to manually edit the compiler flags in the Makefile in the \var{src} folder.\\
\begin{tabular}{ll}
g++ & : GCC C++ compiler \\
icpc & : Intel C++ compiler \\
clang & : CLANG compiler\\
mpic++ & : Compiles the code with MPI. \\
\end{tabular}

As an alternative, if you do not have the MPI compiler wrappers available (mpic++) you can enable MPI by adding the line
\begin{lstlisting}
MPI=1
\end{lstlisting}
to the makefile. However, be warned that you may need to specify the location of the MPI header file, mpi.h, for compilation to succeed. You can do this by adding a line of the following form.
\begin{lstlisting}
INCLUDE = -I/path/to/mpi_header_directory
\end{lstlisting}
Similarly, you may find that the default paths to the MPI library (-lmpi -lmpi++) may need to be changed in src/Makefile. \textbf{Where possible we highly recommend using the MPI compiler wrapper, mpic++.}


\item PYTHON : name of python command-line executable (e.g. python, python2.7)

\item COMPILER\_MODE : Set compiler flags for production or debug runs \\
\begin{tabular}{ll}
DEBUG & : Set all debug compiler options, including flags to use gdb debugger and full warning output (-Wall) \\
STANDARD & : Standard optimisation options (-O3) \\
FAST & : -O3 + fast flag options.  Uses 'potentially unsafe' fast maths optimisation.
\end{tabular}

\item PRECISION : Floating point precision \\
\begin{tabular}{ll}
SINGLE & : 32-bit precision floating point variables \\
DOUBLE & : 64-bit precision floating point variables
\end{tabular}

\item PYSNAP\_PRECISION: Precision when using the python library (see section \ref{S:PYTHONSCRIPT}). The values returned by python will be in either single or double precision, depending on the value selected here. Note that this has nothing to do with the precision of the snapshots written to the disc! Most of the times one wants to match the precision used in the calculation. But to reduce the memory footprint required in the analysis, one might want to read the data in single precision even if it is stored in double precision on the disc.\\
\begin{tabular}{ll}
SINGLE & : 32-bit precision floating point \\
DOUBLE & : 64-bit precision floating point
\end{tabular}

\item OPENMP : Activate OpenMP directives during compilation (0 or 1)

\item OUTPUT\_LEVEL : Amount of output produced by code \\
\begin{tabular}{ll}
0 & : No additional output \\
1 & : Minimal output of main-loop routines \\
2 & : Code routine marker output for all steps
\end{tabular}

\item DEBUG\_LEVEL : Amount of extra debug checking done by code \\
\begin{tabular}{ll}
0 & : No extra debugging computations and output \\
1 & : Activate assert debug statements in code \\
2 & : Activate extra (expensive) debugging computations in code
\end{tabular}

\item FFTW : Include FFTW (Fast Fourier transform) library for initial conditions (0 or 1).  The FFTW\_LIBRARY and FFTW\_INCLUDE variables should contain the library links and include directory if different from the standard Linux directories (otherwise leave blank).

\item GSL : Include GSL (GNU Scientific library) which is required for Ewald forces (0 or 1). The GSL\_LIBRARY and GSL\_INCLUDE variables should contain the library links and include directory if different from the standard Linux directories (otherwise leave blank).

\item PYLIB : Path to directory that includes prefered python libraries.  If left blank, {\var make} will use python routines to locate the libraries automatically.

\item NUMPY : Path to directory that includes prefered numpy libraries.  Ase with PYLIB, If left blank, {\var make} will use python routines to locate the library locations automatically.

\item GTEST : Path to directory containing Google test suite library for unit testing (Currently set witu environment variable)

\end{itemize}

\newpage


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Basic usage}
GANDALF can be run in three principle modes.
\begin{itemize}
\item Command-line mode, where the code is run as a C++ executable via the command line with a parameters file that selects which code options are used,
\item via a python script where the code is run in a python environment with various plotting options available as well as running the code,
\item in interactive mode where all the code and plotting options can be run directly by typing them into an interactive python shell.
\end{itemize}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Command-line mode}

\noindent To run a simulation with the C++ executable on the command line, type : \\
\newline
\var{./gandalf PARAMSFILE} \\
\newline
where PARAMSFILE is the name of the parameters file for that simulation.  If that parameters file does not exist in the current directory, or contains invalid parameter options, then the program will quit citing an error message.  If the parameters file exists, then the code will parse the file, set-up the simulation and then run to completion.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Running python scripts}

\noindent Any python script using the GANDALF python API can be run directly from the command-line as any regular python script.  To run the script PYTHONSCRIPT.py, type \\
\newline
\var{python PYTHONSCRIPT.py} \\
\newline
or, depending on the default version of python on your system (e.g. installing matplotlib with fink on Mac OS X), \\
\newline
\var{python2.7 PYTHONSCRIPT.py} \\
\newline
The code will remain inside the python environment until the script is fully executed.  Python typically checks that the entire file has valid syntax before running each command.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Interactive mode}

\noindent To open the interactive viewer, type : \\
\newline
\var{python analysis/gandalf\_interpreter.py} \\
\newline
or, depending on the default version of python on your system (e.g. installing matplotlib with fink on Mac OS X), \\
\newline
\var{python2.7 analysis/gandalf\_interpreter.py} \\
\newline
\noindent The code should open with a splash screen containing the code title followed by a command prompt of the form '\var{gandalf >'}'.  To run a simulation defined by a parameter file, then type \\
\newline
\var{newsim PARAMSFILE} \\
\newline
\noindent and then \\
\newline
\var{run} \\
\newline
\noindent The current simulation can be plotted at any point by using simple commands such as, for example, \\
\newline
\var{plot x y} \\
\newline
\noindent For a list on available commands, type \var{help} in the command line.  For more detailed information on the functionality of a particular command, type \var{help command}.  For more information on interactive python commands, see Section \ref{S:PYTHONSCRIPT}.

\newpage



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Parameter options} \label{S:PARAMETERS}

\noindent With the exception of a few compiled-time options that are selected in the Makefile, all physics and code feature options are controlled from the parameters file, including the dimensionality, simulation type and SPH options.  We list here all possible parameters, including all possible options for those parameters with a limited choice.  We note that a parameters file need not contain a value for each parameter, in which case the default value is taken.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Main simulation parameters}

\begin{itemize}

\item \var{ndim} : Simulation dimensionality ($1$, $2$ or $3$)

\item \var{sim} : Simulation type \vspace{0.1cm} \\
\begin{tabular}{ll}
sph               & = SPH (+ N-body) algorithm (default : `grad-h' SPH) \\
gradhsph          & = `grad-h' SPH simulation (+ N-body) \\
sm2012sph         & = Saitoh \& Makino (2012) SPH (+ N-body) \\
meshlessfv        & = Meshless Finite-Volume algorithm (default : `mfvmuscl') \\
mfvmuscl          & = Meshless FV MUSCL integration simulation \\
mfvrk             & = Meshless FV Runge-Kutta integration \\
nbody             & = N-body only simulation
\end{tabular}

\item \var{nbody} : Main N-body integration algorithm \\
\begin{tabular}{ll}
lfkdk        & = 2nd-order Leapfrog kick-drift-kick \\
lfdkd        & = 2nd-order Leapfrog drift-kick-drift \\
hermite4     & = 4th-order Hermite scheme \\
hermite4ts   & = Time-symmetric 4th-order Hermite scheme
\end{tabular}

\item \var{ic} : Simulation initial conditions \vspace{0.1cm} \\
\begin{tabular}{ll}
file              & = Load initial conditions from external file \\
bb                & = Boss-Bodenheimer (1979) test \\
binary            & = Simple binary star test \\
binaryacc         & = Binary accretion test \\
blastwave         & = Blastwave test \\
bondi             & = Spherically symmetric Bondi accretion test \\
box               & = Create a uniform box of gas \\
cdiscontinuity    & = Contact discontinuity test \\
disc					& = Accretion disc initial conditions \\
ewaldcylinder     & = Cylinder for 1D Ewald gravity test \\
ewaldsine         & = Sinusoidal density for 3D Ewald gravity test \\
ewaldslab         & = Slab for 2D Ewald gravity test \\
gresho            & = Gresho-Chan vortex test \\
khi               & = Kelvin-Helmholtz instability test \\
noh               & = Noh problem initial conditions \\
plummer           & = Plummer sphere test \\
quadruple         & = Simple hierarchical quadruple star test \\
sedov             & = Sedov blast-wave test \\
shearflow         & = Shear flow test \\
shocktube         & = Shocktube test \\
soundwave         & = 1D soundwave perturbation \\
sphere            & = Uniform density sphere \\
spitzer           & = Ionised bubble Spitzer expansion test \\
triple            & = Simple hierachical triple star test \\
turbcore          & = Turbulent spherical, self-gravitating core \\
python            & = Generate initial conditions from python
\end{tabular}

\item \var{run\_id}  : Simulation run id string

\item \var{in\_file} : Input filename (when \var{ic = file})

\item \var{in\_file\_form} : Format of initial conditions file \\
\begin{tabular}{ll}
column           & = Simple column data format \\
sf/seren\_form   & = SEREN ASCII format \\
su/seren\_unform & = SEREN binary format
\end{tabular}

\item \var{out\_file\_form} : Format of outputted snapshot files \\
\begin{tabular}{ll}
column           & = Simple column data format \\
sf/seren\_form   & = SEREN ASCII format \\
su/seren\_unform & = SEREN binary format
\end{tabular}

\item \var{tend} : Termination time of the simulation (given in {\var tunit}s)

\item \var{tmax\_wallclock} : Maximum allowed wallclock time for simulation before being terminated

\item \var{dt\_snap} : Snapshot time interval (given in {\var tunit}s)

\item \var{tsnapfirst} : Time of first snapshot (given in {\var tunit}s)

\item \var{Nstepsmax} : Maximum no. of steps in simulation before termination

\item \var{noutputstep} : Frequency of screen output (in units of integer steps)

\item \var{ndiagstep} : No. of complete block steps between diagnostic output

\item \var{nrestartstep} : No. of full block steps before producing restart dump

\item \var{litesnap} : Output `lite' snapshots (for generating movies)? ($0$ or $1$)

\item \var{dt\_litesnap} : Lite snapshot time interval (given in {\var tunit}s)

\item \var{tlitesnapfirst} : Time of first lite snapshot (given in {\var tunit}s)

\end{itemize}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Unit parameters} \label{SS:UNITPARAMS}

\begin{itemize}
\item \var{dimensionless} : Are all quantities dimensionless? ($0$ or $1$)

\item \var{routunit} : Position unit \vspace{0.1cm} \\
\begin{tabular}{ll}
pc/kpc/mpc & = parsec/kiloparsec/megaparsec \\
au         & = astronomical unit \\
r\_sun     & = Solar radius \\
r\_earth   & = Earth radius \\
cm/m/km    & = centimetre/metre/kilometre
\end{tabular}

\item \var{moutunit} : Mass unit \vspace{0.1cm} \\
\begin{tabular}{ll}
m\_sun          & = Solar mass \\
m\_jup/m\_earth & = Jupiter mass/Earth mass \\
g/kg            & = gram/kilogram
\end{tabular}

\item \var{toutunit} : Time unit \\
\begin{tabular}{ll}
yr/myr/gyr & = year/megayear/gigayear \\
day        & = day \\
sec        & = second
\end{tabular}

\item \var{voutunit} : Velocity unit \\
\begin{tabular}{ll}
cm\_s/m\_s/km\_s & = centimetres/metres/kilometres per second \\
au\_yr           & = astronomical units per year
\end{tabular}

\item \var{aoutunit} : Acceleration unit \\
\begin{tabular}{ll}
cm\_s2/m\_s2/km\_s2 & = cm/m/km per second squared \\
au\_yr2             & = astronomical units per year squared
\end{tabular}

\item \var{rhooutunit} : Density unit \\
\begin{tabular}{ll}
m\_sun\_pc3 & = Solar masses per parsec cubed \\
kg\_m3      & = kilogrammes per metre cubed \\
g\_cm3      & = grammes per centimetre cubed
\end{tabular}

\item \var{sigmaoutunit} : Column/surface density unit \\
\begin{tabular}{ll}
kg\_m2      & = kilogrammes per meter squared
\end{tabular}

\item \var{pressoutunit} : Pressure unit \\
\begin{tabular}{ll}
Pa    & = pascals/newtons per square metre \\
bar   & = bars
\end{tabular}

\item \var{foutunit} : Force unit \\
\begin{tabular}{ll}
N    & = newtons \\
dyn  & = dynes
\end{tabular}

\item \var{Eoutunit} : Energy unit \\
\begin{tabular}{ll}
J/GJ      & = joules/gigajoules \\
erg       & = ergs \\
eV        & = electron volts \\
10\^40erg & = $10^{40}$ ergs
\end{tabular}

\item \var{momoutunit} : Momentum unit \\
\begin{tabular}{ll}
m\_sunkm\_s  & = Solar masses kilometres per second \\
m\_sunau\_yr & = Solar masses A.U. per year \\
kgm\_s       & = Kilogram metres per second \\
gcm\_s       & = Gram centimetres per second
\end{tabular}

\item \var{angmomoutunit} : Angular momentum unit \\
\begin{tabular}{ll}
m\_sunkm2\_s  & = Solar masses kilometres squared per second \\
m\_sunau2\_yr & = Solar masses A.U. squared per year \\
kgm2\_s       & = Kilogram metres squared per second \\
gcm2\_s       & = Gram centimetres squared per second
\end{tabular}

\item \var{angveloutunit} : Angular velocity unit \\
\begin{tabular}{ll}
rad\_s & = Radians per second
\end{tabular}

\item \var{dmdtoutunit} : Mass (accretion) rate unit \\
\begin{tabular}{ll}
m\_sun\_myr  & = Solar masses per megayear \\
m\_sun\_yr   & = Solar masses per year \\
kg\_s        & = kilogrammes per second \\
g\_s         & = grammes per second
\end{tabular}

\item \var{Loutunit} : Luminosity unit \\
\begin{tabular}{ll}
L\_sun       & = Solar luminosity \\
W            & = watts \\
erg\_s       & = ergs per second
\end{tabular}

\item \var{kappaoutunit} : Mass opacity unit \\
\begin{tabular}{ll}
m2\_kg       & = metre squared per kilogram \\
cm2\_g       & = centimetre squared per gram
\end{tabular}

\item \var{Boutunit} : Magnetic field unit \\
\begin{tabular}{ll}
tesla     & = tesla \\
gauss     & = gauss
\end{tabular}

\item \var{Qoutunit} : Charge unit \\
\begin{tabular}{ll}
C    & = coulomb \\
e    & = electron charge
\end{tabular}

\item \var{Jcuroutunit} : Current density unit \\
\begin{tabular}{ll}
C\_s\_m2    & = coloumb per second per metre squared
\end{tabular}

\item \var{uoutunit} : Specific energy unit \\
\begin{tabular}{ll}
J\_kg  & = Joules per kilogram \\
erg\_g & = ergs per gram
\end{tabular}

\item \var{dudtoutunit} : Heating rate unit \\
\begin{tabular}{ll}
J\_kg\_s  & = Joules per kilogram per second \\
erg\_g\_s & = ergs per gram per second
\end{tabular}

\item \var{tempoutunit} : Temperature unit \\
\begin{tabular}{ll}
K & = Kelvin
\end{tabular}

\end{itemize}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Integration and timestep parameters}

\begin{itemize}

\item \var{accel\_mult} : Acceleration timestep multiplier

\item \var{courant\_mult} : Courant timestep multiplier

\item \var{visc\_mult} : Viscosity timestep multiplier

\item \var{nbody\_mult} : N-body timestep multiplier

\item \var{subsys\_mult} : Sub-system N-body timestep multiplier

\item \var{Nlevels} : No. of initial timestep levels

\item \var{level\_diff\_max} : Maximum allowed SPH neighbour timestep difference

\item \var{sph\_single\_timestep} : Constrain all SPH particles to a single timestep level

\item \var{nbody\_single\_timestep} : Constrain all N-body particles to a single timestep level


\end{itemize}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Hydrodynamical parameters}

\begin{itemize}
\item \var{hydro\_forces} : Compute hydro forces?  ($1$ or $0$)

\item \var{gas\_eos} : Gas particles equation-of-state \\
\begin{tabular}{ll}
energy\_eqn & = Solve energy equation \\
isothermal  & = Isothermal EOS \\
barotropic  & = Barotropic EOS (i.e. for mimicing isothermal + adiabatic phase during protostellar collapse) \\
barotropic2  & = Similar to barotropic, but using discrete power laws rather than smooth change \\
rad\_ws     & = EOS relating to Stamatellos et al. (2007) cooling method \\
disc\_locally\_isothermal & = Locally isothermal equation of state, to be used with the disc setup.
\end{tabular}


\item \var{energy\_integration} : Energy integration scheme (only applicable if solving the energy equation) \\
\begin{tabular}{ll}
null        & = Energy equation not integrated separately \\
rad\_ws     & = Integrate energy terms using Stamatellos et al. (2007) method
\end{tabular}

\item \var{energy\_mult} : Explicit energy integration timestep multiplier

\item \var{gamma\_eos} : Ratio of specific heats for gas

\item \var{temp0} : (Isothermal) temperature (isothermal or barotropic EOS)

\item \var{mu\_bar}    : Mean gas particle mass (in units of hydrogen mass)

\item \var{rho\_bary} : Adiabatic density turnover in barotropic EOS (in g/cm\^3)

\item \var{eta\_eos}   : Polytropic exponent (for barotropic EOS)

\item \var{radws\_table} : Name of EOS file for Stamatellos et al. (2007) cooling method

\item \var{temp\_ambient} : Ambient temperature (for rad\_ws method)

\item \var{lombardi\_method} : Use the Lombardi et al. (2015) metric for rad\_ws method


%\item \var{Acool} : Simple cooling rate factor

%\item \var{u\_eq} : Equilibrium internal energy (for simple cooling law)

\end{itemize}



\subsection{SPH parameters}

\begin{itemize}

\item \var{sph\_integration} : SPH particle integration scheme \\
\begin{tabular}{ll}
lfkdk   & = 2nd-order Leapfrog kick-drift-kick \\
lfdkd   & = 2nd-order Leapfrog drift-kick-drift
\end{tabular}

\item \var{kernel} : SPH kernel function \\
\begin{tabular}{ll}
m4        & = M4 Cubic spline kernel \\
quintic   & = Quintic spline kernel \\
gaussian  & = Gaussian kernel (truncated at 3h)
\end{tabular}

\item \var{tabulated\_kernel} : Tabulate kernel function  ($1$ or $0$)

\item \var{h\_fac} : Particles-per-smoothing length factor (eta in papers)

\item \var{h\_converge} : Smoothing length iteration convergence tolerance

\end{itemize}





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Artificial viscosity parameters}

\begin{itemize}

\item \var{avisc} : Artificial viscosity options \\
\begin{tabular}{ll}
none  & = No artificial viscosity \\
mon97 & = Monaghan (1997) viscosity
\end{tabular}

\item \var{acond} : Artificial conductivity options \\
\begin{tabular}{ll}
none        & = No artificial conductivity \\
price2008   & = Price (2008) conductivity \\
wadsley2008 & = Wadsley et al. (2008) conductivity
\end{tabular}

\item \var{time\_dependent\_avisc} : Aritificial viscosity switch \\
\begin{tabular}{ll}
none     &= No switch, artificial viscoisty always set to \var{alpha\_visc}.\\
mm97     &= Morris \& Monaghan (1997) switch \\
cd2010   &= Advanced switch by Cullen \& Dehnen (2010), using the balsar swith
\end{tabular}

\item \var{alpha\_visc} : (Maximum) value of alpha viscosity parameter

\item \var{alpha\_visc\_min} : Minimum value of alpha for time-dependent viscosity

\item \var{beta\_visc} : Value of beta viscosity as a multiple of alpha


\end{itemize}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Meshless finite-volume parameters}

\begin{itemize}

\item \var{riemann\_solver} : Riemann solver in FV scheme \\
\begin{tabular}{ll}
exact  & = Exact Riemann solver (e.g. Toro 1999) \\
hllc & = HLLC approximate Riemann solver
\end{tabular}

\item \var{slope\_limiter} : Slope limiter for TVD condition \\
\begin{tabular}{ll}
null          & = No limiting \\
zeroslope     & = Set all slopes to zero (effectively 1st order Godunov) \\
balsara2004   & = Balsara (2004) slope-limiter \\
springel2009  & = Original AREPO (Springel 2009) slope limiter \\
tess2011      & = TESS slope limtier \\
gizmo         & = Original GIZMO paper (Hopkins 2015) slope limiter \\
minmod        & = simplified implementation of minmod slope limiter
\end{tabular}

\item \var{zero\_mass\_flux} : Use Meshless-Finite Mass (MFM) scheme to prevent mass-flux between particles?  Should be activated if self-gravity or star/sink particles are used in conjunction with any MFV scheme.  ($1$ or $0$)
\item \var{static\_particles} : Use static particles (Eulerian approach)?  ($1$ or $0$)

\item \var{time\_step\_limiter} : Use an additional limiter on the time-steps \\
\begin{tabular}{ll}
none		 & = No additional limiter \\
conservative & = Predictive global timestep limiter of Springel (2009) \\
simple		 & = Saitoh \& Makino type limiter. Reduces the time step if it is detected to be too large.
\end{tabular}

\item \var{shear\_visc} Enable viscosity with constant kinematic shear viscoity
\item \var{bulk\_visc}  Enable viscosity with constant kinematic bulk  viscoity

\end{itemize}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Gravitational parameters}

\begin{itemize}

\item \var{self\_gravity} : Compute gravitational forces?   ($1$ or $0$)

\item \var{kgrav} : Direction of (external) gravitational acceleration ($0$, $1$ or $2$)

\item \var{grav\_kernel} : Form of gravitational softening \\
\begin{tabular}{ll}
mean\_h   & = Mean smoothing length softening
\end{tabular}

\item \var{external\_potential} : External gravitational potential \\
\begin{tabular}{ll}
none     & = No external potential \\
vertical & = Constant gravitational field \\
plummer  & = Plummer background potential
\end{tabular}

\item \var{avert} : Veritcal (constant) gravitational acceleration

\item \var{rplummer\_extpot} : Background Plummer potential radius

\item \var{mplummer\_extpot} : Background Plummer potential mass

\end{itemize}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Neighbour searching and tree parameters}

\begin{itemize}

\item \var{neib\_search} : Neighbour searching algorithm \vspace{0.1cm} \\
\begin{tabular}{ll}
bruteforce & = Brute-force (i.e. summation over all particles) \\
kdtree     & = Balanced kd-binary tree \\
octtree    & = Barnes-Hut octal tree
\end{tabular}

\item \var{gravity\_mac} : Gravity-tree cell-opening criteria (N.B. always defualts to geometric for now). The eigenmac and gadget2 criteria are more accurate for a given computational cost. \vspace{0.1cm} \\
\begin{tabular} {ll}
geometric & = Standard Barnes-Hut geomtric opening angle criterion. This can fail in pathalogical cases,\\
	      & \quad such as when the forces nearly balance. \\
eigenmac  & = Compute eigenvalues of quadrupole moment tensor for MAC (Hubber et al. 2011) \\
gadget2   & = Relative opening criterion based on comparing estimate of force error to the previous \\
		  & \quad accleleration (Springel 2005).
\end{tabular}

\item \var{multipole} : Multipole expansion for tree-gravity \vspace{0.1cm} \\
\begin{tabular} {ll}
monopole       & = Monopole-only terms for cell gravity \\
quadrupole     & = Include quadrupole moment terms for cell gravity \\
fast\_monopole & = Compute monpoles more efficiently using Taylor expansion about cell COM
\end{tabular}

\item \var{Nleafmax} : Maximum no. of particles allowed in tree leaf cell

\item \var{ntreebuildstep} : Integer steps inbetween tree re-builds

\item \var{ntreestock} : Integer steps inbetween tree re-stocks

\item \var{thetamaxsqd} : Maximum tree gravitational walk opening angle (squared)

\item \var{macerror} : MAC error tolerance for individual cells

\end{itemize}





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{N-body parameters}

\begin{itemize}

\item \var{sub\_systems} : Identify and integrate sub-systems separately? ($0$ or $1$)

\item \var{sub\_system\_integration} : Main N-body integration algorithm \\
\begin{tabular}{ll}
lfkdk        & = 2nd-order Leapfrog kick-drift-kick \\
hermite4     & = 4th-order Hermite scheme \\
hermite4ts   & = Time-symmetric 4th-order Hermite scheme
\end{tabular}

\item \var{Npec} : No. of ${\rm P(EC)^n}$ iterations in time-symmetric scheme (if non time-symmetric scheme is used, automatically sets to $1$)

\item \var{nbody\_softening} : Use SPH kernel-softening between star particles? ($0$ or $1$)

%\item \var{perturbers} : Use perturbers when calculating sub-system quantities ($1$ or $0$)

\item \var{binary\_stats} : Output binary statistics? ($1$ or $0$)

\item \var{nsystembuildstep} : Integer steps inbetween re-building the sub-system tree.

\item \var{gpefrac} : Maximum fraction of total gravitational potential energy from external sources to allow sub-system.

%\item \var{gpesoft} : Maximum fraction of total gravitational potential energy from ...

%\item \var{gpehard} : Maximum fraction of total gravitational potential energy from ...

\end{itemize}





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Sink particle parameters}

\begin{itemize}

\item \var{sink\_particles}   : Do stars/sinks accrete?  ($0$ or $1$)

\item \var{create\_sinks}     : Create new sink particles?  ($0$ or $1$)

\item \var{smooth\_accretion} : Use smooth accretion?  ($0$ or $1$)

\item \var{fixed\_sink\_mass} : Fixed sink mass, even when accreting? ($0$ or $1$)

\item \var{extra\_sink\_output} : Extra output of sink particles? ($0$ or $1$)

\item \var{rho\_sink}         : Sink particle creation density (in cgs units)

\item \var{alpha\_ss}         : Sunyaev-Shakura alpha for smooth disc accretion

\item \var{sink\_radius}      : Sink particle radius (in units of smoothing length)

\item \var{smooth\_accrete\_frac} : Smooth accretion instantaneous accretion mass frac.

\item \var{smooth\_accrete\_dt} : Smooth accretion instantaneous accretion timestep frac.

\item \var{sink\_radius\_mode} : How to calculate new sink radius \\
\begin{tabular}{ll}
hmult   & = sink radius a multiple of SPH particle smoothing length \\
fixed   & = sink radius is fixed for all new sinks \\
\end{tabular}


\end{itemize}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Radiation parameters}

\begin{itemize}

\item \var{radiation} : Main radiation algorithm used \\
\begin{tabular}{ll}
none            &= No radiation field \\
ionisation      &= Multiple source ionising radiation
\end{tabular}

\item \var{Nphoton}   : No. of photon packets (for Monte-Carlo radiation transport)
\item \var{mu\_ion}   : Mean-gas particle mass for ionised gas
\item \var{temp\_ion} : Temperature of ionised gas
\item \var{arecomb}   : Recombination coefficient (in cgs units)
\item \var{Ndotmin}   : No. of ionising photons per second
\item \var{NLyC}      : No. of ionising photons per second

\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Radiative feedback parameters}

\begin{itemize}

\item \var{rad\_fb} : Turn on radiative feedback for supported specific internal energy integration methods.
\item \var{ambient\_heating}   : Add heating from an ambient background defined by \var{temp\_ambient}
\item \var{disc\_heating} : Add heating from a temperature profile with a central star (1) or binary (2)
\item \var{sink\_heating} : Add accretion heating from sinks (except central objects from \var{disc\_heating})
\item \var{sink\_fb} : Type of heating from sinks (except central objects) \\
\begin{tabular}{ll}
continuous    &= Heating from accretion as mass is accreted \\
episodic      &= NOT IMPLEMENTED
\end{tabular}
\item \var{r\_smooth} : Smoothing radius for \var{disc\_heating} temperature profile
\item \var{temp\_q} : Power index for primary in \var{disc\_heating} temperature profile
\item \var{temp\_q\_secondary} : Power index for secondary in \var{disc\_heating} temperature profile
\item \var{temp\_au} : Temperature at 1 AU for primary in \var{disc\_heating} temperature profile
\item \var{temp\_au\_secondary} : Temperature at 1 AU for secondary in \var{disc\_heating} temperature profile
\item \var{f\_acc} : Fraction of energy converted to heating from accretion (use with \var{sink\_heating})
\item \var{r\_star} : Radius of star ($M_{sink} > 80~M_{jup}$) in solar radii (use with \var{sink\_heating})
\item \var{r\_bdwarf} : Radius of brown dwarf ($13 > M_{sink} > 80~M_{jup}$) in solar radii (use with \var{sink\_heating})
\item \var{r\_planet} : Radius of planet ($M_{sink} < 13~M_{jup}$) in solar radii (use with \var{sink\_heating})

\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Boundary parameters}

\begin{itemize}

\item \var{boundary\_lhs[0]} : Boundary conditions for LHS of x-dimension
\item \var{boundary\_rhs[0]} : Boundary conditions for RHS of x-dimension
\item \var{boundary\_lhs[1]} : Boundary conditions for LHS of y-dimension
\item \var{boundary\_rhs[1]} : Boundary conditions for RHS of y-dimension
\item \var{boundary\_lhs[2]} : Boundary conditions for LHS of z-dimension
\item \var{boundary\_rhs[2]} : Boundary conditions for RHS of z-dimension
For all boundaries: \\
\begin{tabular}{ll}
open     & = open boundaries (i.e. extends to infinity) \\
periodic & = periodic wrapping between LHS \& RHS boundary \\
wall     & = wall at boundary (i.e. reflection of particles)
%mirror   & = mirror boundary (i.e. 'ghost' reflections)
\end{tabular}

\item \var{boxmin[0]} : Location of LHS x-boundary
\item \var{boxmax[0]} : Location of RHS x-boundary
\item \var{boxmin[1]} : Location of LHS y-boundary
\item \var{boxmax[1]} : Location of RHS y-boundary
\item \var{boxmin[2]} : Location of LHS z-boundary
\item \var{boxmax[2]} : Location of RHS z-boundary


\end{itemize}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Initial conditions parameters}

\begin{itemize}

\item \var{particle\_distribution} : Particle configuration when generating uniform density fluids on the fly \\
\begin{tabular}{ll}
random         &= Particle positions generated with random number generator \\
cubic\_lattice &= Particle positions generated on a uniform cubic lattice \\
hexagonal\_lattice &= Particle positions generated on a hexagonal closed-packed array
\end{tabular}

\item \var{smooth\_ic} : Smooth any particle quantities around discontinuities
\item \var{com\_frame} : Translate ICs to COM frame before starting simulation
%\item \var{field\_type} : ..
%\item \var{gridsize}    : ..
\item \var{Nhydro}       : No. of hydrodynamical particles
\item \var{Nhydromax}    : Maximum no. of hysrodynamical particles
\item \var{Nstar}        : No. of star particles
\item \var{Nstarmax}     : Maximum no. of star particles

\item \var{Nlattice1[0]} : No. of ptcls on lattice 1 in x-dimension
\item \var{Nlattice1[1]} : No. of ptcls on lattice 1 in y-dimension
\item \var{Nlattice1[2]} : No. of ptcls on lattice 1 in z-dimension
\item \var{Nlattice2[0]} : No. of ptcls on lattice 2 in x-dimension
\item \var{Nlattice2[1]} : No. of ptcls on lattice 2 in y-dimension
\item \var{Nlattice2[2]} : No. of ptcls on lattice 2 in z-dimension

\item \var{vfluid1[0]}   : x-velocity of fluid 1
\item \var{vfluid1[1]}   : y-velocity of fluid 1
\item \var{vfluid1[2]}   : z-velocity of fluid 1
\item \var{vfluid2[0]}   : x-velocity of fluid 2
\item \var{vfluid2[1]}   : y-velocity of fluid 2
\item \var{vfluid2[2]}   : z-velocity of fluid 2

\item \var{rhofluid1}    : Density of fluid 1
\item \var{rhofluid2}    : Density of fluid 2

\item \var{press1}       : Pressure of fluid 1
\item \var{press2}       : Pressure of fluid 2

\item \var{amp}          : Amplitude of applied perturbation
\item \var{lambda}       : Wavelength of applied perturbation

\item \var{kefrac}       : Fraction of energy that is kinetic (Sedov test)
\item \var{radius}       : Radius of cloud
\item \var{angvel}       : Angular velocity of cloud (in radians per second)
\item \var{mcloud}       : Mass of cloud

\item \var{rplummer}     : Plummer radius
\item \var{mplummer}     : Total mass of plummer sphere

\item \var{rstar}        : (Softening) radius of star particles

\item \var{cdmfrac}      : Fraction of mass in cdm particles
\item \var{gasfrac}      : Fraction of mass in gas particles
\item \var{starfrac}     : Fraction of mass in star particles

\item \var{m1}           : Mass of star 1
\item \var{m2}           : Mass of star 2
\item \var{m3}           : Mass of star 3
\item \var{m4}           : Mass of star 4

\item \var{abin}         : Semi-major axis of binary orbit 1
\item \var{abin2}        : Semi-major axis of binary orbit 2
\item \var{ebin}         : Orbital eccentricity of binary orbit 1
\item \var{ebin2}        : Orbital eccentricity of binary orbit 2

\item \var{phirot}       : Phi Euler rotation angle
\item \var{thetarot}     : Theta Euler rotation angle
\item \var{psirot}       : Psi Euler rotation angle
\item \var{vmachbin}     : Speed of binary COM through ambient gas

\item \var{alpha\_turb}  : Turbulent energy (as multiple of gravitational energy)
\item \var{power\_turb}  : Power spectrum slope of initial turbulent velocity field

\item \var{asound}       : Sound speed
\item \var{zmax}         : ??

\item \var{DiscIc}: There are multiple parameters prefixed with this string; they contol the accretion disc setup. Please refer to section \ref{sec:disc} for an extensive explanation.

\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Regularised initial conditions parameters}

\begin{itemize}

\item \var{regularise\_particle\_ics} : Regularise particle initial conditions before main simulation
\item \var{regularise\_smooth\_density} : Regularisation uses the smoothed density field
\item \var{Nreg}         : No. of regularisation iterations
\item \var{alpha\_reg}   : Glass-creation regularisation factor (0 = no glass; $> 0$ = glass-like)
\item \var{rho\_reg}     : Density-field regularisation factor ($> 0$ = iterates to required density field)
\item \var{Nreg}         : No. of regularisation iterations

\end{itemize}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Random number generator parameters}

\begin{itemize}

\item \var{rand\_algorithm} : Random number generator algorithm \\
\begin{tabular}{ll}
none     & = No algorithm selected; use intrinsic generator \\
         & \,\,\,(not recommended since this is system dependent and is not reproducable on different machines) \\
xorshift & = Xorshift generator (see Numerical recipes, Ed 3, Chapter 7 for details)
\end{tabular}

\item \var{randseed}        : Random number seed

\end{itemize}





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{MPI parameters}

\begin{itemize}

\item \var{mpi\_decomposition}  : Mode of MPI decomposition \\
\begin{tabular}{ll}
kdtree & = Use simple KD-tree decomposition
\end{tabular}
\item \var{pruning\_level\_min} : Minimum level to prune exported trees
\item \var{pruning\_level\_max} : Maximum level to prune exported trees

\end{itemize}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Python viewer parameters}

\begin{itemize}

\item \var{dt\_python} : Time interval (in seconds) between view window updates

\end{itemize}


\newpage



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Input and output file formats} \label{S:IO}

GANDALF supports several simple file formats for reading in initial conditions or outputting snapshots for visualisation and analysis.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Column format}
The column format is the simplest snapshot file format in GANDALF and consists of a simple header plus a fixed column-data format with particle data.  It is designed for ease-of-use in quickly generating initial conditions, or perhaps porting initial conditions from other codes/generators.  It is certainly not designed as a long-term format for use in large simulations, partly due to being an ASCII file which can quickly get unfeasibly large for big simulations.

The header comprises of four lines, each with one variable in the following order :
\begin{itemize}
\item \var{Nhydro} : No. of hydro particles
\item \var{Nstar}  : No. of star/sink particles
\item \var{ndim}   : Dimensionality of snapshot data
\item \var{t}      : Time of snapshot
\end{itemize}
Each of the 4 header variables is preceded by a hash (`\#') in order that simple plotting programs (e.g. gnuplot) will ignore the headers allowing the particle data to be simply plotted.  Older versions of GANDALF may not have this hash but nevertheless should still easily be read-in by the code.

The particle data is then split up into \var{Nhydro} + \var{Nstar} rows, hydro particles first and then star particles.  The data columns are ordered as :
\begin{itemize}
\item \var{x [,y, z]} : Particle position vectors (depending on value of \var{ndim})
\item \var{vx [,vy, vz]} : Particle velocity vectors (depending on value of \var{ndim})
\item \var{m} : Particle mass
\item \var{h} : Particle smoothing length (for stars also; radius = \var{2/3\,h} depending on kernel)
\item \var{rho} : Particle density (set to 0 for stars)
\item \var{u} : Particle specific internal energy (set to 0 for stars)
\end{itemize}
Since the column format was only designed for simple IO and plotting, it is a rather minimalistic format that does not necessarily contain all the required data or options for full simulations.  For example, all stars are automatically set to being sinks.  Therefore it is not recommended for long-term use other than simple plotting (e.g. for debugging or developing/transfering ICs).



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{SEREN format}
\label{Sec:SerenForm}
The SEREN format is a legacy format from the original SEREN code which is currently the main format for use in GANDALF.  There are two different version of the SEREN format;
\begin{itemize}
\item \var{seren\_form} or \var{sf} : `Formatted' (i.e. ASCII)
\item \var{seren\_unform} or \var{su} : `Unformatted' (i.e. binary)
\end{itemize}
The format consists of a large multi-part header with different data types, including information about the units used as well as the particle data and star/sink data arrays.  The full header array sizes are given below.  In the formatted/ASCII case, each array element is printed on a separate line but is compressed as a single stream.  For any fields not given below, assume the quantity to be unused and therefore defaulted to zero in the header.  The header reads as follows:
\begin{itemize}
\item \var{format\_id} : A string identifying the exact format type (for verification)
\item \var{pr} : Floating point precision (4 = single precision, 8 = double precision)
\item \var{ndim} : Spatial dimensionality
\item \var{vdim} : Velocity dimensionality
\item \var{bdim} : Magnetic field dimensionality
\item \var{idata[50]} : Integer variables
\begin{itemize}
\item idata[0] = \var{Nhydro} : No. of hydro particles (SEREN = \var{pgas})
\item idata[1] = \var{Nstar}  : No. of star/sink particles (SEREN = \var{stot})
\item idata[2] = N/A in GANDALF (SEREN = \var{pboundary} : No. of static boundary particles)
\item idata[3] = N/A in GANDALF (SEREN = \var{picm} : No. of `inter-cloud medium' particles)
\item idata[4] = \var{Ngas}  : No. of self-gravitating gas particles (SEREN = \var{pgas})
\item idata[5] = \var{Ncdm}  : No. of self-gravitating cdm particles (SEREN = \var{pcdm})
\item idata[6] = \var{Ndust} : No. of dust particles (SEREN = \var{pdust})
\item idata[7] = N/A in GANDALF (SEREN = \var{pion} : No. of  particles)
\item idata[19] = \var{nunit} : No. of unit variables in header
\item idata[20] = \var{ndata} : No. of data arrays in snapshot file
\item idata[29] = N/A in GANDALF (SEREN = \var{dmdt\_range} : No. of accretion rate variables in sink array)
\item idata[30] = N/A in GANDALF (SEREN = \var{pgas\_orig} : Original no. of gas particles)
\item idata[31] = N/A in GANDALF (SEREN = \var{pp\_gather} : Average/exact no. of neighbours)
\item idata[39] = N/A in GANDALF (SEREN = \var{rank} : MPI process that created the file)
\item idata[40] = N/A in GANDALF (SEREN = \var{Nmpi} : Total no. of MPI processes)
\end{itemize}

\item \var{ilpdata[50]} : Long integer variables
\begin{itemize}
\item ilpdata[0] = \var{Noutsnap} : No. of snapshot files created (SEREN = \var{snapshot})
\item ilpdata[1] = \var{Nsteps} : No. of complete integration steps (SEREN = \var{nsteps})
\item ilpdata[2] = N/A in GANDALF (SEREN = \var{ntempnext} : Integer time for next temporary snapshot)
\item ilpdata[3] = N/A in GANDALF (SEREN = \var{ndiagnext} : Integer time for next diagnostic output)
\item ilpdata[4] = N/A in GANDALF (SEREN = \var{nsnapnext} : Integer time for next integer snapshot)
\item ilpdata[5] = N/A in GANDALF (SEREN = \var{nsinknext} : Integer time for next sink output)
\item ilpdata[10] = \var{Noutlitesnap} : No. of lite snapshots (N/A in SEREN)
\end{itemize}
\item \var{rdata[50]} : Standard precision floating-point variables
\begin{itemize}
\item rdata[0] = \var{h\_fac} : `No. of particles per smoothing length' factor in h-rho iteration (SEREN = \var{h\_fac})
\item rdata[1] = N/A in GANDALF (SEREN = \var{gamma} : Ratio of specific heats for gas)
\item rdata[2] = N/A in GANDALF (SEREN = \var{mu\_bar} : Mean gas particle mass (in units of $m_h$))
\item rdata[3] = N/A in GANDALF (SEREN = \var{hmin} : Minimum smoothing length)
\end{itemize}
\item \var{ddata[50]} : Double precision floating-point variables
\begin{itemize}
\item ddata[0] = \var{t} : Simulation time when snapshot was created (SEREN = \var{time})
\item ddata[1] = \var{tsnaplast} : Time that previous snapshot was created (SEREN = \var{lastsnap})
\item ddata[2] = \var{mmean} : Average mass of gas particles (SEREN = \var{mgas\_orig})
\item ddata[10] = \var{tlitesnaplast} : Time when previous lite snapshot was created (N/A in SEREN)
\end{itemize}
\item \var{unit\_data[nunit]} : Strings identifying the chosen input unit for each physical quantity.  Note that this does not necessarily have to be the same as the chosen output units (See Section \ref{S:UNITS}).  In unformatted/binary files, the strings are set to a default maximum length of 20 characters (with white space padding at the end).
\item \var{data\_id[ndata]} : Strings identifying the physical quantities contained in the particle arrays that follow the header.  Note that vector quantities are written grouping together the different components for each component, rathern than different ndim scalar values (e.g., in 2d x and y of a given particle are written one after the other, and after come the x and y of another particle). As with the \var{unit\_data} array using unformatted/binary file, the strings are set to a default maximum length of 20 characters (with white space padding at the end).
\begin{itemize}
\item \var{porig} : Original (unique) particle ids
\item \var{r} : Hydro particle positions
\item \var{m} : Hydro particle masses
\item \var{h} : Hydro particle smoothing lengths
\item \var{v} : Hydro particle velocities
\item \var{rho} : Hydro particle densities
\item \var{u} : Hydro particle specific internal energies
\item \var{sink\_v1} : Star/sink particle data (N.B. a more complicated data structure due to the various star/sink properties in SEREN; recommended to look at the source code for more info)
\end{itemize}
\item \var{typedata} : Information on each particle array, such as the variable type, unit type, etc..  For the data held in array \var{idata}, the data information in each field is :
\begin{itemize}
\item typedata[idata][0] : Dimensionality of array (i.e. 1 to ndim)
\item typedata[idata][1] : i.d. of first particle in array in Fortran (i.e. normally 1 as opposed to 0 in C/C++)
\item typedata[idata][2] : i.d. of last particle in array in Fortran (i.e. normally Nhydro as opposed to Nhydro - 1 as in C/C++)
\item typedata[idata][3] : Variable type of array (1 : bool; 2 : integer; 3 : long integer; 4 : float;  5 : double)
\item typedata[idata][4] : Physical unit of array (same i.d. as \var{unit\_data} array index plus one due to Fortran array convention)
\end{itemize}
\end{itemize}
The headers, in particular \var{data\_id} and \var{typedata}, contain all the information about what type of data is contained in the particle arrays. The particle arrays come immediately after the end of the header.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{SEREN `lite' format}
The SEREN `lite' format is a minimised version of the SEREN format designed for creating many snapshots for movies.  It has various restrictions which do NOT make it useful for other forms of analysis or restarting simulation, such as :
\begin{itemize}
\item No velocity information (only contains position, mass, smoothing length, density and specific internal energy)
\item Hard-wired to single precision to obtain the smallest file-size in binary format
\end{itemize}
It can used simulataneously with the other more complete formats, which can instead be used for the analysis or for restarting.  There are three parameters which control the usage of the lite formats, \var{litesnap}, \var{dt\_litesnap} and \var{tlitesnapfirst} (See parameter tables)


\newpage



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Generating initial conditions}

There are three main ways of generating initial conditions in GANDALF.

%One important difference between SEREN and GANDALF is that initial conditions can be generated on-the-fly in GANDALF, i.e. while running the code, unlike SEREN where initial conditions always had to be generated before the simulation.  While they are still situations where it is appropriate to prepare initial conditions prior to running the simulations, it is more convenient to generate the initial conditions at the same time as running the simulation. It's also possible to generate the initial conditions in python (see next section).


\subsection{`On-the-fly' initial conditions}

It is possible to generate initial conditions `on-the-fly' using internal subroutines in GANDALF.  At present, the following initial conditions are included in the code : \newline

\noindent Hydrodynamical simulations : \\
\newline
\begin{tabular}{ll}
- bb             &: Boss-Bodenheimer test \\
- binaryacc      &: Binary accretion simulation \\
- blastwave      &: 1D blastwave test \\
- bondi          &: Spherically symmetric Bondi accretion test \\
- box            &: Create uniform box \\
- cdiscontinuity &: Contact-discontinuity test \\
- disc &: Accretion disc. See section \ref{sec:disc}. \\
- ewald          &: Several tests of Ewald gravity \\
- gresho         &: Gresho vortex test \\
- khi            &: Kelvin-helmholtz instability \\
- noh            &: Noh shock test \\
- plummer        &: Plummer sphere (stars + gas, or just gas) \\
%- rti            &: Rayleigh-Taylor instability \\
- sedov          &: Sedov blastwave test \\
- shearflow      &: Shearing flow test \\
- shocktube      &: Simple two-fluid shocktube test \\
- soundwave      &: Simple 1D sound-wave perturbation test \\
- sphere         &: Create uniform density sphere \\
- spitzer        &: Spitzer expansion of HII region test \\
- turbcore       &: Create uniform density core with turbulent velocity field
\end{tabular}
\newline
\newline

\noindent N-body simulations : \\
\newline
\begin{tabular}{ll}
- binary    &: Simple circular binary system test \\
- burrau    &: Burrau Pythagorean test \\
- figure8   &: Simple 3-body figure-8 test \\
- plummer   &: Plummer sphere (stars + gas, or just stars) \\
- quadruple &: Simple hierarchical quadruple system test \\
- triple    &: Simple hierarchical triple system test
\end{tabular}
\newline
\newline

\subsubsection{Creating your own initial conditions generators}
The subroutines that create the initial conditions are all contained in the \var{Ic} class.  This class is essentially a container for all subroutines and helper functions in one place.  If the user wishes to create their own initial conditions subroutines inside GANDALF, then there are 3 important files that must be editted.
\begin{itemize}
\item \var{src/Headers/IC.h} : this file contains the class definition for the \var{Ic} class.  The function prototype for any new IC generating function should be placed in this class, e.g. \\
\newline
\var{void MyNewInitialConditions(void);})
\item \var{src/Common/Ic.cpp} : this file contains all the functions contained in the Ic class, including the helper functions.  The full code for generating the initial conditions should be placed in this file, e.g. \\
\newline
\var{void Ic::MyNewInitialConditions(void) \{....\} }
\item \var{src/Common/SimulationIC.hpp} : Interface file that links (and calls) the IC class functions from the main simulation class.  In this file, you will need to add your own function call (together with the `if' statement to call the function), e.g. \\
\newline
\var{else if (ic == ``myic'') \{icGenerator.MyNewInitialConditions();\}}
\end{itemize}
To select your own initial conditions when running GANDALF, you simply set the \var{ic} parameter in your parameters file to the chosen string, e.g. \var{ic = myic}\,.


\subsection{Description of initial conditions}

\subsubsection{Accretion disc}

\label{sec:disc}

If ``ic'' is set to ``disc'', a standard thin accretion disc (e.g., \citealt{Pringle1981}) will be generated. See table \ref{table:disc} for a summary of the parameters. This setup can be used in 2-d (that is, vertically integrated) or in 3-d. Currently the setup only works in dimensionless units. A star with a mass of 1 is created at the origin of the coordinate system. The disc will extend in radius from ``DiscIcRin'' to ``DiscIcRout''. The disc has an initial power-law surface density profile $\Sigma \propto r^{-p}$, where the value of the exponent $p$ is given by the parameter ``DiscIcP''. The normalization is fixed by the parameter ``DiscIcMass'', which is the total mass of the disc (relative to the star). This setup is designed to use a locally isothermal equation of state, where the sound speed follows a power-law with radius: $c_s \propto r^{-q}$, where the exponent $q$ is given by the parameter ``DiscIcQ''. This equation of state is selected if ``gas\_eos'' has the value ``disc\_locally\_isothermal''. The setup emits a warning if one uses a different equation of state. When running in 3-d, the disc is assumed to be vertically isothermal and in hydrostatic equilibrium in the vertical direction, that is, $\rho \propto \exp (-z^2/2H^2)$, where $H=c_s/\Omega$ and $\Omega$ is the keplerian angular frequency. The normalization of the sound speed is set by the parameter ``DiscIcHr'', which sets the aspect ratio $H/r$ at the inner radius of the disc. Notice that we use random placement to initialise the disc; the initial density structure will thus be noisy. This in general not an issue since the particles typically rearrange themselves over one orbit. However, please keep this into account if you a need a structure with little noise already in the initial conditions.

The setup also supports the option of embedding a planet in the disc. This is selected if the parameter ``DiscIcPlanet'' takes a value of 1. The planet will be on a keplerian orbit; see table \ref{sec:disc} for setting the orbital parameters and planet mass. The setup does not currently support multiple planets, but it would be very straightforward to extend it - please contact the authors if you are interested in it.

As explained in chapter \ref{sec:dust}, GANDALF can solve for the dust dynamics in addition to the gas dynamics. If you wish to follow also the dust dynamics, the parameter ``NDust'' controls the number of dust particles and the parameter ``DustGasRatio'' gives the mass ratio between the dust and the gas. Notice that when using dust the parameter ``DiscIcMass'' refers to the \textit{total} (i.e., gas+dust) mass of the disc. Other choices of the dust algorithm are described in chapter \ref{sec:dust}.

Finally, the star and the planet will accrete mass from the disc if sink particles are switched on in the disc. This is controlled by the normal parameter ``sink\_particles''. The accretion radius of the star is always fixed to the inner radius of the disc; the accretion radius of the planet instead is controlled by the parameter ``DiscIcPlanetAccretionRadiusHill'', which sets the accretion radius in units of the Hill radius of the planet. Notice that, if the planet accretes, effectively its gravity is not smoothed, since we smooth the gravitational potential only inside the accretion radius. If you are not using sink particles, instead, the parameter ``DiscIcPlanetAccretionRadiusHill'' expresses the smoothing length of the gravitational potential. Currently is not possible to let only the star accrete and not the planet (or viceversa); please contact the authors if you require this behaviour.

Table \ref{table:disc} lists all the parameters discussed so far and their default values. Most of the default values come from the standard setup of \citet{deValBorro2006}; note however that we use a surface density $\propto r^{-1}$ rather than a flat one. We also do not employ ghost particles at the inner and outer boundary.

\begin{table}


\begin{tabular}{llc}
\textbf{Parameter name} & \textbf{Description} & \textbf{Default value}\\
DiscIcMass & Disc mass & 0.01\\
DiscIcP & Power-law exponent of the surface density & 1\\
DiscIcQ & Power-law exponent of the sound speed & 0.5 \\
DiscIcRin & Inner radius of the disc & 0.4\\
DiscIcRout & Outer radius of the disc & 2.5\\
DiscIcHr & Aspect ratio of the disc at the inner radius\\
DiscIcPlanet & Wheter the setup contains a planet (0 or 1) & 1\\
DiscIcPlanetRadius & Semi-major axis of the planet & 1\\
DiscIcPlanetEccen & Eccentricity of the orbit & 0\\
DiscIcPlanetIncl & Inclination of the orbit (in degrees) & 0\\
DiscIcPlanetMass & Mass of the planet & 1e-3\\
DiscIcPlanetAccretionRadiusHill & Accretion radius of the planet (in units of the Hill radius) & 0.4\\
NDust & Number of dust particles & 0\\
DustGasRatio & Dust to gas ratio when using dust & 0.01
\end{tabular}
\caption{Summary of the parameters controlling the disc setup and their default values}
\label{table:disc}
\end{table}





\subsection{Load from external file}
GANDALF can load two main file formats; simple ASCII column format and SEREN format.  These are described in detail in Section \ref{S:IO}.  In order to read from file, you must set the following parameters :
\begin{itemize}
\item \var{ic = file} : Tell GANDALF to ignore initial conditions generators and read from file instead
\item \var{in\_file = ..}$\;$ : The name of the initial conditions file
\item \var{in\_file\_form = column/sf/su} : The format of the initial conditions file
\end{itemize}

\subsection{Generate inside python script}
GANDALF can generate initial conditions if using python scripts.  The initial conditions can be first set in numpy arrays and then imported into the C++ code in order to perform the simulation.  More information is given about this method, including detailed examples, in Section \ref{SS:PYTHONTUTORIAL}.


\newpage



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{GANDALF Python library} \label{S:PYTHONSCRIPT}

The GANDALF Python library can be used to control a simulation via a Python script. In addition, it also provides routines for performing analysis tasks. Combined with the power of the Python language, this allows you to use it for many purposes. Some examples are listed below:
\begin{itemize}
\item Load in previously run simulations for analysis and producing visualisation
\item Load in multiple previously run simulations for comparisons
\item Prepare often-run simulations including the analysis in a single script
\item Run batches of simulations (e.g. a parameter study) with a single controlling script
\item Run or analyse simulations interactively through the interpreter
\item Generate initial conditions directly in python (instead of via a file or in C++) and run the simulation
\end{itemize}
It also provides an easier entry into using SPH and N-body than other codes which require more investigation of the code mechanics and file formats before even basic simulations can be run.

There are two main ways to use the Python side of GANDALF; (a) an interpreter, which works similarly to a shell (e.g. bash, csh) (b) a python script.  The interpreter understands a specific set of commands to load a simulation and plot quantities. You can use it for example to quickly read the output of a simulation and do a simple plot to check what is happening. Using a Python script instead is more powerful since, in addition to the same commands supported by the interpreter, you have the full power of Python at your disposal. For example you can access the raw particle data, compute additional particle properties, generate the initial conditions, and much more. The library will take care of the boring details such as reading the snapshots and extracting the data; you can concentrate on what to do with the data and on the science, and forget about the details.

The Python files for both the library and the interpreter are contained in the \var{gandalf/analysis} sub-directory. We explain briefly in the next sections how to run the interpreter or import the GANDALF library (if you are writing a script) and then an extensive tutorial follows which should clarify how to use the library.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Interpreter}

The interpreter is located in the analysis/ folder in the main directory of GANDALF.
To start the interpreter, type from the main GANDALF folder:\\
\newline
\noindent \var{python analysis/gandalf\_interpreter.py}. \\
\newline
A list of commands available for the python interpreter can be printed by typing the 'help' command.  Furthermore, typing 'help command' gives more information on the chosen command.  The commands have a streamlined syntax so that you can perform the tasks you need with as little typing as possible. For example, to render the density of the sph particles in a 2d simulation, you can just issue the command \var{render x y rho}. Under the hood, this is translated in the python code \lstinline[columns=fixed]{render('x','y','rho')}.


%Like in a shell, you can use the interpreter interactively or you can write little scripts for it (that is, texts files with a list of the commands that you want to be executed). \textbf{TODO: must actually implement reading from a script}. Any script can be read in and processed by typing : \var{read SCRIPTNAME}.

%The interpreter provides commands for creating new simulations, loading a simulation and do the plotting. The next section contains some examples to show how to use it.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Python script - overview of the libraries to import}
The GANDALF Python library contains a number of modules that can be imported in to provide the desired functionality.
\begin{itemize}
\item \lstinline[columns=fixed]{gandalf.analysis.facade} \\
\noindent This module is the main front-end to the GANDALF Python library and contains the Python-wrapped C++ executable and all functions required to set-up and run simulations.  This module must always be imported in GANDALF Python scripts.
\item \lstinline{gandalf.analysis.compute} \\
\noindent This module contains all extra and user-defined routines for computing important quantities from the simulation snapshot data.  Currently contains routines for computing centre-of-mass properties, L1 error norm (when provided with an analytical solution) and the Lagrangian radius.
\item \lstinline{gandalf.analysis.data_fetcher} \\
This module contains routines for exporting data from the C++ code to Python, generating custom data (e.g. time evolution) from snapshot data and generating user-defined data quantities from raw particle data.
\end{itemize}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\subsection{Write your Python script}

%All the commands available in the interpreter can also be used in your Python script. Just note the slight difference in syntax (since they are function, you need to put the arguments inside parentheses and they must be comma-separated). To use them, you need to import the module facade.py in the analysis/ folder. We recommend to add it to your PYTHONPATH to make the import simpler. Each function in facade comes with documentation, so you can inspect that to know what a function does.

%Here are the same examples from the previous section, rewritten to show how to use them in a Python script:
%\textbf{add examples}

%In addition, in your script you can also perform other tasks that are not possible from the interpreter, since you have a full programming language at your disposal. For example:
%\begin{itemize}

%\item You can modify from Python the parameters of the simulation before running it. This is very useful to do parameter studies, without generating in advance many different parameter files. Just remember to change the run\_id for each one of the simulations that you want to run, otherwise all the snapshots will be overwritten.

%\textbf{Example}

%\item It is possible to generate the initial conditions from Python. You need to create numpy arrays containing the quantities for the SPH particles (position, velocity, mass, internal energy) and you can import them in the code. An example follows:
%\textbf{Write the example!}

%\item You can define your own quantities to add to the plots in addition to the built-in ones. \textbf{Not implemented yet}

%\item It is possible to get access to the raw data of the simulation, in order to compute whatever quantity you need. \textbf{Actually, this is not possible yet in an easy way... I need to write a proper wrapper around the C++ method which is too low level. Remember -> give 'snap' attribute to the array to save reference counting. Implement it!}

%\item Everything else that you can think. The point is, it's GANDALF that takes care of reading the snapshots, getting the data, and all this boring parts; you can concentrate on what to do with the data and on the science, and forget about the details.

%\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\section{Structure of the python code}

%If you just want to use the python code, the previous section is enough. But if you want to develop some new feature, or you are just curious to know how we implemented it, you can go on reading.

%The functions defined in facade act mostly as wrappers around lower levels of abstraction. Sometimes they also make two different portions of the code communicate. The three main lower-level portions of the code are:
%\begin{itemize}

%\item SimBuffer. This class is responsible to keep track of all the simulations you have loaded, and the snapshots that are available for each one of them, either on disk or in memory (if you are running it). To eliminate the overhead of reading from the disk as much as possible (which is particularly important with network mounted filesystems), data is cached in memory once read. When all the cache is used, the buffer will look for objects that are no longer used since a long time and deallocate them to make space. The algorithm used is LRU (least recently used), which unfortunately is not scan resistant. However, we think that scans along all the snapshots are not something frequently enough to be a problem. The size in memory of the cache is a parameter that you can set; by default it is 1 GB. \textbf{Give to the user the possibility to change it}. Always use the functions in the buffer when you want to have the raw simulation/snapshot object, rather than reading them manually; in addition to saving you a lot of coding, this ensures that you are using caching and speeding up the reading.

%\item Plotting process. In order to have the figures responsive while the simulation is running, the plotting part is done in another process. A queue is used to make the main process (the one that responds to the user commands, and where the simulation runs) communicate with the plotting one. Two kinds of objects are transferred by this queue: the data to plot, and the commands to execute (see next element). The plotting process executes a loop: it reads what is there in the queue, and executes the commands if there is any. Then sleeps for a while, letting the user interact with the plots.

%There is a potential performance improvement in the transferring of the data along the queue. Internally, python uses the pickle module to serialize the array, which is quite slow. Using zeromq, it would be possible to copy the binary data, which would be just a memcopy and therefore very fast compared to the plotting itself. Another possible improvement is to store the data used, and avoid transferring it again each time a plot is updated.

%\item Commands. The functions in facade that interact with the plotting process create a command object and send them through the queue. A command is a class that contains a method that gets called by the plotting process when the object is received. Therefore, it's the command itself that knows what to do when it's called, rather than the plotting process. This separates the roles of the plotting process from the ones of the commands. Some of the commands also contain methods for getting the data from the simulation (this happens in the main process, where the data is stored).


%\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{GANDALF Python library tutorial} \label{SS:PYTHONTUTORIAL}
The GANDALF Python library contains a variety of commands that are used to create and run new simulations, load in old simulations from the disk, analyse and plot results and more.  A full list of all of these commands is given in Appendix \ref{reference_facade}.  The same information can be found in the source code file \var{gandalf/analysis/facade.py}.  Also, when running the interpreter, information on interactive commands can be obtained by typing \var{help} for a full list of available commands, or \var{help command} for information on that particular command. Finally, as common in Python, if running through a script one can get the documentation string of a function in the special variable called \lstinline{__doc__}.

We provide here a short tutorial demonstrating from the most basic to the more advanced functionality of the Python library.  These examples are also contained in the \var{gandalf/examples} sub-directory and should run as typed, although it is also useful for the user to write the examples themselves if unfamiliar with python syntax. \\
\newline


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example 1 - Creating and running a new simulation from a parameter file}
This first example demonstrates how to set-up and run a new simulation from a given parameters file (nominally called \var{params.dat} in these examples).  \\

%\begin{tabular}{p{14.0cm}}
%\var{from gandalf.analysis.facade import *} \\
%\var{sim = newsim(`adsod.dat')}            \\
%\var{setupsim()} \\
%\var{run()} \\
%\end{tabular} \\
%\newline

\lstinputlisting{../examples/example01.py}

\noindent The first command \singlecommand{from gandalf.analysis.facade import *} loads in all definitions and functions from the GANDALF python frontend, \var{facade.py}.  Note that this assumes that the \var{PYTHONPATH} environment variable has been correctly set-up.  Otherwise, the full absolute path of the \var{facade.py} file must be given.  The second command \singlecommand{sim = newsim(`adsod.dat')} reads in the parameters file \var{'adsod.dat'} and creates a new simulation from those parameters and returns a Python object, \var{sim}, which can be used to refer to that simulation in Python (N.B. the user can use any name for this object, as long as it does not clash with another GANDALF or Python object).  In this case, it is assumed the parameters file is selecting an initial conditions generator in the main code to generate the particles in the simulation.  In order to fully set-up the simulation so it is ready to be run, the command \singlecommand{setupsim()} must be run.  Finally, to run the simulation to completion, we must execute the command \singlecommand{run()}  The code will run until the specified endtime and then exit the python environment.

If run through the interpreter, the same example would look like:\\
\newline
\begin{tabular}{p{14.0cm}}
\var{newsim adsod.dat}            \\
\var{setupsim} \\
\var{run} \\
\end{tabular} \\
\newline




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example 2 - Creating a simulation and modifying a parameter before running}
This example demonstrates how to modify parameters in Python that have been loaded in from a given parameters file. \\

%\begin{tabular}{p{14.0cm}}
%\var{from gandalf.analysis.facade import *} \\
%\var{sim = newsim(`adsod.dat')} \\
%\var{sim.SetParam(`tend',1.5)} \\
%\var{setupsim()} \\
%\var{run()} \\
%\end{tabular} \\

\lstinputlisting{../examples/example02.py}



\noindent The first two lines of this example are the as the first example, where we import the GANDALF Python frontend and then create a new simulation from the parameters file \var{`adsod.dat'}.  Once loaded in, we can now modify parameters with the command \singlecommand{sim.SetParam(parameter\_name,new\_value)}  In this case, we set the parameter \var{tend} (the total simulation run-time) to the value \var{2.0} with \lstinline{sim.SetParam(`tend',2.0)}.  Once all parameters have been modified, then we can set-up the simulation fully with \lstinline{setupsim()} and run it with \lstinline{run()}.  Note that if we attempt to modify any parameters AFTER calling \lstinline{setupsim()}, then the code will terminate with an exception. It is not possible to change the parameters of a simulation through the interpreter.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example 3 - Generating initial conditions using internal routines without a parameters file}
In this example, we will create and run a simulation without needing to load-in an external parameters file. \\

%\begin{tabular}{p{14.0cm}}
%\var{from gandalf.analysis.facade import *} \\
%\var{sim = newsim(ndim=2,sim=`sph')} \\
%\var{sim.SetParam(`sph',`gradh')} \\
%\var{sim.SetParam(`ic',`sedov')} \\
%\var{sim.SetParam(`run\_id',`SEDOV1')} \\
%\var{sim.SetParam(`Nlattice[0]',64)} \\
%\var{sim.SetParam(`Nlattice[1]',64)} \\
%\var{sim.SetParam(`boxmin[0]',-1.0)} \\
%\var{sim.SetParam(`boxmin[1]',-1.0)} \\
%\var{sim.SetParam(`boxmax[0]',1.0)} \\
%\var{sim.SetParam(`boxmax[1]',1.0)} \\
%\var{sim.SetParam(`dimensionless',1)} \\
%\var{sim.SetParam(`Nlevels',10)} \\
%\var{sim.SetParam(`tend',0.5)} \\
%\var{sim.SetParam(`tsnapfirst',0.0)} \\
%\var{sim.SetParam(`dt\_snap',0.1)} \\
%\var{setupsim()} \\
%\var{run()} \\
%\end{tabular} \\

\lstinputlisting{../examples/example03.py}



\noindent In this example, we call the \var{newsim} function passing the dimensionality of the simulation instead of a parameter file.  Therefore the line \lstinline{sim = newsim(ndim=2,sim=`sph')} sets up a 2-dimensional SPH simulation, but is otherwise unspecified because no other parameters have been set.  Next, we can set all the parameters that would have otherwise been set in the parameters file.  In this example, we set up an SPH simulation with initial conditions to perform the Sedov blast-wave test (\lstinline{sim.SetParam(`ic',`sedov')}); (See Section \ref{S:PARAMETERS} for a full description of all parameters).  Once all parameters have been set, we can fully set-up and run the simulation with \lstinline{setupsim()} and \lstinline{run()}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example 4 - Running a simulation and plotting results to screen and file}
In this example, we show how to plot particle data interactively during a live simulation. \\

%\begin{tabular}{p{14.0cm}}
%\var{from gandalf.analysis.facade import *} \\
%\var{sim = newsim(`glass.dat')} \\
%\var{sim.SetParam(`dt\_python',2.0)} \\
%\var{setupsim()} \\
%\var{plot(`x',`y')} \\
%\var{run()} \\
%\var{savefig(`figure.eps')} \\
%\var{block()} \\
%\end{tabular} \\

\lstinputlisting{../examples/example04.py}



\noindent In this simulation, we will run a live simulation while also plotting particle data in a matplotlib window as the simulation is still running.  We first create the simulation as in previous examples with \var{newsim}.  We then set the value of the \var{dt\_python} which is the time between python/matplotlib window updates.  After the \var{setupsim()} command, we can then plot the initial conditions from the simulation using the command \singlecommand{plot(x,y)} where \var{x} and \var{y} are the x- and y-axis plotting variables (N.B. these can be any variables, not necessarily x and y. e.g. \var{plot(`y',`z')}).  If we now run the simulation with \var{run()}, the simulation will compute as normal while updating the plot window every \var{dt\_python} seconds (in this casem every 2 seconds).  After the simulation has finished, we can save the final plot to file with the command \singlecommand{savefig(filename)} where \var{filename} is the intended filename {\it including} the file extension.  matplotlib automatically determines the correct file format from the extension  (e.g. typing \lstinline{savefig(`fig1.png')} will automatically save the window to a png file without need for additional arguments). The use of \var{block()} in the last line is needed because otherwise python would exit immediately from the script. Instead, the function pauses the script so that you can look at the plot. Pressing Enter (or any other key) causes to script to continue (and so in this case to terminate as there is no other instruction afterwards).

This same example can also be run through the interpreter (apart for changing the \var{dt\_python} parameter). Note that there is no need to call block as the interpreter already pauses after each command in order to wait for the next one:\\
\newline
\begin{tabular}{p{14.0cm}}
\var{newsim glass.dat} \\
\var{setupsim} \\
\var{plot x y} \\
\var{run} \\
\var{savefig figure.eps } \\
\end{tabular} \\




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example 5 - Reading a simulation from disc and plotting several snapshots}
In this example, we load in a previously run simulation from the disk and generate various plots. \\

%\begin{tabular}{p{14.0cm}}
%\var{from gandalf.analysis.facade import *} \\
%\var{loadsim('ADSOD1')} \\
%\var{plot(`x',`rho',snap=0)} \\
%\var{savefig(`snap1.eps')} \\
%\var{plot(`x',`rho',snap=1)} \\
%\var{savefig(`snap2.eps')} \\
%\var{plot(`x',`rho',snap=0)} \\
%\var{addplot(`x',`rho',snap=1)} \\
%\var{savefig(`snap12.eps')} \\
%\var{block()} \\
%\end{tabular} \\

\lstinputlisting{../examples/example05.py}



\noindent Instead of a creating a new simulation with the \var{newsim} function, we will load in a previous run simulation with the command \singlecommand{loadsim(runid)} where \var{runid} is the simulation run identification string.  Note, this assumes that all files from that simulation are present in the folder, including the \var{runid.param} file which contains all the parameters used to run the simulation.  In this case, we load a simulation with the runid \var{ADSOD1} with the command \lstinline{loadsim(`ADSOD1')}.  We can then plot particular snapshots of the simulation using the \var{plot} command with the optional \var{snap} argument.  Therefore to plot the first snapshot, type \var{plot(`x',`y',snap=0)}, to plot the second snapshot, type \var{plot(`x','y',snap=1)}, etc..  Note that Python uses C-style indexing, i.e. starting from zero, unlike other conventions like in FORTRAN where indexing starts from one.  If we wish to overplot one snapshot over another, we first plot the first snapshot and then use the command \singlecommand{addplot(`x',`y',snap=1)} to add the second snapshot plot over the first one.  This behaviour can also be replicated by adding the \var{overplot} optional boolean argument to \var{plot}, i.e. \lstinline{plot(`x',`y',snap=1,overplot=True)}.

It is possible to run this example too through the interpreter:\\
\newline
\begin{tabular}{p{14.0cm}}
\var{loadsim ADSOD1 } \\
\var{plot x rho snap=0} \\
\var{savefig snap1.eps} \\
\var{plot x rho snap=1} \\
\var{savefig snap2.eps} \\
\var{plot x rho snap=0} \\
\var{addplot x rho snap=1} \\
\var{savefig snap12.eps} \\
\end{tabular} \\
\newline
From now on, we stop providing also the code for the interpreter. As you have seen, the conversion is trivial to do.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example 6 - Reading  and plotting multiple simulations}
In this example, we show how to read in multiple simulations to make comparison plots. \\

%\begin{tabular}{p{14.0cm}}
%\var{from gandalf.analysis.facade import *} \\
%\var{sim0 = loadsim(`ADSOD1')} \\
%\var{sim1 = newsim(`adsod.dat')} \\
%\var{sim1.SetParam(`run\_id',`ADSOD2')} \\
%\var{sim1.SetParam(`avisc',`none')} \\
%\var{setupsim()} \\
%\var{run()} \\
%\var{plot(`x',`rho',sim=0,snap=2)} \\
%\var{addplot(`x',`rho',sim=1,snap=2)} \\
%\var{savefig(`sod12.eps')} \\
%\var{block()} \\
%\end{tabular} \\

\lstinputlisting{../examples/example06.py}




\noindent Multiple simulations can be loaded into GANDALF using the \var{loadsim} command.  However, they must be returned to objects with different names; in this case, we name them \var{sim0} and \var{sim1} (to follow the C-style numbering convention).  A snapshot from a particular simulation can be plotted by adding the optional argument \var{sim} to the plot command, i.e. \lstinline{plot(`x',`rho',sim=0,snap=0)} to plot the first snapshot in the first simulation.  The second simulation is then over-plotted with the command \lstinline{addplot(`x',`rho',sim=1,snap=0)}.  In this example, we load in an old simulation (as run in example 1) and then run a new simulation (the same initial conditions but without any artificial viscosity) and then plot the results on the same window and to an eps file.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example  7 - Overplotting the analytical solution with the simulation results}
In this example, we demonstrate how to overplot analytical solutions over particle plots. \\

%\begin{tabular}{p{14.0cm}}
%\var{from gandalf.analysis.facade import *} \\
%\var{sim = newsim(`adsod.dat')} \\
%\var{sim.SetParam(`dt\_python',2.0)} \\
%\var{setupsim()} \\
%\var{plot(`x',`rho')} \\
%\var{plotanalytical(`x',`rho')} \\
%\var{limit(`x',-1.1,1.1)} \\
%\var{run()} \\
%\var{block()} \\
%\end{tabular} \\

\lstinputlisting{../examples/example07.py}



\noindent For some sets of initial conditions (e.g. shocktubes), the analytical solution is provided in one of the python modules (\var{gandalf.analysis.analytical}).  The python module will read in the values specified in the simulation parameters file and return the correct solution for that given snapshot time.  To overplot a specified quantity (e.g. \var{x} vs density), then use the command \singlecommand{plotanalytical(`x',`rho')}  Another useful command to set the limits of any plot quantity is \singlecommand{limit(quantity,min,max)} where \var{quantity} is a string of the quantity and \var{min} and \var{max} define the range of the plot.  In this example, \lstinline{limit(`x',-1.1,1.1)}
%\newpage


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example 8 - Creating initial conditions directly in the python script}
In this example, we demonstrate how to generate some simple initial conditions inside a python script.  For simplicity, we start with a simple 1D shocktube example. \\

%\begin{tabular}{p{14.0cm}}
%\var{from gandalf.analysis.facade import *} \\
%\var{import numpy as np} \\
%\var{import time} \\
%\var{Nsph = 200} \\
%\var{vfluid = 4.0} \\
%\var{xmin = -1.5} \\
%\var{xmax = 1.5} \\
%\var{deltax = (xmax - xmin) / Nsph} \\
%\var{x = np.linspace(xmin + 0.5*deltax,xmax - 0.5*deltax,num=Nsph)} \\
%\var{m = np.ones(Nsph)*(xmax - xmin)/Nsph} \\
%\var{vx = np.ones(Nsph)*vfluid}  \\
%\var{vx[x > 0] = -vfluid} \\
%\var{sim = newsim(ndim=1,sim=`sph')} \\
%\var{sim.SetParam(`ic',`python')} \\
%\var{sim.SetParam(`gas\_eos',`isothermal')}  \\
%\var{sim.SetParam(`Nsph',200)}  \\
%\var{sim.SetParam(`tend',0.2)}  \\
%\var{sim.SetParam(`dt\_snap',0.05)} \\
%\var{sim.SetParam(`dimensionless',1)} \\
%\var{sim.SetParam(`vfluid1[0]',vfluid)} \\
%\var{sim.SetParam(`vfluid2[0]',-vfluid)} \\
%\var{sim.SetParam(`boxmin[0]',xmin)} \\
%\var{sim.SetParam(`boxmax[0]',xmax)} \\
%\var{sim.SetParam(`run\_id',`SHOCKTUBE1')} \\
%\var{sim.PreSetupForPython()} \\
%\var{sim.ImportArray(x,`x')} \\
%\var{sim.ImportArray(vx,`vx')} \\
%\var{sim.ImportArray(m,`m')}   \\
%\var{sim.SetupSimulation()}  \\
%\var{plot(`x',`rho')} \\
%\var{plotanalytical(`x',`rho',ic=`shocktube')} \\
%\var{limit(`x',-0.17,0.17,window=`all')} \\
%\var{limit(`rho',0.0,21.0,window=`all')} \\
%\var{run()} \\
%\var{savefig(`shocktube.png')} \\
%\var{block()} \\
%\end{tabular} \\

\lstinputlisting{../examples/example08.py}



\noindent As well as importing the GANDALF Python library, we must also import the NUMPY library, which is a popular mathematical library in Python for handling arrays, with the command \singlecommand{import numpy as np}
We first create the initial conditions in Python before exporting the data to the C++ part of the code.  This can be done by creating NUMPY arrays for each variable, e.g. \var{x}, \var{m}, etc..  Note that even for vector quantities (e.g. position, velocity), we must allocate separate 1D NUMPY arrays for each component, for example the commands \\
\newline
\indent \lstinline{a1 = np.zeros(N)}, \\
\indent \lstinline{a2 = np.ones(N)}, \\
\indent \lstinline{a3 = np.linspace(0.0,1.0,num=N)} \\
\newline
return NUMPY arrays of size \var{N} with all zeros, all ones, and uniform values between 0 and 1 respectively  (See NUMPY website and documentation for more information on NUMPY functions).  For this simple shocktube test, we set all velocities left of the shock interface to $4.0$ and right of the shock interface to $-4.0$ with the commands \lstinline{vx = np.ones(Nsph)*4.0; vx[x > 0.0] = -4.0}.

We then create an undefined 1D simulation object with \lstinline{sim = newsim(ndim=1)} and set various parameters determing the simulation type, integration scheme, equation of state, etc..  We must then inform the code that we intend to generate initial conditions in the python script (and not say from an external file or one of the internal C++ routines) by the command \lstinline{sim.SetParam('ic','python')}.  We must also inform the code of the number of SPH particles to be used with the command \lstinline{sim.SetParam('Nsph',Nsph)}, in order to allocate enough memory for the simulation (N.B. \var{Nsph} is simply a local Python variable in this example.  This can be any other variable, or even simply an integer number).

The last step before importing the arrays to C++, we must initialise various things in the code, including importantly allocating memory for the particles, by calling the function \singlecommand{sim.PreSetupforPython()}
Next, we can import the array into the C++ arrays by the command \singlecommand{sim.ImportArray(numpyarray,varname)}
where \var{numpyarray} is the local NUMPY array and \var{varname} is a string containing the C++ variable name.  For example, to import the x-positions, we call the command \lstinline{sim.ImportArray(x,'x')}.  Once all arrays have been called, we can finally call the \var{setupsim} function and plot and run the simulation to completion.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example  9 - Creating initial conditions for N-body simulation in python script}
In this example, we demonstrate how to generate some simple initial conditions for pure N-body simulation inside a python script. \\

%\begin{tabular}{p{14.0cm}}
%\var{from gandalf.analysis.facade import *} \\
%\var{import numpy as np} \\
%\var{import time} \\
%\var{Nstar = 3} \\
%\var{x = np.zeros(Nstar)} \\
%\var{y = np.zeros(Nstar)} \\
%\var{vx = np.zeros(Nstar)} \\
%\var{vy = np.zeros(Nstar)} \\
%\var{m = np.zeros(Nstar)} \\
%\var{h = 0.000001*np.ones(Nstar)} \\
%\var{m[0] = 3.0;}\,\,\, \var{x[0] = 1.0;}\,\,\,\var{y[0] = 3.0} \\
%\var{m[0] = 4.0;}\,\,\, \var{x[0] = -2.0;}\,\,\,\var{y[0] = -1.0} \\
%\var{m[0] = 5.0;}\,\,\, \var{x[0] = 1.0;}\,\,\,\var{y[0] = -1.0} \\
%\var{sim = newsim(ndim=2,sim=`nbody')} \\
%\var{sim.SetParam(`ic',`python')} \\
%\var{sim.SetParam(`nbody',`hermite4ts')} \\
%\var{sim.SetParam(`sub\_systems',0)}  \\
%\var{sim.SetParam(`Npec',3)} \\
%\var{sim.SetParam(`Nlevels',1)} \\
%\var{sim.SetParam(`tend',80.0)} \\
%\var{sim.SetParam(`dt\_snap',1.0)} \\
%\var{sim.SetParam(`noutputstep',128)} \\
%\var{sim.SetParam(`ndiagstep',2048)} \\
%\var{sim.SetParam(`dimensionless',1)} \\
%\var{sim.SetParam(`run\_id',`BURRAU1')} \\
%\var{sim.PreSetupForPython()} \\
%\var{sim.ImportArray(x,`x',`star')} \\
%\var{sim.ImportArray(y,`y',`star')} \\
%\var{sim.ImportArray(vx,`vx',`star')} \\
%\var{sim.ImportArray(vx,`vy',`star')} \\
%\var{sim.ImportArray(m,`m',`star')}   \\
%\var{sim.ImportArray(h,`h',`star')}   \\
%\var{sim.SetupSimulation()}  \\
%\var{plot(`x',`y',type=`star')} \\
%\var{limit(`x',-30.0,30.0,window=`all')} \\
%\var{limit(`y',-20.0,40.0,window=`all')} \\
%\var{run()} \\
%\var{block()} \\
%\end{tabular} \\

\lstinputlisting{../examples/example09.py}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example  10 - Generating rendered images from SPH simulations}
In this example, we show how to generate rendered images of some quantity.  Note that rendered plots can only be created for simulations using 2 or 3 dimensions. \\

%\begin{tabular}{p{14.0cm}}
%\var{from gandalf.analysis.facade import *} \\
%\var{sim = newsim(`khi.dat')} \\
%\var{setupsim()} \\
%\var{render(`x',`y','rho',res=128)} \\
%\var{limit(`x',-0.5,0.5)} \\
%\var{limit(`y',-0.5,0.5)} \\
%\var{limit(`rho',1.0,2.0)} \\
%\var{run()} \\
%\var{block()} \\
%\end{tabular} \\

\lstinputlisting{../examples/example10.py}



\noindent A rendered plot of the density on the x-y plane can be created using the command \singlecommand{render(`x',`y',`rho')}  The first argument specifies the x-coordinate, the second argument the y-coordinate and the third argument the rendered quantity.  An important optional argument is \var{res} which specifies the resolution (in pixels) of the rendered image in each dimension.  If the image is not square, then the user may pass the resolution in each dimension in parenthesis, e.g. \lstinline{render(`x',`y','rho',res=(128,64))}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example  11 - Plotting in alternative coordinate systems}
In this example, we show how to plot in alternative coordinate systems (i.e. other than Cartesian coordinates). \\

%\begin{tabular}{p{14.0cm}}
%\var{from gandalf.analysis.facade import *} \\
%\var{sim = loadsim(`SEDOV1')} \\
%\var{plot(`r',`rho')} \\
%\var{block()} \\
%\end{tabular} \\

\lstinputlisting{../examples/example11.py}




\noindent Once the simulation has been loaded or set-up as usual, then it is possible to plot in 3 pre-defined coordinate systems, Cartesian, spherical polar and cylindrical polar coordinates, or even a mixture.  Cartesian coordinates are represented by the strings, \var{x}, \var{y}, \var{z}; sphericial polar coordinates by the strings \var{r}, \var{theta}, \var{phi}; cylindrical coordinates by the strings \var{R}, \var{phi}, \var{z}.  For example, to plot the density as a function of radius (for spherically symmetric simulations), then we can use the command \singlecommand{plot(`r',`rho')}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example  12 - Changing the plotting units}
In this example, we show how to plot using different units to those used in the simulation or provided in the snapshot files. \\

%\begin{tabular}{p{14.0cm}}
%\var{from gandalf.analysis.facade import *} \\
%\var{sim = newsim(`bossbodenheimer.dat')} \\
%\var{sim.SetParam(`tend',0.02)} \\
%\var{setupsim()} \\
%\var{plot(`x',`y')} \\
%\var{run()} \\
%\var{block()} \\
%\var{plot(`x',`y',xunit='au',yunit='au')} \\
%\end{tabular} \\

\lstinputlisting{../examples/example12.py}


\noindent In order to plot quantities using a different unit to the default, the optional arguments \var{xunit} and \var{yunit} must be appended to the argument list for the relevant plotting command.  These arguments contain strings of the required unit (See section \ref{SS:UNITPARAMS} for the list of available units).


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example  13 - Creating and plotting user-defined quantities}
In this example, we show how to create new quantities to be plotted as regular quantities by the python library. \\

%\begin{tabular}{p{14.0cm}}
%\var{from gandalf.analysis.facade import *} \\
%\var{sim = loadsim(`TEST1')} \\
%\var{CreateUserQuantity(`ke',`0.5*m*vx*vx',scaling\_factor=`u',label=`frac\{1\}\{2\}mv\^2')} \\
%\var{plot(`x',`ke',snap=2)} \\
%\var{limit(`ke',0.0,2.7)} \\
%\var{addplot(`x',`u',snap=2)} \\
%\var{block()} \\
%\end{tabular} \\

\lstinputlisting{../examples/example13.py}



\noindent We can create new quantities from existing particle properties using the command \singlecommand{CreateUserQuantity(quantity\_name,quantity\_formula)} where \var{quantity\_name} is the new string name of the derived quantity and \var{quantity\_formula} is string containing an algebraic formula that defines the new quantity.  In our example, we create a new quantity to compute the kinetic energy of a particle with the algebraic formula \var{0.5*m*vx*vx}.  For 2D and 3D simulations, this would become \var{0.5*m*(vx*vx + vy*vy)} and \var{0.5*m*(vx*vx + vy*vy + vz*vz)} respectively.  Note that the algebraic formula can also use parameters by using the same string as defined in the parameters file.

The new quantity can then be plotted as any other regular quantity, e.g. \var{plot(`x',`ke')}.  Two important optional arguments are \var{scaling\_factor}, which defines the string of the unit (e.g. in this case \var{E} since \var{ke} is an energy), and \var{label}, which is a latex string used on the plots (e.g. \var{frac\{1\}\{2\}mv\^2}).



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example  14 - Plotting time series of particle properties}
In this example, we show how to plot a time series of particle quantities from each snapshot recorded in the simulation buffer.  This can be used to plot useful quantities versus time (e.g. mass of a sink versus time), or x-y tracks of a particle as it moves in space. \\

%\begin{tabular}{p{14.0cm}}
%\var{from gandalf.analysis.facade import *} \\
%\var{sim = loadsim(`TEST1')} \\
%\var{time\_plot(`t',`rho',id=0,linestyle=`-')} \\
%\var{block()} \\
%\var{time\_plot(`x',`y',id=0,linestyle=`-')} \\
%\var{block()} \\
%\end{tabular} \\

\lstinputlisting{../examples/example14.py}



\noindent In order to plot a time series of some particle property, we use the command \singlecommand{time\_plot(`t',`rho',id=0)} where \var{id} is the unique/original id of the particle to be plotted.  If no id is given, then an error is returned and nothing will be plotted on screen.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\subsubsection{Example  16 - Creating and plotting multiple sub-plots on a single window}
%To be done



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example  15 - Creating an animation from simulation snapshot files}
In this example, we show how to generate a movie from a series of snapshot files loaded into memory. \\


%\begin{tabular}{p{14.0cm}}
%\var{from gandalf.analysis.facade import *} \\
%\var{sim = loadsim(`TEST1')} \\
%\var{plot(`x',`y')} \\
%\var{make\_movie(`MOVIE1.mp4')} \\
%\var{block()} \\
%\end{tabular} \\
%\newline

\lstinputlisting{../examples/example15.py}




\noindent Before generating a movie, we must first specify the plot commands that are to be used to generate the movie, for example with the regular plot command \lstinline{plot('x','rho')}.  Once the commands to generate the plot window have been completed, we then generate the movie with the command \singlecommand{make\_movie(`MOVIE1.mp4')} where the only required argument is the filename of the movie.  The python library will then step through each snapshot in turn, generate the image as a temporary png file, and then attach all images together using ffmpeg.  If you do not have ffmpeg installed, the command will fail. All temporary png files should then be deleted automatically from the disc.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example  16 - Retrieving data from the simulation}
In this example, we show how to save some data from the snapshots in a variable. You can even save the rendered image. \\



\lstinputlisting{../examples/example16.py}




\noindent The syntax of the functions get\_data and get\_render\_data should be self-explanatory. These routines just return the requested quantity as a numpy array, that you can save in a variable for further processing. Alternatively, every time you do a plot you can also save the return value in a variable as shown further on in the example. Particle plots save the data used on the x-axis inside the \lstinline{x_data} field and on the y-axis inside the \lstinline{y_data} field. Render instead saves the image in the \lstinline{render_data} field. In this way, you can use the rendered image to grid your SPH data. At the moment render only renders to a 2d grid; however, we plan in the future to extend its capabilities to render to a 3d grid.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example  17 - Creating and plotting user-defined quantities from a function given by the user}
In this example, we show how to create and plot an user-defined quantity. However, rather than building it from a mathematical formula as shown in the example 13, we compute it using a function provided by the user. This allows you to have much more flexibility in performing the computations you want with the data of the simulation. \\



\lstinputlisting{../examples/example17.py}




\noindent Refer to the example 13 for the first use of \lstinline{CreateUserQuantity}. For what concerns the second, we first need to define the function that we want to use for computing the data. In this case, the function just retrieves the array with the velocities (see the previous example) and computes the kinetic energy from them. The function needs to return the computed value. Note that we pass the parameters specifying the snapshot, the type of particle (e.g., gas or star) and the unit to \lstinline{get_data}. This is very important to do; if you forget to do it, you might get data from the wrong snapshot, the wrong particle type, or in the wrong units.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Example  18 - Creating and plotting time series from a function given by the user}
In this example, we show how to plot time series from functions defined by the user.  \\



\lstinputlisting{../examples/example18.py}




\noindent For accomplishing this task we need to use the function \lstinline{CreateTimeData}. Note that the way this function works is very similar to \lstinline{CreateUserQuantity} shown in the previous example. In the first example we use the function \lstinline{lagrangian_radii}, which is defined in GANDALF in the module \lstinline{compute}. Note that we can pass additional parameters to the function using keywords; for example in this case we specify that we want the half-mass radius using the keyword \lstinline{mfrac}. In the second example, we show how to define your own function for computing a quantity from a snapshot. Just remember that in this case (differently from the previous example) you need to return a scalar, which will be plotted as a function of time. The example retrieves the array with the masses (see example 16) and sums up all of them to get the total one.


\subsection{Tips and tricks}

To plot the results we use the well-known library matplotlib. This means that if you are running from a script you have several ways to customize the plots if you are not satisfied by their graphical appearance:
\begin{itemize}
\item Calling functions in matplotlib. Example: to logscale the y axis, just call \lstinline{plt.yscale('log')} (assuming you have imported matplotlib).
\item Passing keyword arguments to the plotting functions defined by GANDALF. Example: if you are doing a particle plot and want to change the marker size, you can call \lstinline{plot('x','y',ms=2)}, and the ms keyword will be passed to matplotlib.
\item Modifying the defaults for matplotlib using the .matplotlibrc files (see matplotlib documentation at \url{http://matplotlib.org/users/customizing.html}).
\end{itemize}

You might notice that the matplotlib window freezes when you are running a simulation live. This is not GANDALF's fault unfortunately, but a flaw in the design of the Python interpreter (if you want to know more, read here: \url{https://wiki.python.org/moin/GlobalInterpreterLock}). If you do require interaction with the matplotlib window (for example to be able to pan/zoom while the simulation is running), you can run the plots in another process. In order to do this, open the file defaults.py in the analysis folder and change the value of \var{parallel} to \lstinline{True}. There is a drawback in this case: given that the plots are on another process, you cannot call matplotlib functions to customise your plots (but the other two ways we suggest for customising the plots still work).

\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Dust Dynamics}

\label{sec:dust}

\subsection{Theory}
The equations of motion for the mutual dynamics of gas and dust coupled by drag forces are
\begin{align}
\rho_g \frac{\mathrm{D}_g \vec{v}_g}{\mathrm{D} t} &= -\nabla P  + \rho_g \vec{a}_{g, \mathrm{ext}}
                                                     + K (\vec{v}_d - \vec{v}_g),  \\
\rho_d \frac{\mathrm{D}_d \vec{v}_d}{\mathrm{D} t} &= \rho_d \vec{a}_{d, \mathrm{ext}}
                                                     - K (\vec{v}_d - \vec{v}_g),
\end{align}
where $\frac{\mathrm{D}_g}{\mathrm{D} t}$ is the Lagrangian derivative with respect to the gas. The external accelerations are included to take into account additional forces, which may either be the same for both the gas and dust particles (such as gravity) or affect just one species (e.g. viscosity). These equations can be re-written in terms of an equation of motion for the relative and barocentric velocities, $\vec{v} = (\rho_g \vec{v}_g + \rho_d \vec{v}_d)/(\rho_g + \rho_d)$ and $\Delta \vec{v} = \vec{v}_d - \vec{v}_g$, (Youdin \& Goodmann, 2005; Laibe \& Price, 2014; Lor\'en-Aguilar \& Bate, 2015),
\begin{align}
\frac{\mathrm{D} \vec{v}}{\mathrm{D} t} &= \frac{-\nabla P + \rho_g \vec{a}_{g, \mathrm{ext}} + \rho_d \vec{a}_{d, \mathrm{ext}}}{\rho_g + \rho_d} - \vec{F}, \\
\frac{\mathrm{D} \Delta \vec{v}}{\mathrm{D} t} &= -  \frac{\Delta \vec{v}}{t_s}+  \vec{a}_{d, \mathrm{ext}} + \frac{\nabla P}{\rho_g}  - \vec{a}_{g, \mathrm{ext}} - (\Delta \vec{v} \cdot \nabla) \vec{v} - \vec{G},
\end{align}
where $t_s = \rho_g \rho_d / K(\rho_g + \rho_d)$. The functions $\vec{F}$ and $\vec{G}$ are
\begin{align}
\vec{F} &= \frac{1}{\rho_d + \rho_g} \nabla \cdot \left(\frac{\rho_d \rho_d}{\rho_d + \rho_d} \Delta \vec{v} \Delta \vec{v} \right), \\
\vec{G} &= \frac{1}{\rho_d + \rho_g} \left[
 \rho_g \Delta \vec{v} \cdot \nabla \left(\frac{\rho_g}{\rho_g + \rho_d} \Delta \vec{v} \right) -
 \rho_d \Delta \vec{v} \cdot \nabla \left(\frac{\rho_d}{\rho_g + \rho_d} \Delta \vec{v} \right) \right].
\end{align}
These two formulations naturally lead to two different numerical approaches (in the same spirit as Eulerian vs. Lagrangian hydro codes), with equations (1) and (2) naturally represented by two different particle types representing fluid (gas) particles and dust particles. Equations (3) -- (6) are more naturally represented as a single fluid that moves with the barycentric velocity, on top of which the dust fraction is advected relative to the barycentre at velocity $\Delta \vec{v}$ (see Laibe \& Price 2014 for more details).

In GANDALF we solve directly equations (1) and (2) using a multi-fluid approach, using the semi-implicit time-stepping following Booth, Sijacki \& Clarke (2015) and Lor\'en-Aguilar \& Bate (2015). In these approaches we use of $\langle \Delta \vec{v} \rangle$, the average of the relative velocity in the time-step $\left[ t,t + \Delta t\right]$. This is found by solving equation (4) under the approximation that the densities, accelerations and $t_s$ are constant in both time and space. We find
\begin{equation}
\langle \Delta \vec{v} \rangle \frac{\Delta t}{t_s} = \Delta \vec{v}(t+\Delta t) \left(1 - \exp(-\Delta t/ t_s\right) -
\left(\vec{a}_{d, \mathrm{ext}} + \frac{\nabla P}{\rho_g}  - \vec{a}_{g, \mathrm{ext}}\right)\left[(\Delta t + t_s)\left(1 - \exp(-\Delta t/ t_s)\right) - \Delta t \right],
\end{equation}
where $ \Delta \vec{v}(t+\Delta t) = \Delta \vec{v} (t) + \left(\vec{a}_{d, \mathrm{ext}} + \nabla P/\rho_g  - \vec{a}_{g, \mathrm{ext}}\right) \Delta t$. From equation (7) the gas and dust accelerations can be found to be
\begin{align}
\vec{a}_{g,\mathrm{drag}} &= + \frac{\rho_d}{\rho_d + \rho_g} \frac{\langle \Delta \vec{v} \rangle}{t_s}, \\
\vec{a}_{d,\mathrm{drag}} &= - \frac{\rho_g}{\rho_d + \rho_g} \frac{\langle \Delta \vec{v} \rangle}{t_s}.
\end{align}
When the approximations made above hold exactly the resulting accelerations are exact since equations (1) and (2) are linear as long as the drag coefficient does not depend on $\Delta v$. The change in internal energy is computed directly from the change in kinetic energy to ensure conservation.

\subsection{Implementations}

GANDALF now includes two implementations of combined dust-gas dynamics, based upon equations (1) \& (2) and  (7) -- (9). Currently these are only implemented for the grad-h SPH algorithm. In both implementations the dust and gas are represented by separate particle types (which are identified in the code via the \lstinline{ptype} particle attribute). In the first implementation the back-reaction of the dust on the gas is neglected (equivalent to  $\rho_d \rightarrow 0$) -- the gas dynamics are not modified when dust is included (test-particle limit). While in the second implementation the back-reaction on the gas is included (full two-fluid). The advantages and disadvantages of the methods are discussed below. These implementations can be found in src/Common/Dust.cpp

\subsubsection{Test-particle implementation}

The test-particle implementation closely follows Booth et al (2015). Since in general the gas particles and dust particles will not be at the same locations, the gas properties needed for evaluating equation (7) must be interpolated to the location of the dust particle. This is done in the standard SPH way. For each dust particle we compute a `gas smoothing length' set by the local number density of gas particles, in a similar way to that of the gas. The gas properties are then interpolated to the location of the dust particles using a the standard kernel sum and the interpolated values are used in evaluating the drag acceleration.

\subsubsection{Full two-fluid algorithm}

The full two-fluid algorithm follows the improved semi-implicit algorithm of Lor\'en-Aguilar \& Bate (2015). In order to ensure conservation of both energy and angular momentum the forces must: (1) change sign under the exchange of particles, (2) be directed along the lines between interacting particles (Laibe \& Price 2012). (1) can be achieved simply by interpolation, but (2) requires forces to be directed between particles rather than in the direction of the velocity difference. Following Laibe \& Price (2012) and Lor\'en-Aguilar \& Bate (2015) we project the drag force between pairs of particles and evaluate the total acceleration using a kernel sum:
\begin{align*}
 \vec{a}_{d,\mathrm{drag}} &= - \sum_i^\mathrm{gas} \frac{m_i}{\rho_{g,i}} ( \vec{S}_{id}  \cdot \vec{\hat r}_{id}) \vec{\hat r}_{id} D(\vec{r}_{id}, h_i) \\
 \vec{a}_{g,\mathrm{drag}} &= - \sum_i^\mathrm{dust} \frac{m_i}{\rho_{d,i}} ( \vec{S}_{ig} \cdot \vec{\hat r}_{ig}) \vec{\hat r}_{ig} D(\vec{r}_{ig}, h_g),
\end{align*}
where $\vec{\hat r}_{ij} = (\vec{r}_i - \vec{r}_j) / |\vec{r}_i - \vec{r}_j|$ and
\begin{equation}
 \vec{S}_{ij} = \frac{\rho_i}{\rho_i + \rho_j} \frac{\langle \Delta \vec{v}_{ij} \rangle}{t_s}.
\end{equation}
For the kernel $D$ we use the double hump version of the kernel used for the SPH forces $D(r,h) \propto (r/h)^2 W(r,h)$, which reduces the bias in the direction of the force estimate from nearby neighbours when compared with bell-shaped kernels (Laibe \& Price, 2012).

\subsubsection{Advantages, limitations and reasons for caution}

\begin{itemize}
\item \emph{Conservation} -- The full fluid has the advantage that the total energy, momentum and angular momentum
are conserved explicitly, while in the test-particle limit these are conserved for the gas particles, but not for the dust. In some problems, such as when radial drift is important the explicit conservation can help maintain the correct solution.

\item \emph{Dissipation} -- The full two-fluid system is a dissipative system, with the drag forces converting kinetic energy into thermal energy. If the smoothing length, $h < c_s t_s$, this can lead to massive numerical dissipation especially when $\rho_d \approx \rho_g$. The semi-implicit algorithms help to reduce the dissipation (see Lor\'en-Aguilar \& Bate, 2014), but this can still be an issue. For $\rho_d \ll \rho_g$ the dissipation is weak, with no dissipation occurring in the test particle limit.

\item \emph{Force Errors} -- The standard SPH interpolation used in the test-particle implementation is generally accurate to $\sim 0.1$ per cent for the cubic spline, and even more accurate for smoother kernels (e.g. the Quintic spline). However, the projected interpolation is considerably more sensitive to the particle distribution. The double hump kernels help reduce the errors due to the particle distribution, but even for regular particle distributions the error can be $\sim 1$ per cent. When strong density gradients are present (e.g. the vertical structure of discs), the errors can be much larger (c.f the settling test, Lor\'en-Aguilar \& Bate, 2014). For this reason smoother kernels are highly recommended when using the full-two fluid implementation.

\item \emph{Individual particle time-steps} -- Individual particle time-steps are implemented for both dust algorithms, but typically the full-two fluid algorithm should not be used with individual particle time-steps, since the time averaged drag force is only correct once the time-step has been completed. Order unity errors can arise even the simplest dustybox tests when the time-step for the gas and dust differ for this reason, making it clear the exact conservation is necessary for accurate simulations with this algorithm. Development of a robust multi-step scheme for the full two-fluid scheme is being considered. Conversely, individual particle time-steps have been used in Booth, Sijacki \& Clarke (2015) and Booth \& Clarke (2016) with the test-particle limit algorithm and prove to be robust because the predicted velocities are not needed for force calculations. Individual particle time-steps have not been tested extensively in the GANDALF test-particle implementation, but preliminary tests have been successful.

\item \emph{Sinks} -- The sink routines have not been modified to support dust, but in principal should work. However, I would recommend checking this when both dust and sink particles are needed in the same simulation. Accretion using  the `vacuum cleaner' type sinks \emph{should} work properly, but has not been rigourously tested. The smoothed accretion mode is more likely to give issues since it removes mass from the particle closest to the sink only so it is possible that mass would end up being removed from the gas or dust particles only.

\end{itemize}

\subsection{Drag Laws}
Several drag laws have been implemented in via functions form of $t_s \equiv t_s (\rho_g , \rho_d , c_s)$. The different functional forms currently available can be chosen via the parameter options and are described below and are defined in src/Headers/DragLaws.h. Most of these drag laws are phenomelogical in form, and only work with dimensionionless form. The one exception is the Epstein drag law, in which the stopping time is given by:
\begin{equation}
t_s = \frac{3}{4} \sqrt{\frac{\pi}{8}} \frac{m}{A(\rho_g + \rho_d)c_s},
\end{equation}
where $m$ is the grain mass and $A$ is its grain area ($\pi a^2$ for spherical grains). In the code, absorb the area-to-mass ratio into a single coefficient.


\subsection{IO and initial conditions}
Dust particles are currently supported only in the SEREN formats (both formatted and unformatted), with the dust particles following the gas particles on the disc. The number of dust particles is stored in the header, as described in section \ref{Sec:SerenForm}. The column format is \emph{not} formally supported, but if used it will currently still contain the dust particles after the gas particles as long as MPI is not used. However, the number of dust particles is not currently saved in the column format.

Currently the only systematic way to generate initial conditions that include dust is to provide an IC file in the SEREN format that includes dust. For the test problems currently included, when dust is turned  on (via setting the drag\_forces parameter) the initial conditions generated are identical to those for the gas, although dust-to-gas ratio can be varied. Currently, generating multi-species initial conditions in python is also not supported (with the exception of star particles), although we hope to support this soon.

\subsection{Test Problems}
Some simple two-fluid test problems have been included in the directory tests/dust tests. Most of these are straight-forward modifications of classic hydrodynamics test problems to include dust particles. A description of a number of these can be found in Laibe \& Price (2012 a,b), including the \textsc{dustybox}, \textsc{dustywave}, \textsc{dustyshock} and \textsc{dustysedov} problems. The dust to gas ratio can be verified, but otherwise the dust initial conditions are identical to the those in the gas in all the problems provided.

\subsection{Dust parameter options}
\begin{itemize}

\item \var{drag\_forces} : Select the drag force algorithm

\begin{tabular}{ll}
none              & = No drag forces \\
test particle     & = Test particle limit drag regime \\
full twofluid     & = Full two fluid drag regime
\end{tabular}

\item \var{drag\_law} : Select the drag law

\begin{tabular}{ll}
fixed             & $t_s = 1 / K_D$ \\
density           & $t_s = 1 / (K_D(\rho_d + \rho_g))$ \\
epstein           & $t_s = \frac{(3/4) \sqrt{\pi/8}}{K_D(\rho_d + \rho_g)c_s}$, where $K_D$ is the area-to-mass ratio \\
LB2012            & $t_s = \rho_d \rho_g / (K_D(\rho_d + \rho_g))$ \\
\end{tabular}

\item \var{drag\_coeff} : Sets $K_D$

\item \var{dust\_mass\_factor} : Sets the dust-to-gas mass ratio for the test problems

\end{itemize}

Including additional drag forces is straight forward. First add a new class to src/Headers/DragLaws.h. The new drag law can then be added to DustFactory::ProcessParameters (at the bottom of src/Common/Dust.cpp), which is responsible for selecting the desired drag law.

\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Developer notes}
GANDALF makes heavy use of the object-oriented features in C++ in order to write a flexible code that can used in many different ways while maintaining good performance.  For developers that wish to modify sections of GANDALF, or add new classes or functions, it is recommended to read this section in order to understand some of the ...


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Class hieararchy}
GANDALF uses a hierarchy of classes using inheritance to structure the code in a logical and maintainable manner.  Although classes permeate the code at almost all levels, we briefly discuss the main classes here.

\subsubsection{Simulation class}
The \var{Simulation} class is the class that contains every other class, including SPH and N-body classes, in order to fully run the simulation.  The \var{Simulation} class is split into three levels :
\begin{itemize}
\item \var{SimulationBase} - Non-templated base class for generating the simulation class.  Provides binding layer between C++ code and Python (since Python cannot link to templated classes).  Contains all subroutines that can be called by Python.
\item \var{Simulation} - Main parent simulation class that contains all common functions for all possible simulations.  Inherits from \var{SimulationBase}.
\item \var{AlgorithmSimulation} - Child class containing specific algorithmic implementation.  For example, `grad-h' SPH simulations are run by the \var{GradhSphSimulation} class; pure N-body simulations are run by the \var{NbodySimulation} class, etc..
\end{itemize}

\subsubsection{Hydrodynamics class}

\subsubsection{SmoothingKernel class}

\subsubsection{NeighbourSearch class}

\subsubsection{EOS class}

\subsubsection{EnergyEquation class}

\subsubsection{Nbody class}

\subsubsection{Sinks class}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Templates}
Templates are used in C++ to write generic code that can be used by different data types throughout the code without needing to write multiple versions of essentially the same code.  It can also be used as an optimisation technique for performance, which we explain later.

The simplest case of using templates in GANDALF is for dimensionality.  For almost all of the code, we use the template parameter `\var{ndim}' to generate a version of the code for each dimensionality used, i.e. $1$, $2$ and $3$.  In principle, \var{ndim} could simply be a variable which is defined by the user depending on which dimensionality is used for that simulation.  However, the value of \var{ndim} would then need to be accessed in memory very frequently.  By defining \var{ndim} as a template parameter, \var{ndim} will be defined statically and therefore does not require the overhead of accessing the memory address each time.

Throughout the code, the variables that are used as template parameters are :
\begin{itemize}
\item \var{ndim} : Number of dimensions
\item \var{Kernel} : Smoothing kernel function used
\item \var{ParticleType} : Particle data structure for employed hydrodynamics method
\item \var{CellType} : Radiation transport tree cell-type, depending on employed RT method
\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Particle data structures}

To optimise the speed of the code, GANDALF uses data structures that contain only the quantities for that particular hydrodynamics algorithm.  For example, `grad-h' SPH in 1D requires a very different number of variables to Godunov SPH in 3D.  Therefore, a different version of each sub-routine that uses the SPH particle data directly is required using templates.

The various particle data structures are defined in the header file \var{`SphParticle.h'}.  All particles are derived from a base data-type called \var{`Particle}

%Need a figure showing the particle data structure inheritance hierarchy here in the future.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Particle array pointers}

If a class which does not know the exact particle data structure wishes to access an array, or a single array element.


\newpage


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Units and scaling} \label{S:UNITS}
In GANDALF, all physical calculations are done in dimensionless units for precision and accuracy reasons.  However, the scaling factors for converting the initial conditions and parameters (in physical units) are calculated automatically, once the user has specified the units of the quantities in the parameters file.


\subsection{Calculating scaling factors}
Most physical quantities in GANDALF are some combination of the three basic physical unit types, length, mass and time.  In principle, we are free to choose any scaling factor in order to convert to dimensionless units
\begin{equation}
r' = \frac{r}{R_0} \;\;\;\;\;\;\;\;
m' = \frac{m}{M_0} \;\;\;\;\;\;\;\;
t' = \frac{t}{T_0}
\end{equation}
where $r, m, t$ are the physical quantities, $r', m', t'$ are the dimensionless analogues and $R_0, M_0, T_0$ are the scaling factors.  The internal scaling factors for most quantities is some combination of the scaling factors of these three basic quantities.  For example, the scaling factor for the dimensionless velocity, $v' = v/V_0$ is equal to $V_0 = R_0/T_0$.


\subsection{Scaling factors for \var{G=1}}
If we take some equation using \var{G} (e.g. Newton's law of gravity) and then substitute in the scaling factors and rearrange, we obtain :
\begin{eqnarray}
a = \frac{G\,m}{r^2} \;\;\;\; \Rightarrow \;\;\;\;
\frac{R_0}{T_0^2}\,a' = \frac{G\,M_0\,m'}{{r'}^2\,R_0^2} \;\;\;\; \Rightarrow \;\;\;\; a' = \underbrace{\left\{\frac{G\,M_0\,T_0^2}{R_0^3} \right\}}_{G'}\,\frac{m'}{{r'}^2}\,.
\end{eqnarray}
The final dimensionless equation has the same form as the original equation, with all constants and scaling factors grouped together in the `dimensionless gravitational constant', \var{G'}.  In N-body codes, it is common to adjust one of the scaling factors in order to set this new constant equal to unity.  By convention, it is the time scaling factor that is adjusted :
\begin{eqnarray}
\frac{G\,M_0\,T_0^2}{R_0^3} = 1 \;\;\;\; \Rightarrow \;\;\;\;T_0 = \left( \frac{R_0^3}{G\,M_0} \right)^{1/2}
\end{eqnarray}
For example, typical N-body units (at least in star cluser/formation simulations) of $R_0 = 1\,{\rm pc}$ and $M_0 = 1\,{\rm M}_{\odot}$ give a time scaling factor of $T_0 = 14.91\,{\rm Myr}$.  It is important to realise that this then has a knock-on effect on any other physical quantity that contains the dimensions of time.  For example, the unit of velocity becomes $V_0 = R_0/T_0 = 0.065\,{\rm km\,s^{-1}}$.


\subsection{Temperature scaling factor}
In thermodynamics, one additional unit that is often used (more often in the output rather than in internal quantities) is the thermodynamic temperature.  This quantity represents an additional unit that cannot be represented in any combination of length, mass and time.  The conversion between the three principle quantities is achieved via the Boltzmann constant.  Just as with setting $G = 1$, we can also combine and convert quantities involving temperature to an appropriate dimensionless unit, $T' = T/\theta_0$, by effectively setting the quantity $k_b/m_h$ equal to unity, i.e.
\begin{equation}
c^2 = \gamma \frac{k_b\,T}{\bar{m}} \;\;\;\; \Rightarrow \;\;\;\;
V_0^2\,{c'}^2 = \gamma\,\frac{k_b\,\theta_0\,T'}{m_h\,\bar{\mu}}
\;\;\;\; \Rightarrow \;\;\;\;
{c'}^2 = \gamma \; \underbrace{ \left\{ \frac{k_b\,\theta_0}{m_h\, V_0^2} \right\} }_{1} \; \frac{T'}{\bar{\mu}}
\end{equation}
\begin{equation}
\frac{k_b\,\theta_0}{m_h\, V_0^2} = 1 \;\;\;\; \Rightarrow \;\;\;\; \theta_0 = V_0^2\,\frac{m_h}{k_b}
\end{equation}




\subsection{Computing scaling variables in GANDALF}
For consistency, GANDALF computes all scaling variables internally in SI units.  However, this is not useful for astrophysical applications where the natural units to choose may be parsecs or megayears.  Therefore, the internal scaling variables are split into two parts for convenience in chosing an appropriate.  For a given unit \var{X}, the scaling factor is split into $X_0 = X_{\rm outscale} \, X_{\rm outSI}$, where
\begin{itemize}
\item $X_{\rm outscale}$ : scale factor to convert $X'$ to $X$ in the user-requested units;
\item $X_{\rm outSI}$ : the requested unit in SI units.
\end{itemize}
To convert from the code units to the user-requested units, $X_{\rm user} = X_{\rm outscale}\,X'$.  Alternatively, if you wish to convert directly to SI units, then you must also multiply by $X_{\rm outSI}$, i.e. $X_{\rm user} = X_{\rm outscale}\,X_{\rm outSI}\,X'$.  One alternative possibility is if the user wishes to convert from code units to c.g.s. units; in which case, they must use $X_{\rm outcgs}$ in place of $X_{\rm outSI}$, i.e. $X_{\rm user} = X_{\rm outscale}\,X_{\rm outcgs}\,X'$

\subsection{Converting initial conditions to code units}
Converting from physical units to code units (e.g. when setting up initial conditions) is trivial once the scaling factors have been computed.  Assuming the initial conditions for all variables are in the same units as specified in the parameters file, then we must simply divide by the $X_{\rm outscale}$, i.e.
\begin{equation}
X' = \frac{X_{\rm user}}{X_{\rm outscale}}\,.
\end{equation}
If the variable (or parameter) is in either SI or cgs units, then they can be converted directly also using the SI or cgs scaling factors.  For example, to convert a cgs quantity to code units,
\begin{equation}
X' = \frac{X_{\rm user}}{X_{\rm outscale}\,X_{\rm outcgs}}\,.
\end{equation}

Inside the code, each unit is its own separate class defined in the src/Headers/SimUnits.h and src/Common/SimUnits.cpp files.  All units are then stored in a single container class, called `SimUnits'.  This is then passed around the simulation as an object with the name \var{simunits}.  The object name for each individual unit itself is usually a single character (e.g. \var{r} for length, \var{m} for mass, \var{t} for time, etc..  See the src/Headers/SimUnits.h files for a full list).  For example, the following line could be used to convert the particle mass from physical (user) units to code units,  \\
\newline
\noindent \var{part.m /= simunits.m.outscale;} \\

\noindent For converting from say cgs units to code units, we could write \\
\newline
\noindent \var{part.m /= (simunits.m.outscale*simunits.m.outcgs);} \\





\subsection{Input and output units}
GANDALF does support the ability to read in initial conditions in a different set of units to use specified in the parameters file.  At the moment, this is only implemented in the \var{sf/seren\_formatted} and \var{su/seren\_unformatted}.  In practice, this is implemented by noting that the input and output scaling factors should be exacty the same, i.e.
\begin{equation}
X_0 = X_{\rm outscale} \, X_{\rm outSI} = X_{\rm inscale} \, X_{\rm inSI}
\end{equation}




\newpage


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{To-do list}

\subsection{Known bugs}
List of known bugs as of version \VERNO.

\begin{itemize}
\item Saitoh \& Makino (2012) SPH has not been fully updated since recent code refactoring.
\item For very short simulations, plots may not be updated correctly since the simulation process finishes before the plotting process has recieved its commands.
\item If running a simulation in interactive mode and a different simulation is loaded into memory, then it is no longer possible to continue running that simulation.
\item Rendered images are technically not done correct if smoothing lengths are smaller than the grid size (which can often be the case, although the images are fine for viewing/movie purposes).
\item Rendered images do not correctly calculate automatic limits; must be inserted manually using the limit command.
\end{itemize}


\subsection{Proposed features}
List of possible new features for future versions.

\begin{itemize}
\item Finish implementation of MPI
\item Allow ability to add titles to graphs in matplotlib
\item Include analytical solutions for Sedov blast-wave test and freefall collapse test
\item Plot star/sink-based statistics, e.g. sink mass-functions, binary statistics
\item Include analysis routines for python for N-body specific statistics, e.g. Q-parameter, $\lambda$-parameter, etc..
\item Strict energy-error checking for N-body simulations
\item More sanity-checking, error-trapping and assert statements (to prevent crashing on erroneous input and to help debugging purposes).
\end{itemize}

\newpage

\appendix

%\part*{Appendix}
%\phantomsection
%\addcontentsline{toc}{part}{Appendix}
%\appendix
%\def\maketitle{}
\def\tableofcontents{}
\section{Command reference for the python functions}
\label{reference_facade}
\input{../analysis/_build/latex/GANDALF}
\bibliographystyle{mnras}
\bibliography{biblio}

\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
