Adding a new geometry class [TUTORIAL]

Illustration above: exponential disk; density cut along $xz$-plane (left) and edge-on surface brightness (right).

Illustration above: isothermal disk; density cut along $xz$-plane (left) and edge-on surface brightness (right).

In this tutorial, you will extend the SKIRT code with a new geometry that can be used for the spatial distribution of sources or media. Specifically, you will implement the isothermal disk as illustrated above.

Getting ready

Before starting this tutorial, you should know how to use SKIRT as described in the User Guide and its various tutorials. You should be ready to use a text editor or an integrated development environment to modify portions of the SKIRT source code and rebuild the project as described in the Installation Guide.

Background information

As a first step, please review The ins and outs of simulation items, which provides the background information required to follow along with this tutorial.

In summary, a SKIRT simulation is represented in memory as run-time hierarchy of objects, called simulation items, reflecting the structure of a ski file. Each simulation item is an instance of a C++ class that inherits the SimulationItem class. Thus, simulation items form a compile-time class hierarchy through inheritance (with the SimulationItem class at the top), and a run-time object hierarchy through aggregation (with an instance of a Simulation subclass at the top).

The SKIRT code includes a mechanism that allows all information about a simulation item class to be defined in a single place. All relevant metadata, such as property definitions and human readable strings to describe these properties, are provided in the class definition itself rather than being scattered across the code.

The mathematics

The new geometry is an isothermal disk with the axisymmetric density profile

$$\rho (\mathbf{r})\equiv \rho (R,z)={\rho }_0\exp (-\frac{R}{h_R}){\text{sech}}^2\left(\frac{z}{h_z}\right)$$

with the scalelength $h_R$ and the scaleheight $h_z$ as two free parameters.

Density

The Geometry class documentation states that all geometries in SKIRT must be normalized so that the total mass is equal to one. Given this normalization, we can determine the central density for the new geometry as a function of $h_R$ and $h_z$. We find

$$M=\iiint \rho (\mathbf{r})\text{d}\mathbf{r}=2\pi {\rho }_0{\int }_0^{\mathrm{\infty }}\exp (-\frac{R}{h_R})R\text{d}R{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\text{sech}}^2\left(\frac{z}{h_z}\right)\text{d}z.$$

Using the equalities

$${\int }_0^{\mathrm{\infty }}x{\text{e}}^{-x}\text{d}x={\int }_0^{\mathrm{\infty }}{\text{sech}}^{2}x\text{d}x=1$$

we find that

$$M=4\pi h_R^2{h}_z{\rho }_0$$

and hence (using $M=1$)

$${\rho }_0=\frac{1}{4\pi h_R^2{h}_z}.$$

The density $\rho (R,z)$ is then easily calculated using its original definition above.

Surface densities

We also need to calculate the edge-on and face-on surface densities ${\mathrm{\Sigma }}_R$ and ${\mathrm{\Sigma }}_Z$. The documentation for the AxGeometry::SigmaR() function defines the edge-on or radial surface density as the integration of the density along a line in the equatorial plane starting at the centre of the coordinate system,

$${\mathrm{\Sigma }}_R={\int }_0^{\mathrm{\infty }}\rho (R,0)\text{d}R.$$

The documentation for the Geometry::SigmaZ() function defines the face-on or Z-axis surface density as the integration of the density along the entire Z-axis,

