Adding New Classes

Code Directives

Galacticus is designed to be flexible and extensible, allowing you to add new classes and functionality without having to hack the code extensively. To achieve this it makes much use of embedded code directives which, for example, explain to the build system how a particular subroutine or function connects into the Galacticus code. Such code directives are enclosed in blocks delimited by !![ and !!], and take the form of short blocks of XML. For example, a typical code directive might look like:

!![
<accretionDisks>
  <unitName>Accretion_Disks_Shakura_Sunyaev_Initialize</unitName>
</accretionDisks>
!!]

This directive would typically appear just prior to a subroutine which initializes the Shakura-Sunyaev accretion disk module (it could appear anywhere throughout that module, but it makes sense to keep it close to the subroutine that it references). The accretionDisks tag explains to the Galacticus build system that this module contains an implementation of black hole accretion disks. The unitName tag specifies the name of a program unit which (in this case) should be called to initialize this accretion disk implementation. The build system will then insert appropriate use and call statements into the Galacticus code such that this routine will be called if and when accretion disks are required by Galacticus.

Identifying Components and Mass Types

Many functions can be applied to different components or groups of components and to different types of mass within a node. In general, these functions make use of a set of label defined in the Galactic_Structure_Options module. Components are identified by a componentType label which can take on the following values:

componentTypeAll

All components are matched;

componentTypeDisk

Only disk components are matched;

componentTypeSpheroid

Only spheroid components are matched.

componentTypeBlackHole

Only black hole components are matched.

componentTypeHotHalo

Only hot halo components are matched.

componentTypeDarkHalo

Only dark matter halo components are matched.

Types of mass are identified by a massType which can take one of the following values:

massTypeAll

All mass is included;

massTypeDark

Only dark matter is included;

massTypeBaryonic

Only baryonic mass is included;

massTypeGalactic

Only galactic mass is included.

massTypeGaseous

Only gaseous mass is included.

massTypeStellar

Only stellar mass is included.

massTypeBlackHole

Only black hole mass is included.

Components

This section describes the internal structure of node components, and how a component is implemented.

Component Structure

Each node in the merger tree consists of an arbitrary number of “components”, each of which can actually be an array, allowing multiple components of a given class. Each component represents a specific class of object, which could be a dark matter halo, a galactic disk or a black hole etc. A component of each class may be of one or more different implementations of that component class. Component classes are extensions of the nodeComponent base class, while each implementation is an extension of its component class (or, sometimes, of another implementation of that same class). Each component implementation type consists of a set of data[1], representing the properties (mass, size etc.) of the component, along with the rates of change (and ODE solver tolerances) for any properties which are evolvable. Additionally, each component contains a large number of methods (functions) which can be used to access its properties, query its interfaces and which are used internally to perform ODE evolution, output etc. The nodeComponent base class and all classes derived from it are built automatically by Galacticus::Build::Components at compile time (take a look in work/build/objects.nodes.components.Inc if you want to see the generated code).

Extending Components

It is possible to create a component which extends an existing component (see the discussion of the extends element in Section Component Definition). This capability is intended to allow new properties to be added to a component without having to create a whole new copy of the component. It is not intended to allow changes in the way in which the component is evolved through the halo hierarchy. (With the exception that rules to describe how the newly added properties will evolve through the halo hierarchy can be added of course.)

A simple example of this extension capability can be found in the scaleShape dark matter profile component, which extends the scale dark matter profile component. In this case, the scaleShape component adds a new property, shape, and specifies how it is to be initialized, evolved, output, and change by node promotion events. It does not affect how the scale property, inherited from the scale dark matter profile component, is evolved.

Implementing a New Component

Implementing a new component involves writing some modules and functions which contain a definition of the component and, if necessary, handle initialization, creation, evolution, and responses to any events. Frequently, the easiest way to make a new component is to copy a previously existing one and modify it as needed. Details of the various functions that component modules must perform are given below. By convention, a component implementation lives in its own directory, objects/nodes/components/<component>/<implementation>/ (with <component> and <implementation> acting as placeholders for the component class and implementation names, e.g. disk/standard), and is split into three or four files, although some components might not need all of these files. These files are named as follows:

objects/nodes/components/<component>/<implementation>/_class.F90

The primary file which describes the component and its properties, and which contains functions that manipulate the component as it evolves through a merger tree (ODE rates, behavior during mergers, etc.);

objects/nodes/components/<component>/<implementation>/bound_functions.Inc

Contains functions which will be bound to the component object (i.e. the nodeComponent<Class><Implementation> class), and so will be available as type bound procedures. Generally, these functions will include any which get or set values of properties in the component, those which return information about its internal state (such as a massDistributionClass object describing the structure of the component), and any other functions which we may want to be overridden by extensions to the component.

objects/nodes/components/<component>/<implementation>/data.F90

Contains any data which may need to be shared between the above two files. This might contain parameters which control some property of the component that is the same for all instances (e.g. if spheroids are modeled as Sérsic profiles all with the same value of the Sérsic index, that value might be placed into this file).

