Introduction
SMILE is a software for orbit analysis and Schwarzschild modelling of stellar systems. The primary applications are:
- Exploring orbits in various potentials, either given by analytical formulae or several approximate expansions;
- Studying properties of orbits, computed by the builtin orbit integrator or obtained from external data;
- Constructing self-consistent equilibrium models of triaxial stellar systems with a given density profile using the Schwarzschild method.
The software comes in two versions – the GUI interactive program (implemented in CSmileGUI class), which forms an interface over the main business logic layer implemented in CSmileCore class, and a console tool whose functionality is implemented in CSmileConsole class. In addition, there are several standalone utility programs:
- mkspherical.cpp (create spherical mass models from a text table or from N-body snapshots)
- renderdensity.cpp (create a visual representation of a non-spherical density model specified by a BSE/Spline potential)
- measureshape.cpp (measure the axis ratio of an N-body snapshot as a function of radius)
- testaccuracy.cpp (check the accuracy of potential approximation constructed from an N-body snapshot)
- snaporbits.cpp (perform orbit analysis on trajectories of particles from an N-body snapshot)
- energydiff.cpp (measure the rate of energy diffusion in N-body simulation)
The user's guide and the mathematical background behind several algorithms employed in the software can be found here, while the scientific description of Schwarzschild modelling is given in the code paper.
Obtaining and compiling the software
The software can be obtained at http://td.lpi.ru/~eugvas/smile/. Compiled versions for several platforms are available; if you need to compile it from source then read this section.
The following third-party libraries are required:
The following libraries or external programs are optional:
Check smile.pri for the correct paths to libraries and set the compilation options HAVE_*** according to the available optional components, then run qmake smile.pro
, then make
; repeat for qmake smilec.pro
.
Version history
v1.0 (March 2011): initial release (not advertised anywhere)
v1.1 (July 2011): implemented a broader class of basis-set expansion, fixed some bugs.
v2.0 (1 August 2013): first major public release, accompanies the paper [Vasiliev, 2013, MNRAS]
- program structure substantially reworked towards modularity; it should be possible to use only potential or orbit integrator/analysis in other programs by including a subset of source files.
- New class of general-purpose potential approximation (Spline spherical-harmonic expansion).
- Choice of symmetry degree for BSE/Spline (from 'None' through standard 'Triaxial' to 'Spherical').
- Schwarzschild modelling: separated solver (several possible variants) from data objects and the model itself;
- different choices for Schwarzschild model type (classical grid-based, basis set and spherical harmonic expansion).
- Different method of generating initial conditions for Schwarzschild model (based on distribution function for an equivalent spherical model, constructed by Eddington inversion formula in CMassModel class).
- A number of features were scrapped: initial conditions for SM are only generated randomly, no stationary/principal-plane start space anymore (it did not produce good models anyway), no 'orbit bunches' (or 'dithering', which averages over several nearby orbits: just increase the total number of orbits, and quadratic optimization will drive the orbit weights towards more uniformity). Lucy iterations optimization removed – was too slow on realistic problem sizes and never outperformed quadratic optimization.
v2.1 (late 2013): non-public, accompanies the paper [Vasiliev,Antonini,Merritt, 2014, ApJ]
- added an option for initial conditions refinement at the center.
- added Schwarzschild data object for constraining the relative angular momentum distribution at a given energy.
v2.2 (March 2014): non-public, accompanies the papers [Vasiliev, 2014, Clas.Quant.Grav. and Vasiliev, 2014, MNRAS]
- several new orbit integrators using
odeint
library from boost
.
- added velocity perturbations (for the non-spherical Monte-Carlo code
RAGA
).
- improved robustness of handling of spherical mass models.
v2.5 (1 February 2015): second major public release, accompanies the paper [Vasiliev&Athanassoula 2015] substantially expanded feature set:
- added figure rotation to the potential.
- added several new potential and density models: Ferrers, Miyamoto-Nagai, Sersic, ExpDisk; potential expansion for finite-extent models based on spherical harmonics and Bessel functions; cylindrical spline potential expansion for highly flattened models;
- new variants of Schwarzschild density data object for cylindrical geometry.
- reworked the interface for initial conditions generation: better sampling of density profile (based on self-tuning Monte-Carlo sampling), spherical and axisymmetric Jeans models for velocity sampling.
- support for multi-component potential and Schwarzschild models.
- new orbit integrators: IAS15 and Hermite.
- characteristic diagram of periodic orbits (presently not used for generating initial conditions), routine for locating periodic orbits in the Poincare surface of section.
- SMILEPOT library exposing the potential machinery for external applications.
Overview of the software structure
The class hierarchy of SMILE is designed to have reasonably abstract interfaces between different parts of the code.
- Classes derived from CDensity and CPotential are used to represent a density model or a density-potential pair; in the latter case functions are provided to compute potential and forces at a given point, which are necessary for orbit integration. These classes are 'standalone' and may be used in external software without carrying all the dependencies from the entire SMILE machinery.
- COrbit is an interface class between orbit integrators (derived from CBasicOdeIntegrator) and orbit analysis routines (derived from CBasicOrbitRuntimeFnc), the exact nature of their interaction is described below. It contains the initial data for the orbit (COrbitInitData) and the results of integration and classification.
- The orbit integrators are advancing the solution with their internal, automatically chosen timestep; after each timestep they call one or more orbit runtime functions that are attached to the orbit. These functions, in turn, may obtain position and velocity at any moment of time within the given timestep (for example, at regular intervals of time independent of the actual integration timestep).
- The runtime functions (CBasicOrbitRuntimeFnc) thus serve dual purpose: they collect data during the orbit integration, and then may perform some kind of analysis upon completion. The results of this analysis are passed around as instances of objects derived from CBasicInformation (generally, each runtime function class has a corresponding information class). The runtime functions attached to COrbit persist after finishing the integration, as long as the COrbit object exists. (This makes sense only in the GUI for studying a single orbit, but not in the context of creating the orbit library).
- The abstractly defined data in the array of CBasicInformation objects is attached to an orbit after the integration, and may be stored or loaded from disk by the I/O routines in iodata.h
- COrbitDesc is another intermediate class that represents the orbit initial data and the same array of information objects as COrbit itself. Internally, it creates an instance of COrbit class to do the integration and analysis; however, it deletes this object upon finishing the integration and thus not retain runtime functions attached to COrbit; only the information resulting from analysis is kept.
- COrbitLibrary is a collection of COrbitDesc objects, and functions to assign initial conditions to them. Before the integration of the library, only the initial conditions are assigned; after completion each COrbitDesc keeps the analysis information.
- Another intermediate class family is derived from CBasicOrbitRuntimeFncCreator. Each class of COrbitRuntimeFnc has its own creator class, whose sole purpose is to create an instance of runtime function for each orbit. An array of creator classes is passed to the COrbitLibrary and down to COrbitDesc and COrbit, so each orbit creates the array of runtime functions for itself. (Such an intermediate step is required because runtime functions know only about the particular orbit, not the entire orbit library).
- The classes for performing the Schwarzschild modelling are also split into several groups. There are classes for dealing with various Schwarzschild Data, classes for performing the modelling, and interface classes for various third-party optimization solvers. Let's consider them in turn. We define the process of Schwarzschild modelling as matching the source data (represented by the orbit library and the associated orbit analysis information) to the target model, which consists of one or more data objects describing the desired properties of the galaxy.
- Abstract CBasicSchwData class defines a general way of providing the required properties of the model, e.g. the density distribution parametrized as the masses of cells of spatial grid in the classic Schwarzschild model, or some sort of kinematic data coming from the observations (although this is not implemented in the current version). The data objects may take the potential/density class as the input parameter and compute the required quantities, which we refer to as constraints. In addition, they provide methods to compute the same quantities for each orbit. This is mediated by a corresponding orbit runtime function; for each kind of Schwarzschild data class there is a runtime function that computes the necessary data for an orbit and stores it in the CSchwInformation class in the form of a one-dimensional array.
- The modelling classes, derived from CBasicSchwModel, take the array of CBasicSchwData objects as the target, and the COrbitLibrary, in which all orbits have necessary kind of CSchwInformation for every CBasicSchwData object in the model, as the source. The modelling itself amounts to fitting the source to the target, i.e. finding the orbit weights that minimize or nullify deviations of constraints combined from all orbits from their required values given in schw.data objects. The precise method of finding the solution, and additional conditions imposed on the distribution of orbit weights, are left for specification in the derived classes. In the present version, there is only one class, CSchwModelQuadOpt, which uses the linear or quadratic optimization as the way of finding the solution. The task of the modelling object is to combine all data objects and orbits into a single matrix equation to be passed to the solver. Meanwhile, to save memory by not copying large amount of data in matrices, this matrix equation is represented as an abstract matrix interface (derived from CMatrix) rather than an actual two-dimensional array of numbers. Additional preferences about what orbits should have more weight in the model are provided by a CBasicOrbitFilteringFnc (see below).
- The quadratic optimization solver used in the modelling is represented by an instance of a class derived from CBasicQuadraticOptimizationSolver. This class has a single method to call the solver for a given matrix equation, which is taken from the modelling object via CMatrix interface. The derived classes present a unified interface to a particular third-party linear or quadratic solver used as the back-end number cruncher.
- Of the classes that we have not considered so far, there are several simple ones that represent single point masses and point mass arrays (CPointMassSet). They are used to read and save Nbody snapshots. Reading snapshots is used to construct CPotential objects of one of the general-purpose potential solvers from a mass model specified by a discrete point mass set. Writing snapshots is used in converting the Schwarzschild model to an Nbody model, which is performed by COrbitLibrary.
- Reading and writing Nbody snapshots is done through classes derived from CBasicIOSnapshot, which handle the representation of snapshots in various formats (text, NEMO, GADGET).
- Classes derived from CBasicOrbitFilteringFnc provide a general interface for evaluating an orbit according to some criterion and obtain a numerical value for it. One example is the weighting of orbits in the Schwarzschild model; another is used in the context of generation of the Nbody model with unequal masses of particles (mass refinement). The refinement strategy is provided by a CMassRefinementFnc.
- Input/output of the orbit library and all data associated with orbits is provided by a pair of classes, CFileInputText / CFileOutputText.
- Finally, to put everything together, there is CSmileCore class that provides all basic functionality for both GUI and console versions of the program. This is (almost) the only place where Qt framework is used, for creating helper threads to perform various tasks in parallel with the main program, passing messages between threads, handling the configuration information stored in INI files, etc.
- The access to the CSmileCore object is provided by CSmileGUI or CSmileConsole, depending on the variant of the main program.
Ufff.. that's the basic outline of the relationships between various classes in SMILE; for details refer to the documentation of particular modules and classes. It looks quite complex, but probably is inevitable to achieve the necessary level of generality and abstraction.
Coding conventions
- The program makes extensive use of C++ concepts such as classes, templates, STL algorithms, etc. However, some parts that interface with C code are inevitable (in particular, GSL routines often require wrappers to deal with data structures belonging to object instances; also the old C code is used in the DOP853 orbit integrator). All memory allocation is done in terms of
new
/delete
operators or by std::vector
s; we don't use C-style malloc
/free
or char* strings, nor the more advanced "smart pointers" (or anything else from boost
, except for the optional odeint
library).
- Exceptions are not used consistently throughout the software. Most often, routines that may generate an out-of-memory or other notable exception, will wrap the internal code in a try-catch block and will safely return an error code instead of passing out the exception.
- There are several class hierarchies derived from abstract classes (named CBasic***). Often the exact nature of the instance of a derived class needs to be known by a non-class function; in these cases, the class type is provided by some sort of getType() member functions; these return a enum value from a pre-defined list, which is declared in the base class. Thus we do not rely on RTTI to obtain the typename of a class, but rather use this contrived mechanism that needs to be manually extended in the case that a new derived class is implemented.
- const correctness is maintained throughout the code. Member functions that do not change the internal state of the object have
const
qualifiers. Input parameters are labelled as const myType& myVar
or const myType* myVarPtr
in the case that myType
is not a simple type. Output parameters are passed as myType* outVar
regardless of whether myType
is a simple or complex type. (Of course, if there is only one output variable, it is returned as the return type).
- On several occasions, numerical data may be provided both in
float
or double
variants. For these cases, we use template classes with template parameter NumT; often they have templated constructors to allow a seamless conversion from a different data type.