$${\mathrm{\Sigma }}_Z={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\rho (0,z)\text{d}z.$$

In our present case, we find

$${\mathrm{\Sigma }}_R={\rho }_0{\int }_0^{\mathrm{\infty }}\exp (-\frac{R}{h_R})\text{d}R={\rho }_0{h}_R$$

and

$${\mathrm{\Sigma }}_Z={\rho }_0{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\text{sech}}^2\left(\frac{z}{h_z}\right)\text{d}z=2{\rho }_0{h}_z.$$

Random position

Finally a geometry must be able to draw a random position from its three-dimensional probability density distribution. For an axisymmetric geometry with separable density function, the problem is reduced to drawing two values, each from a one-dimensional probability distribution, in the radial and vertical directions.

The documentation for the SepAxGeometry::randomCylRadius() function states that it returns a random $R$ from the one-dimensional probability distribution

$$p(R)\text{d}R=2\pi {\rho }_R(R)R\text{d}R.$$

In the present case we find

$$p(R)\text{d}R=\frac{\rho (R,z)R\text{d}R}{{\int }_0^{\mathrm{\infty }}\rho (R^{\prime },z)R^{\prime }\text{d}{R}^{\prime }}=\frac{1}{h_R^2}\exp (-\frac{R}{h_R})R\text{d}R$$

with cumulative distribution

$$P(R)={\int }_0^{R}p(R^{\prime })\text{d}{R}^{\prime }=1-(1+\frac{R}{h_R})\exp (-\frac{R}{h_R}).$$

Solving $\mathcal{X}=P(R)$ for $R$, where $\mathcal{X}$ is a uniform random deviate, we obtain

$$R=h_R[-1-W_{-1}\left(\frac{\mathcal{X}-1}{\text{e}}\right)]$$

where $W_{-1}$ is the Lambert function of order $-1$, which can be calculated in SKIRT using the function SpecialFunctions::LambertW1().

The documentation for the SepAxGeometry::randomZ() function states that it returns a random $z$ from the one-dimensional probability distribution

$$p(z)\text{d}z={\rho }_z(z)\text{d}z.$$

In the present case we find

$$p(z)\text{d}z=\frac{\rho (R,z)\text{d}z}{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\rho (R,z^{\prime })\text{d}{z}^{\prime }}=\frac{1}{2h_z}{\text{sech}}^2\left(\frac{z}{h_z}\right)\text{d}z$$

with cumulative distribution

$$P(z)={\int }_{-\mathrm{\infty }}^{z}p(z^{\prime })\text{d}{z}^{\prime }=\frac{1}{2}[1+\text{tanh}\left(\frac{z}{h_z}\right)].$$

Solving $\mathcal{X}=P(z)$ for $z$, where $\mathcal{X}$ is a uniform random deviate, we obtain

$$z=h_z\text{artanh}(2\mathcal{X}-1)=\frac{1}{2}{h}_z\ln \frac{\mathcal{X}}{1-\mathcal{X}},$$

where we used the identity

$$\text{artanh}x=\frac{1}{2}\ln \frac{1+x}{1-x}.$$

Scoping out the work

To add a new geometry, you need to create a new class and tie it into the simulation item class hierarchy. Rather than reinventing the wheel, browse the SKIRT documentation and look for an existing class with properties similar to your needs. It so happens that the ExpDiskGeometry class also defines an axisymmetric disk. To examine it more closely, open the ExpDiskGeometry.hpp header file in the code editor.

It turns out that ExpDiskGeometry inherits from SepAxGeometry, an abstract class representing axisymmetric geometries where the density is a separable function of $R$ and $z$. The PROPERTY and ATTRIBUTE macros define properties for the scale length $h_R$ and the scale height $h_z$, corresponding to the free parameters you need, plus several properties related to inner and outer truncation of the disk, which you don't need. These macros actually provide the metadata used by the SMILE library, in addition to defining appropriate public getter functions with the same name as the corresponding property.

Furthermore, ExpDiskGeometry offers the following functions:

• the setupSelfBefore() function, which calculates the central density ${\rho }_0$ and stores the result in a data member for later use;
• the functions density(), randomCylRadius(), randomZ(), SigmaR(), and SigmaZ(), which implement various aspects of the geometry modeled by this class.

You can use the same structure for your new class. The main work will consist of adjusting the key functions for the new geometry.

Creating the new source and header files

Locate the files ExpDiskGeometry.hpp and ExpDiskGeometry.cpp in your local working copy of the SKIRT source code (in the SKIRT/core directory), and copy them to a set of new files named IsoThermalDiskGeometry.hpp and IsoThermalDiskGeometry.cpp in the same directory. Do not remove the original files!

In the new header file (.hpp), replace the two occurrences of "EXPDISKGEOMETRY_HPP" by "ISOTHERMALDISKGEOMETRY_HPP". This is the standard header guard to avoid multiple includes of the same class definition.

In both new files (.hpp and .cpp), replace all occurrences of "ExpDiskGeometry" by "IsoThermalDiskGeometry".

Build the SKIRT project with these changes. Depending on the build procedure in use, you may need to explicitly rerun CMake for the build system to recognize the newly introduced source files. To make sure that the new files are actually being compiled, introduce a syntax error and verify that the compiler catches it. After correcting these intentionally introduced issues, there should be no compilation or build errors left at this point.

The class header

Because the class definition for IsoThermalDiskGeometry is very similar to that of ExpDiskGeometry, the header file IsoThermalDiskGeometry.hpp for your new class does not need a lot of work.

Adjust the description of the class in the ITEM_CONCRETE macro. For example, replace "an exponential disk geometry" by "an iso-thermal disk geometry". This text will be used for identifying your new geometry to a user, such as when constructing a new .ski file during an interactive Q&A session.

Remove the PROPERTY and ATTRIBUTE macros for the truncation parameters, and the corresponding data members, since these parameters are not needed in the new class.

While updating the documentation, take care to follow the rules for comments in "/ **"-style documentation blocks. These comments are processed by the Doxygen utility to produce the SKIRT manual in HTML format. It is really important to write decent comments for every function. When in doubt about the syntax for formatting comments, look for examples elsewhere in the SKIRT source code or check the Doxygen manual on the web.

The class implementation

The source file IsoThermalDiskGeometry.cpp needs more extensive edits beyond replacing the class name.

The implementation of the function setupSelfBefore() must always start with a call to the same function in the base class (SepAxGeometry in this case). This call should remain unchanged. Remove the code that verifies the values of the truncation-related properties, and update the code that calculates the central density:

void IsoThermalDiskGeometry::setupSelfBefore()
{
    SepAxGeometry::setupSelfBefore();
 
    // calculate central density
    _rho0 = 1.0/(4.0*M_PI*_hz*_hR*_hR);
}

Update the implementation of the functions density(), sigmaR(), and sigmaZ() to reflect the (surface) density formulas derived before:

double IsoThermalDiskGeometry::density(double R, double z) const
{
    double sechz = 1.0 / cosh(z/_hz);
    return _rho0 * exp(-R/_hR) * sechz*sechz;
}
 
double IsoThermalDiskGeometry::sigmaR() const
{
    return _rho0*_hR;
}
 
double IsoThermalDiskGeometry::sigmaZ() const
{
    return 2.0*_rho0*_hz;
}

Update the implementation of the functions randomCylRadius() and randomZ() to reflect the formulas derived before for generating the appropriate random values:

double IsoThermalDiskGeometry::randomCylRadius() const
{
    double X = random()->uniform();
    return _hR * (-1.0-SpecialFunctions::LambertW1((X-1.0)/M_E));
}
 
double IsoThermalDiskGeometry::randomZ() const
{
    double X = random()->uniform();
    return 0.5 * _hz * log(X/(1.0-X));
}

Save your changes and perform a build. There should be no errors.

Registering the new class

Having declared and defined the new class, you must make its existence known to the SMILE library so that it will be properly listed when constructing a new .ski file during an interactive Q&A session. Note that this registration requirement applies only to simulation item classes, i.e. classes that inherit directly or indirectly from the SimulationItem class, as is the case for your new IsoThermalDiskGeometry class.

Locate and open the file SimulationItemRegistry.cpp (also in the SKIRT/core directory). You need to add two lines in this file: an #include directive for the new class header, and an add<> statement to register the new class.

Add the #include directive to the already long list of #include directives, in alpahabetical order. For example:

//...
#include "IsoThermalDiskGeometry.hpp"
#include "IsotropicAngularDistribution.hpp"
//...

Insert a new add<> statement to the existing list, in the location you want it to appear in the corresponding choice list while creating a new .ski file in the interactive Q&A session. In this case, locate the list of add<> statements for geometry classes and insert a new statement (for example) right after the one for the ExpDiskGeometry class:

//...
ItemRegistry::add<ExpDiskGeometry>();
ItemRegistry::add<IsoThermalDiskGeometry>();
//...

Again save your changes and perform a build. Assuming there are no errors, SKIRT is now ready for use with the new geometry.

Using the new class

Launch the newly built SKIRT version without any command line arguments. Provide appropriate responses to the interactive questions for a simulation with a source or medium component using a built-in geometry. The list of choices will include a new item for the iso-thermal disk geometry:

   Possible choices for the geometry of the spatial luminosity distribution for the source:
      1. A Plummer geometry
      2. A gamma geometry
         ...
      8. An exponential disk geometry
      9. An iso-thermal disk geometry
         ...
 ? Enter one of these numbers [1,21] (1): 9

Enter the appropriate choice to select your new geometry and complete the Q&A session to save the ski file. Alternatively, you can edit a copy of a ski file you made earlier to use your new geometry.

To verify the results for a primary source distribution, run a simulation without media and with some appropriate FrameInstrument instances to generate images of the surface brightness observed from a few relevant viewing angles.

To verify the results for a medium density distribution, include an instance of the ConvergenceCutsProbe probe to generate cuts through the medium along the coordinate axes.

Congratulations, you made it to the end of this tutorial!