objects/nodes/components/<component>/<implementation>/structure.F90

Contains any functions which implement the structure (e.g. density, rotation curve) of the component and which cannot be placed in objects/nodes/components/<component>/<implementation>/bound_functions.Inc due to dependencies on modules which in turn depend on the Galacticus_Nodes module.

In general, the _class.F90 file is the place for the component definition and functions which process the component during tree evolution (including output), while bound_functions.Inc is intended for functions which record or report the internal state of the component.

Component Definition

Component definition itself takes the form of an embedded XML document. The following example illustrates such a document:

!![
<component>
 <class>disk</class>
 <name>exponential</name>
 <isDefault>yes</isDefault>
 <properties>
  <property>
    <name>isInitialized</name>
    <type>logical</type>
    <rank>0</rank>
    <attributes isSettable="true" isGettable="true" isEvolvable="false" />
  </property>
  <property>
    <name>massStellar</name>
    <type>real</type>
    <rank>0</rank>
    <attributes isSettable="true" isGettable="true" isEvolvable="true" isNonNegative="true" />
    <output unitsInSI="massSolar" unitsDescription="Solar masses" unitsQuantity="solMass" comment="Mass of stars in the exponential disk."/>
  </property>
  <property>
    <name>abundancesStellar</name>
    <type>abundances</type>
    <rank>0</rank>
    <attributes isSettable="true" isGettable="true" isEvolvable="true" isNonNegative="true" />
    <output unitsInSI="massSolar" unitsDescription="Solar masses" unitsQuantity="solMass" comment="Mass of metals in the stellar phase of the exponential disk."/>
  </property>
  <property>
    <name>massGas</name>
    <type>real</type>
    <rank>0</rank>
    <attributes isSettable="true" isGettable="true" isEvolvable="true" createIfNeeded="true" isNonNegative="true" />
    <output unitsInSI="massSolar" unitsDescription="Solar masses" unitsQuantity="solMass" comment="Mass of gas in the exponential disk."/>
  </property>
  <property>
    <name>coolingMass</name>
    <attributes isSettable="false" isGettable="false" isEvolvable="true" isDeferred="rate" isVirtual="true" />
    <type>real</type>
    <rank>0</rank>
  </property>
  <property>
    <name>halfMassRadius</name>
    <attributes isSettable="false" isGettable="true" isEvolvable="false" isVirtual="true" />
    <type>real</type>
    <rank>0</rank>
    <getFunction>Node_Component_Disk_Exponential_Half_Mass_Radius</getFunction>
  </property>
  <property>
    <name>luminositiesStellar</name>
    <type>real</type>
    <rank>1</rank>
    <attributes isSettable="true" isGettable="true" isEvolvable="true" isNonNegative="true" />
    <classDefault modules="Stellar_Population_Properties_Luminosities" count="Stellar_Population_Luminosities_Count()">0.0d0</classDefault>
    <output labels="':'//Stellar_Population_Luminosities_Name({i})" count="Stellar_Population_Luminosities_Count()" condition="Stellar_Population_Luminosities_Output({i},time)" modules="Stellar_Population_Properties_Luminosities" unitsInSI="luminosityZeroPointAB" unitsDescription="AB-magnitude zero point" unitsQuantity="4.465920e17 W/Hz" comment="Luminosity of disk stars."/>
  </property>
 </properties>
 <bindings>
  <binding method="attachPipes" function="Node_Component_Disk_Exponential_Attach_Pipes" type="void" />
 </bindings>
 <functions>objects/nodes/components/disk/exponential/custom_methods.inc</functions>
</component>
!!]

The elements of this document have the following meaning:

class

[Required] Specifies the component class of which this is an implementation.

name

[Required] Specifies the name of this specific implementation.

extends

[Optional] If present, this element must contain class and name elements which specify the type of component which should be extended. The component then automatically inherits all properties and type-bound functions of the extended type.

isDefault

[Required] Specifies whether or not this should be the default implementation of this class. Note that only one implementation of each class can be declared to be the default. If no implementation of a given class is declared to be the default then the (automatically generated) null implementation will be made the default.

properties

[Optional] Contains an array of property elements which specify the properties of this implementation. Each member property has the following structure:

name

[Required] The name of the property.

type

[Required] The type (one of real, integer, logical, history, abundances, chemicals, or keplerOrbit at present) of the property.

rank

[Required] The rank of this property (currently 0 for a scalar or 1 for a 1-D array).

attributes

[Required] Attributes of this property:

isSettable

If true then the value of this property can be set directly.

isGettable

If true then the value of this property can be got directly.

isEvolvable

If true this property evolves as part of the Galacticus ODE system.

createIfNeeded

If true then any attempt to get, set, or adjust the rate of this property will cause the component to be created if it does not already exist. This is useful if the component should be created in response to mass transfer from some other component for example.

isDeferred

Contains a “:” separated list which can contain get, set, and rate. The methods present in this list will not have functions bound to them at compile time. Instead a function will be created which allows a function to be bound to these methods at run time. For example:

call myComponent%massFunction    (My_Component_Mass_Get_Function)
call myComponent%massSetFunction (My_Component_Mass_Set_Function)
call myComponent%massRateFunction(My_Component_Mass_Rate_Function)

Additionally, a method is created which returns true or false depending on whether the method has been attached to a function yet, e.g.

myComponent%massIsAttached    ()
myComponent%massSetIsAttached ()
myComponent%massRateIsAttached()

output

[Optional] If present, the property will be included in the Galacticus output file. The following attributes control the details of that output:

unitsInSI

The units of the output quantity in the SI system.

unitsDescription

A human-readable description of the units of the output quantity.

unitsQuantity

An astropy.units-parseable description of the units of the output quantity.

comment

A comment to be included with the HDF5 dataset for this property.

condition

A statement which must evaluate to true or false and which will be used to determine if the property will be output. The present output time for is available as time. In the case of an array property the construct “{i}” can be used to pass the index of the element for which the condition should be evaluated.

modules

A comma-separated list of any modules required to perform the output (e.g. modules which contain functions or values that are used).

Additional attributes are required for array properties:

labels

This can be an array, declared as “[L_1,…,L_N]”, specifying the suffix to be added to the property name for each component of the array in the output, or a function which returns the suffix. In the case of a function the construct “{i}” can be used to pass the index of the element for which the suffix is required.

count

A statement which evaluates the the number of elements to be output (i.e. the length of the array).

isVirtual

[Optional] If present and set to “true”, this property is a virtual property. A virtual property has no data associated with it and must supply its own functions for getting, setting and adjusting its rate of change (if allowed by the property’s attributes). Virtual properties are used for quantities which are derived from actual properties of the component implementation (for example, a star formation rate could be a virtual property if it is derived from an actual gas mass property) or for adjusting the rates of several actual properties simultaneously.

isNonNegative

[Optional] If present and set to “true”, this property should always be non-negative (i.e. zero or positive). This is typically the case for quantities such as masses for example. This attribute is informative only—it may or may not be taking into account by the class responsible for evolving the component properties. For example, the mergerTreeNodeEvolverStandard class will evolve properties marked as non-negative in such a way as to ensure they remain non-negative, but only if its parameter [enforceNonNegativity]=true.

getFunction

[Optional] Specifies the function to be used for getting the value of the property, overriding the default get function. The function must be included in the Galacticus_Nodes module by use of the functions element described below. Note that this function, by virtue of its privileged access to the internal structure of node components, can access the value of the data associated with the property using:

myComponent%<property>Data%value

setFunction

[Optional] The same as getFunction but defines a function to set the value of the property.

classDefault

[Optional] Specifies the default value for this property if the component class has not been created (i.e. has no specific implementation yet). The content of this element gives the default value (which can be a scalar, an array, a function, etc.). Additional, optional attributes control the use of this element:

modules

Specifies a comma-separated list of modules which are required to set the default values (e.g. modules which contain the value or function to be used).

count

For array properties whose size is not known at compile-time, it is possible to specify a function which will return the appropriate size of the array at run-time. The scalar default value given in the classDefault element will then be replicated the appropriate number of times.

bindings

[Optional] Contains an array of binding elements which specify functions to bind to this implementation. Each member binding has the following structure:

method

The name of the bound method, such that the function can be accessed using

myComponent%<method>(...)

function

The function to which the method should be bound. (This function must be included in the Galacticus_Nodes module by use of the functions element described below.

type

The type of function.

functions

[Optional] Contains the name of a file which will be included into the Galacticus_Nodes module. This file can contain functions which will be bound to this implementation. By virtue of being included in the Galacticus_Nodes module these functions have privileged access to the internal structure of all node component objects.

Component Initialization

Initialization of a component module (if necessary, for example, to read parameters or allocate workspace) can occur at a number of different points in the execution of Galacticus. Providing initialization occurs in advance of any calculations then any point is acceptable. One possibility is simply to call an initialization function at the head of all functions defined in the component module. This initialization function should return immediately if it has already been called (to avoid duplicate initialization). Another option is to use a mergerTreeOperator to perform initialization just before merger trees are constructed (the initialization function must again return immediately if it has been previously called).

Optionally, a component may include a mergerTreeEvolveThreadInitialize directive, which gives the name of a subroutine in its unitName element. The routine specified by mergerTreeEvolveThreadInitialize is called by all threads prior to merger tree evolution, and can therefore be used to perform any “per thread” initialization. Note that this routine will be called many times during a given Galacticus run—it is the responsibility of the routine to ensure that it performs any initialization only once.

Component Access, Creation and Destruction

When a node is created, it initially contains no components. A component must therefore create itself on the fly as needed. Typically, a component is first created when an attempt is made to set a property value, or to adjust the rate of change of a property value or in response to some event (e.g. a satellite component may be created in response to a node merging with a larger node). Requests for property values frequently do not require that the component exist, as a zero value can often be returned instead[2].

To access a component from a node, use:

myComponent => thisNode%<class>([instance=<N>,autoCreate=<create>])

where class is the component class required, the optional instance argument requests a specific instance of the component (relevant if the node contains more than one of a particular component, e.g. if it contains two supermassive black holes for example; if no instance is specified the first instance will be returned), and the autoCreate option specifies whether or not the component should be automatically created (assuming it does not already exist). autoCreate\(=\)true should be used to create components initially.

A component of a node can be destroyed using:

call thisNode%<class>Destroy()

Component Implementations

Component implementations optionally provide functions to get and set their properties (and to set the rate of change of evolvable properties) so that other components and functions within Galacticus to can interact with them in a way that is independent of the specific component implementation chosen. To permit this, Galacticus creates functions for each property to access it in all permitted ways. For example, the exponential implementation of the disk component class has a “massStellar” property defined by:

<method>
  <name>massStellar</name>
  <type>real</type>
  <rank>0</rank>
  <attributes isSettable="true" isGettable="true" isEvolvable="true" />
</method>

This causes Galacticus to define several functions bound to the nodeComponentDisk class:

massStellarIsSettable

Returns true if this property is settable;

massStellarIsGettable

Returns true if this property is gettable;

massStellarSet

Sets the value of this property to the supplied argument;

massStellarGet

Gets the value of this property;

massStellarRate

Cumulates its argument to the rate of change of this property;

massStellarScale

Sets the absolute scale for this property used in ODE error control;

along with several others used internally for output, serialization etc.

Component Evolution

All component properties which have an isEvolvable attribute set to true are included in Galacticus’s ODE solver as the node is evolved forward in time. As described in Section Component Implementations, Galacticus will create two functions that permit the rate of change of a property adjusted and for the absolute scale used in ODE error control to be set.

When evolving ODEs the ODE solver aims to keep the error on property \(i\) below

\[D_i = \epsilon_\mathrm{abs} s_i + \epsilon_\mathrm{rel} |y_i|,\]

where \(epsilon_\mathrm{abs}=\)[odeToleranceAbsolute], \(epsilon_\mathrm{rel}=\)[odeToleranceRelative], \(y_i\) is the value of property \(i\) and \(s_i\) is a scaling factor which controls the absolute tolerance for this property. By default, \(s_i=1\), but this can be changed for a component utilizing the scaleSetTask directive. This allows a function to be called in which the component sets suitable scale factors for each of its properties prior to any ODE evolution being carried out. This can be very useful, for example, in cases where two components are coupled. Consider a case where a disk is transferring material to a spheroid via a bar instability. If the disk is orders of magnitude more massive that the spheroid then the rate of mass transfer can be very high (i.e. \(\dot{y}/y\) for the spheroid will be large). With just a relative tolerance (i.e. the \(\epsilon_\mathrm{rel} |y_i|\) term) this would require very short timesteps for the spheroid. However, in such cases we don’t care about such tiny tolerances for the spheroid (since it will grow to be substantially more massive). Therefore, it may be appropriate to set \(s_i\) to be equal to the sum of the disk and spheroid properties for example. The scale set directive and associated subroutine should follow this template:

!![
<scaleSetTask>
  <unitName>Node_Component_Disk_Exponential_Scale_Set</unitName>
</scaleSetTask>
!!]
subroutine Node_Component_Disk_Exponential_Scale_Set(thisNode)
  implicit none
  type (treeNode         ), pointer, intent(inout) :: thisNode
  class(nodeComponentDisk), pointer                :: disk

  ! Get the disk component.
  disk => thisNode%disk()
  ! Check if an exponential disk component exists.
  select type (disk)
  class is (nodeComponentDiskExponential)
    ...
    call disk%massStellarScale(massScale)
    ...
  end select
  return
end subroutine Node_Component_Disk_Exponential_Scale_Set

Sensible choices for the \(s_i\) factors can significantly speed-up execution of Galacticus.

Evolution Interrupts

It is often necessary to interrupt the smooth ODE evolution of a node in Galacticus. This can happen if, for example, a galaxy mergers with another galaxy (in which case the merger must be processed prior to further evolution) or if a component must be created before evolution can continue. The rate adjust and rate compute subroutines allow for interrupts to be flagged via their interrupt and interruptProcedure arguments. If an interrupt is required then interrupt should be set to true, while interruptProcedure should be set to point to a procedure which will handle the interrupt. Then, providing no other interrupt occurred earlier, the evolution will be stopped and the interrupt procedure called before evolution is continued.

An interrupt procedure should have the form:

subroutine My_Interrupt_Procedure(thisNode)
  implicit none
  type(treeNode), pointer, intent(inout) :: thisNode

  ! Do whatever needs to be done to handle the interrupt.

  return
end subroutine My_Interrupt_Procedure

Existing Classes

Function Classes

Functions implement basic calculations (e.g. computing the power spectrum). Additional implementations of a functionClass are added using the a directive with the name of that functionClass. The implementation should be placed in a file containing the directive. For example, for the mergerTreeTimestep functionClass the file must containing a directive of the form

!![
<mergerTreeTimestep name="mergerTreeTimestepMyImplementation">
  <description>A short description of the implementation.</description>
</mergerTreeTimestep>
!!]

where MyImplementation is an appropriate name for the implementation. This file should be treated as a regular Fortran submodule, but without the initial submodule and final end submodule lines. That is, it may contain use statements and variable declarations prior to the contains line, and should contain all functions required by the implementation after that line. It is standard, but not required, that function names should be prefixed with the name of the implementation (e.g. “MyImplementation” in the above example). The file must define a type that extends the named functionClass class (or extends another type which is itself an extension of the named functionClass class), containing any data needed by the implementation along with type-bound functions required by the implementation.

Events

Events are triggered during merger tree evolution. Examples are when a node needs to be promoted to its parent node, or when a minor node merges with its parent.

Node Promotion Events

Additional methods for node promotion (i.e. when a primary progenitor reaches its parent halo) can be added using the nodePromotionTask directive. The directive should contain a single argument, giving the name of a subroutine to be called to initialize the method. For example, the basic tree node method uses this directive as follows:

!![
<nodePromotionTask>
 <unitName>Tree_Node_Basic_Promote</unitName>
</nodePromotionTask>
!!]

Here, Tree_Node_Basic_Promote is the name of a subroutine which will be called to perform whatever tasks are required prior to the promotion. The subroutine must have the following form:

 subroutine Node_Promotion_Task(thisNode)
  implicit none
  type(treeNode), pointer, intent(inout) :: thisNode
  .
  .
  .
  return
end subroutine Node_Promotion_Task

where thisNode is the node about to be promoted.

Tasks

Tasks are any processing which must be performed on a node as a result of some specific event (such as a merger).

Analytic Solvers

Galacticus typically solves the system of ODE which describe the evolution of a galaxy using a numerical solver. However, in some situations analytic solutions may be available. Galacticus allows for the use of analytic solvers in such cases. Currently available analytic solvers are described below.

Simple Disk Satellites

This solver, which can be activated by setting [diskVerySimpleUseAnalyticSolver]\(=\)true is applicable to satellite systems in the very simple disk component. In particular, the following conditions must be met for it to be valid:

• The hot halo component must be of type verySimple or verySimpleDelayed;

• The satellite component must be of type verySimple;

• The disk component must be of type verySimple;

• Properties of all other components must not change for satellite galaxies (with the exception of the time property of the basic component), or else the component must be null;

• Disk star formation and outflow rates must scale as \(M_\mathrm{gas}/\tau\) where \(M_\mathrm{gas}\) is the instantaneous mass of gas in the galaxy disk, and \(\tau\) is a fixed timescale for any given satellite galaxy.

Under these conditions, the flow of mass between gas, stellar, and outflowed phases is analytically solvable.

Halo Formation Events

Tasks to be performed when a halo is deemed to have “formed” (or reformed) can be registered using the haloFormationTask directive. For example, the Tree_Node_Methods_Hot_Halo module registers a task using

!![
<haloFormationTask>
  <unitName>Hot_Halo_Formation_Task</unitName>
</haloFormationTask>
!!]

The contents of <unitName> should give the name of the subroutine to be called on halo formation. The subroutine should have a single argument, thisNode, which is the node that has (re)formed.

HDF5 File Close

Tasks to be performed just prior to closing the Galacticus output HDF5 file (typically involving writing accumulated data to that file) can be registered using the hdfPreCloseTask directive. For example, the Merger_Tree_Timesteps_History module registers a task using

!![
<hdfPreCloseTask>
  <unitName>Merger_Tree_History_Write</unitName>
</hdfPreCloseTask>
!!]

The contents of <unitName> should give the name of the subroutine to be called prior to HDF5 file closure. The subroutine should have no arguments.

Merger Tree Extra Output Tasks

Extra outputs for merger trees (i.e. those which do not involve output of a fixed number of properties for every node—examples might be star formation histories for a subset of galaxies) can be added using the directive: mergerTreeExtraOutputTask. The directive should give the name of the subroutine to be called to perform the task. A template for this task is:

!![
<mergerTreeExtraOutputTask>
 <unitName>Galacticus_Extra_Output_Example</unitName>
</mergerTreeExtraOutputTask>
!!]
subroutine Galacticus_Extra_Output_Example(thisNode,iOutput,treeIndex,nodePassesFilter)
  implicit none
  type(treeNode),          intent(inout), pointer :: thisNode
  integer,                 intent(in)             :: iOutput
  integer(kind=kind_int8), intent(in)             :: treeIndex
  logical,                 intent(in)             :: nodePassesFilter
  .
  .
  .
  return
end subroutine Galacticus_Extra_Output_Example

The subroutine will be called for each node in each merger tree at each output, and should perform whatever extra output related to thisNode. The index of the output and tree are provided as iOutput and treeIndex for reference, and may be used in organizing output. The nodePassesFilter flag will be set to true if thisNode passed all active output filters (see galacticFilter). If it is false then typically no output should occur (although other tasks may still be undertaken).

Merger Tree Output Tasks

Additional outputs for merger trees can be added using three directives: mergerTreeOutputPropertyCount, mergerTreeOutputNames and mergerTreeOutputTask. Each directive should give the name of the subroutine to be called to perform the task and, additionally, a name for sorting (this should be the same for all three directives and ensures that output tasks are always called in the correct order). Templates for these tasks are:

!![
<mergerTreeOutputNames>
 <unitName>Galacticus_Output_Tree_Example_Names</unitName>
 <sortName>Galacticus_Output_Tree_Example</sortName>
</mergerTreeOutputNames>
!!]
subroutine Galacticus_Output_Tree_Example_Names(integerProperty,integerPropertyNames,integerPropertyComments,integerPropertyUnitsSI &
     &,doubleProperty,doublePropertyNames,doublePropertyComments,doublePropertyUnitsSI,time)
  implicit none
  double precision, intent(in)                  :: time
  integer,          intent(inout)               :: integerProperty,doubleProperty
  character(len=*), intent(inout), dimension(:) :: integerPropertyNames,integerPropertyComments,doublePropertyNames &
       &,doublePropertyComments
  double precision, intent(inout), dimension(:) :: integerPropertyUnitsSI,doublePropertyUnitsSI
  .
  .
  .
  return
end subroutine Galacticus_Output_Tree_Example_Names

!![
<mergerTreeOutputPropertyCount>
 <unitName>Galacticus_Output_Tree_Example_Property_Count</unitName>
 <sortName>Galacticus_Output_Tree_Example</sortName>
</mergerTreeOutputPropertyCount>
!!]
subroutine Galacticus_Output_Tree_Example_Property_Count(integerPropertyCount,doublePropertyCount)
  implicit none
  integer, intent(inout) :: integerPropertyCount,doublePropertyCount
  .
  .
  .
  return
end subroutine Galacticus_Output_Tree_Example_Property_Count

!![
<mergerTreeOutputTask>
 <unitName>Galacticus_Output_Tree_Example</unitName>
 <sortName>Galacticus_Output_Tree_Example</sortName>
</mergerTreeOutputTask>
!!]
subroutine Galacticus_Output_Tree_Example(thisNode,integerProperty,integerBufferCount,integerBuffer,doubleProperty&
     &,doubleBufferCount,doubleBuffer)
  implicit none
  type(treeNode),          intent(inout), pointer :: thisNode
  integer,                 intent(inout)          :: integerProperty,integerBufferCount,doubleProperty,doubleBufferCount
  integer(kind=kind_int8), intent(inout)          :: integerBuffer(:,:)
  double precision,        intent(inout)          :: doubleBuffer(:,:)
  .
  .
  .
  return
end subroutine Galacticus_Output_Tree_Example

The mergerTreeOutputPropertyCount subroutine must simply increment integerPropertyCount and doublePropertyCount by the number of integer and double precision properties that will be output respectively. The mergerTreeOutputNames subroutine must store the dataset names, comments and units in the SI system[3] for each integer and double precision property in the supplied arrays. The value of integerProperty and doubleProperty should be incremented by 1 before each property name/comment is set—these then supply the position within the input arrays in which to store the name. The mergerTreeOutputTask subroutine must similarly place the desired property values for thisNode into the supplied arrays. The value of integerProperty and doubleProperty should be incremented by 1 before each property value is set. The value can then be stored in, for example, integerBuffer(integerBufferCount,integerProperty).

Merger Tree Initialization Tasks

Additional tasks to be performed during merger tree initialization can be added using the mergerTreeInitializeTask directive. The directive should contain a single argument, giving the name of a subroutine to be called to initialize the method. For example, the standard basic component method uses this directive as follows:

!![
<mergerTreeInitializeTask>
 <unitName>Halo_Mass_Accretion_Rate</unitName>
</mergerTreeInitializeTask>
!!]

Here, Halo_Mass_Accretion_Rate is the name of a subroutine which will be called to perform whatever tasks are required as a result of the merger. The subroutine must have the following form:

 subroutine Merger_Tree_Initialize_Task(thisNode)
  implicit none
  type(treeNode), pointer, intent(inout) :: thisNode
  .
  .
  .
  return
end subroutine Merger_Tree_Initialize_Task

where thisNode is the node to be initialized. The subroutine will be called once for each node in the tree.

Post-step Tasks

Additional methods for post-step tasks (i.e. things that should be done after each ODE solver step when evolving a node differentially) can be added using the postStepTask directive. The directive should contain a single argument, giving the name of a subroutine to be called to perform the task. For example, the standard hot halo component adds a task as follows:

!![
<postStepTask>
 <unitName>Tree_Node_Hot_Halo_Poststep_Standard</unitName>
</postStepTask>
!!]

Here, Tree_Node_Hot_Halo_Poststep_Standard is the name of a subroutine which will be called to perform whatever tasks are required. The subroutine must have the following form:

 subroutine Post_Step_Task(node,status)
  implicit none
  type   (treeNode), intent(inout) :: node
  integer          , intent(inout) :: status
  .
  .
  .
  return
end subroutine Post_Step_Task

where node is the node for which tasks should be performed. If any change is made to the state of the node then status should be set equal to GSL_Failure. Tasks typically involve cleaning up after differential evolution.

Radius Solver Tasks

Galactic radii solver functions (see galacticStructureSolver) need to be able to interact with the components of a tree node to

1. Determine which components want a radius to be solved for;

2. Get and set the properties of those components.

The radiusSolverPlausibility and radiusSolverTask directives facilitate this. A component which has a radius to be solved for should include directives of the form:

!![
<radiusSolverTask>
 <unitName>Component_Radius_Solver_Plausibility</unitName>
</radiusSolverTask>
 !!]

and

!![
<radiusSolverTask>
 <unitName>Component_Radius_Solver</unitName>
</radiusSolverTask>
 !!]

where Component_Radius_Solver_Plausibility is the name of a subroutine which will specify whether or not the component is physically plausible for radius solving (e.g. has non-negative mass) and should have the following form:

subroutine Component_Radius_Solver_Plausibility(thisNode,galaxyIsPhysicallyPlausible)
   implicit none
   type(treeNode), pointer, intent(inout) :: thisNode
   logical,                 intent(inout) :: galaxyIsPhysicallyPlausible
   .
   .
   .
   return
end subroutine Component_Radius_Solver_Plausibility

which should set galaxyIsPhysicallyPlausible to false if the component is not physically plausible, but should otherwise leave galaxyIsPhysicallyPlausible unchanged. Additionally, Component_Radius_Solver is the name of a subroutine which will supply the necessary information about the node, and which should have the following form:

subroutine Component_Radius_Solver(thisNode,componentActive,specificAngularMomentum,Radius_Get,Radius_Set,Velocity_Get,Velocity_Set)
   implicit none
   type(treeNode),   pointer, intent(inout) :: thisNode
   logical,                   intent(out)   :: componentActive
   double precision,          intent(out)   :: specificAngularMomentum
   procedure(),      pointer, intent(out)   :: Radius_Get,Velocity_Get
   procedure(),      pointer, intent(out)   :: Radius_Set,Velocity_Set
   .
   .
   .
   return
end subroutine Component_Radius_Solver

When called, the subroutine should set componentActive to indicate whether or not this nod contains an active component of the type. If it does, it should also set specificAngularMomentum to reflect the specific angular momentum (in km s\(^{-1}\) Mpc) of the component (at whatever point in its profile the radius is required) and should point the four procedure pointers to routines which get and set the radius and circular velocity properties of the component (which should have the standard form for component get and set methods). It is acceptable for the set procedures to point to dummy routines.

The galactic structure radii solver routines will use this information to determine (and set) the radius and circular velocity of the component. An advantage of this approach is that different radii solver methods can all use this same system, ensuring that just a single interface is needed in each component.

Satellite Host Change Tasks

Additional methods for satellite host change events (i.e. when a satellite node moves to a new host) can be added using the satelliteHostChangeTask directive. The directive should contain a single argument, giving the name of a subroutine to be called to initialize the method. For example, the simple satellite orbits components uses this directive as follows:

!![
<satelliteHostChangeTask>
 <unitName>Satellite_Orbit_New_Host</unitName>
</satelliteHostChangeTask>
!!]

Here, Satellite_Orbit_New_Host is the name of a subroutine which will be called to perform whatever tasks are required as a result of the host change. The subroutine must have the following form:

 subroutine New_Host_Task(thisNode)
  implicit none
  type(treeNode), pointer, intent(inout) :: thisNode
  .
  .
  .
  return
end subroutine New_Host_Task

where thisNode is the node which has changed host (the new host halo is thisNode%parentNode).

Subsystems

This section describes some of the subsystems within Galacticus that support various physical entities or processes.

Kepler Orbits

The keplerOrbit object (provided by the Kepler_Orbits_Structure module) stores the parameters of a single Keplerian orbit. It internally handles computation of additional/alternate orbital parameters once an orbit has been fully defined. Currently, the orientation of orbits (i.e. the unit vector normal to the orbital plane and the argument of periapsis) is not tracked. As such, orbits are fully defined by three parameters (in addition to the masses of the orbiting bodies). The following limitations presently apply to the keplerOrbit object:

• If an orbit is overdefined (i.e. if more than three parameters are set manually) no checking is performed to ensure that the parameters are consistent with a Keplerian orbit;

• Not all interconversions between parameters are supported[4]. If a conversion cannot be performed, an error message will be given.

A keplerOrbit object can be reset by calling the reset() method, and its defined/undefined status can be tested with the isDefined() method or asserted with the assertIsDefined() method. The following orbital parameters are supported, each method returning the value of the parameter and a corresponding method suffixed with Set can be used to set the parameter: radius, velocityRadial, velocityTangential, energy, angularMomentum, eccentricity, semiMajorAxis, radiusPericenter, radiusApocenter. Additionally, the masses of the orbiting bodies are provided by the hostMass() and reducedMassSpecific()\(=M_\mathrm{host}/(M_\mathrm{host}+M_\mathrm{satellite})\) methods. Finally, the velocityScale() method returns \(\mathrm{G}M_\mathrm{host}/r\) where \(r\) is the radius of the orbit.

Chemicals

The chemicals subsystem provides both a interface to a database of known chemicals (allowing their physical properties to be queried) and a structure to store abundances/masses/etc. of the set of chemicals being tracked in Galacticus. The name “chemicals” is used to denote any chemical species that might be involved in reactions, including molecules, atoms, atomic and molecular ions and electrons.

Chemical Database

The file data/Chemical_Database.cml contains a database of chemicals that can currently be used by Galacticus. It uses a simplified version of the Chemical Markup Language to describe chemicals. An excerpt from the database is shown below:

<list>
 <chemical>
  <id>MolecularHydrogenAnion</id>
  <formalCharge>-1</formalCharge>
  <atomArray>
    <atom>
     <id>1</id>
     <elementType>H</elementType>
    </atom>
    <atom>
     <id>2</id>
     <elementType>H</elementType>
    </atom>
  </atomArray>
  <bondArray>
    <bond>
     <atomRefs2>1 2</atomRefs2>
     <order>1</order>
    </bond>
  </bondArray>
 </chemical>
 .
 .
 .
</list>

The database contains a list of chemicals, each contained within a chemical element. The id element provides a label for the chemical (usually a descriptive name with no white space). The formalCharge element gives the charge of the chemical in units of the elementary charge. The chemical is then describe by a list of atoms and bonds inside atomArray and bondArray elements respectively. The atomArray can contain any number of atom elements, which should describe each atom in the chemical giving it a unique id number and an elementType, which is the short one or two letter label for the element (e.g. H, Ni, etc.). The bondArray should contain a bond entry for each atomic bond, which itself contains a atomRefs2 element giving the IDs of the two atoms participating in the bond and an order element which gives the order of the bond (e.g. “1” for a single bond).

Chemical Structure

Within Galacticus a chemical is represented using the chemicalStructure type which is provided by the Chemical_Structures module. A chemicalStructure object can be assigned a particular chemical by retrieving that chemical from the database using:

call myChemical%retrieve("chemicalID")

where chemicalID is the ID of the chemical in the database. Any chemical can be exported to a CML file using

call myChemical%export(fileName)

where fileName gives the name of the file to which to export.

Once assigned a chemical, basic properties such as mass and charge (in atomic units) can be accessed using myChemical%mass and myChemical%charge respectively. The mass is computed from the known atomic masses of the constituent atoms of the chemical.

Chemical Abundances

Within Galacticus a set of abundances (or masses, or densities…) for all chemicals being tracked, as specified by the [chemicalsToTrack] input parameter, is stored within a chemicalAbundancesStructure type, as provided by the Chemical_Abundances_Structure module. The structure provides interfaces for setting and retrieving the abundance of a given chemical species, to pack/unpack all chemicals to/from an array, to convert from mass-weighted to number-weighted quantities and to multiply and divide the chemicals abundances by a given amount. Additionally, the Chemical_Abundances_Structure module provides functions which provide a count of the number of chemicals tracked, to look up the index of a chemical array representation from its name, and to retrieve the name of a given chemical.

Coordinates

The coordinate class describes a position in three-dimensional space. Each extension of this class (currently, coordinateCartesian, coordinateCylindrical, and coordinateSpherical) supply methods to convert to and from Cartesian coordinates. The assignment operator (=) is overloaded such that coordinate objects of any class can be assigned to any other class and conversion to the appropriate coordinate system will happen automatically. A function accepting a coordinate object can therefore convert it to, for example, spherical coordinates simply using

class(coordinate         ), intent(in) :: coordinates
type (coordinateSpherical)             :: coordinatesSpherical
coordinatesSpherical=coordinates

and thereby allow a position to be passed to it in any coordinate system.

[1]

Data objects in components can be real, integer, boolean or of derived type. For derived types, currently history, abundances, chemicals, and keplerOrbit are supported. Adding additional derived types is possible, providing that the type supports the required methods for output, serialization, etc. Data objects can currently be scalar or rank-1 arrays.

[2]

Or some other value if a classDefault has been specified (see Section Implementing a New Component).

[3]

For dimensionless quantities, the units may be set to zero. In such cases, the unitsInSI attribute for the dataset will not be written to the Galacticus output file.

[4]

The keplerOrbit object works by trying to convert to the combination radius, radial and tangential velocities. Once these are defined, all other parameters can be computed. However, for orbits defined in terms of other parameters, the keplerOrbit object does not know how to convert from every such combination of parameters.