TARDIS is a Monte Carlo radiative-transfer spectral synthesis code for 1D models of supernova ejecta. It is designed for rapid spectral modelling of supernovae. The code is described in this documentation and Kerzendorf & Sim 2014.
Note
This documentation is currently under construction and does not describe all of the modes of operations available for TARDIS.
Warning
Currently TARDIS only works on 64-bit python installations. We’re working on making it work on 32-bit python distributions.
TARDIS has the following requirements:
numpy==1.10.0
scipy==0.17.0
pandas==0.16.2
pytables==3.2.2
h5py==2.5
matplotlib==1.4.3
astropy==1.0.5
PyYAML==3.11
numexpr==2.4.4
Cython==0.21
jinja2
networkx==1.10
pytest==2.6.3
Most of these requirements are easy to install using package managers like OS X’s macports or linux package managers or using the Anaconda python distribution.
TARDIS is using astropy’s excellent installation helpers and thus uses similar instructions to astropy.
If you are interested in doing some development for TARDIS please read Developer Workflow.
Once you are done you can run the simple example listed in Running TARDIS.
We highly recommend using the Anaconda python environment to install TARDIS (or
any other scientific packages for that matter). Anaconda has the advantage of
being an isolated environment that can be set to be the default one, but
by no means will mess with your other environments. It will also work on
computers where root-rights are not available. Use these
instructions to install
Anaconda on your machine. The next step is to create an environment for tardis
that contains all of the necessary packages (this ensures that TARDIS
requirements won’t clash with any other python installs on disc:
conda create -n tardis --file https://raw.githubusercontent.com/tardis-sn/tardis/master/conda-requirements python=2
This command fails on some systems. If that is the case, you can download the file manually from the official repository or, if you have a local copy of the tardis repository, use the local file. The command in that case would be:
conda create -n tardis --file conda-requirements python=2
Then to activate this environment simply do:
source activate tardis
and after you are done with TARDIS you can deactivate:
source deactivate
One does not need to recreate the environment, but simply activate it everytime TARDIS is used.
To install the latest stable version of TARDIS simply do:
pip install tardis-sn
or to use the development version:
pip install git+https://github.com/tardis-sn/tardis
A virtual environment is python’s way to ensure that the versions of third-party libraries that TARDIS requires do not interfere with the system-wide installation. This is also the way that the majority of core developers for TARDIS operate.
It is nevertheless recommended to install a number of python packages using the system packagemanager. This ensures that third-party non-python libraries like libhdf5, lapack, etc. are installed.
For OS X we recommend the macports package manager:
sudo port install python27 py27-scipy py27-numpy py27-virtualenv py27-astropy py27-h5py py27-yaml py27-pandas py27-pip py27-tables
In Ubuntu 14.04 the pre-requesite install works with this:
sudo apt-get install python-virtualenv python-numpy python-pandas python-scipy python-h5py python-yaml ipython python-matplotlib cython git
First use the following guide to install virtualenv and the TARDIS requirements Python Environment for TARDIS.
After the virtualenv and the requirements are installed, there are a few options of how to proceed.
If one just wants to use TARDIS you can install the latest stable version:
pip install tardis-sn
the latest development version can be installed using:
pip install git+https://github.com/tardis-sn/tardis
If you are interested in doing some development for TARDIS please read Developer Workflow.
Once you are done you can run the simple example listed in Running TARDIS.
We use a clean install of Ubuntu 14.04 as one of our testing grounds. Here’s how we get TARDIS to run:
sudo apt-get install python-dev python-pip python-numpy python-scipy python-h5py python-pandas python-yaml
We now need to install the newest astropy and we will install everything into our users directory:
pip install astropy --user
Once astropy is installed, install TARDIS:
pip install tardis-sn
Note
On a clean install of Mountain Lion, here’s how we get TARDIS running:
First install macports
Use macports to install:
sudo port install python27 py27-astropy py27-h5py py27-yaml py27-pandas py27-pip
Then install TARDIS:
pip-2.7 install tardis-sn --user --pre
Before running, ensure that the directory ~/Library/Python/2.7/bin is in the appropriate path.
Note
This has also been successfully tested on a clean MAC OS 10.9.1 (Mavericks) install.
We highly encourage with any installation problems to try the recommended install method because this often fix problems. Here are some common problems when installing and their fixes:
Problem: While building tardis via python setup.py build you
may encounter the following error:
error: tardis/montecarlo/montecarlo.c: Could not find C file tardis/montecarlo/montecarlo.c for Cython file tardis/montecarlo/montecarlo.pyx when building extension tardis.montecarlo.montecarlo. Cython must be installed to build from a git checkout.
Solution: There are several solutions to this problem. A clean checkout will
help. To clean up your repository please try python setup.py clean and
then git clean -dfx (WARNING will delete any non tardis file in that directory)
This will often clean this problem. If it still persists:
Go into the tardis/montecarlo directory and build montecarlo.c by hand:
cython montecarlo.pyx
Then, python setup.py build should run without problems.
To run TARDIS requires two files. The atomic database (for more info refer to Atomic Data for TARDIS) and a configuration file (more info at Configuration File).
After installing TARDIS just download the example directory https://www.dropbox.com/s/svvyr5i7m8ouzdt/tardis_example.tar.gz and run TARDIS with:
tar zxvf tardis_example.tar.gz
cd tardis_example
tardis tardis_example.yml output_spectrum.dat
Then plot the output_spectrum.dat with your favourite plotting program. Here’s an example how to do this with python. (The only thing you need to install is ipython and matplotlib - in addition to TARDIS’s requirements)
ipython --pylab
>>> tardis_spec = loadtxt('output_spectrum.dat')
>>> plot(tardis_spec[:,0], tardis_spec[:,1])
More atomic datasets can be downloaded from Atomic Data for TARDIS.
To get more information from each run of TARDIS one can run it interactively and have full access to the model properties
>>> from tardis import run_tardis
>>> model = run_tardis('myconfig.yml')
This model can then be used for inspecting the run as described Accessing Physical Quantities
To get a detailed explanation on gui layout go to Graphical User Interface .
To setup and run the GUI(under development) follow these steps:
The gui can use one of two python bindings for qt, namely PyQt4 and PySide. You can choose which binding is used by setting the environment variable QT_API in your bash.
1. Installing required packages
source activate tardis
conda install ipython=3.0.0 pyside=1.2.1 shiboken=1.2.1
2. Choosing between PySide and PyQt4
#To choose PySide
export QT_API=pyside
#To choose PyQt
export QT_API=pyqt
3. An example of creating a model and GUI
To show the gui from the ipython shell use the following commands.
ipython --pylab=qt4
>>> from tardis import run_tardis
>>> mdl = run_tardis('yamlconfigfile.yml', 'atomdatafile.h5')
>>> from tardis.gui import interface
>>> interface.show(mdl)
If you just want to run from a configuration file and show the results, you can do that outside the ipython shell.To do this navigate to the folder where you installed tardis and go to tardis/tardis/gui, and use the following command.
python interface.py path-to-yaml-configuration-file path-to-atomic-data-file
TARDIS is designed to carry out calculations of synthetic supernova spectra (see Kerzendorf & Sim 2014).
A TARDIS calculation requires
These inputs are specified via the input file (yaml file). As a starting point, examine the input file in the example download (above). The example file sets up a very simple model in which the density distribution is based on the well-known W7 model (Nomoto et al. 1984; fit following Branch et al. 1985) and in which uniform abundances are specified. To get you started, here is a short list of some settings in that example that can be experimented with to get a feel for using TARDIS. Please contact us if you have questions or would like information on more advanced possibilities!
In the example file, you can alter the elemental abundances (mass fractions; specified in the “abundances” section). Use this to explore the sensitivity of features to composition.
You can alter the (emergent) luminosity and see how this affects the synthetic spectrum. Do this by varying the “luminosity requested” in the “supernova” section.
In the example file, you can change the epoch for which the synthetic spectrum is calculated (change the value of “time explosion”; specified in the “supernova” section). When doing this you might also change the inner boundary velocity (“start” value in the “velocity” section), and probably the luminosity (see above)!
In the “plasma” section you can change the “line interaction type” between “scatter”, “downbranch” and “macroatom” - investigate how important fluorescence is in your synthetic spectrum. (See Kerzendorf & Sim and references therein for the meaning of these settings.)
Despite our best efforts TARDIS contains bugs. We want to make sure that fixing these is quick and painless for users and developers.
When reporting bugs please use the mailing list (tardis-sn-user) or create an issue for TARDIS on github (??? ;preferred).
If there are very long ouputs needed to debug the problem, please use a website like pastebin that way it is easier for us to look through it.
One of the major problems that we see are mismatches in installed versions of third-party libraries (that’s where virtualenvs really help - see Python Environment for TARDIS). To make sure that we can test this properly please tell us the installed version of third party packages. The easiest way is:
$ pip freeze
Cython==0.20.2
Jinja2==2.7.2
MarkupSafe==0.21
PyYAML==3.11
Pygments==1.6
Sphinx==1.2.2
astropy==1.0.dev10065
astropy-helpers==0.4.2
backports.ssl-match-hostname==3.4.0.2
cov-core==1.11
coverage==3.7.1
docutils==0.11
execnet==1.2.0
h5py==2.4.0a0
ipdb==0.8
ipython==1.2.1
latexcodec==0.3.2
matplotlib==1.4.0
mock==1.0.1
nose==1.3.4
numexpr==2.4
numpy==1.9.0
numpydoc==0.4
oset==0.1.3
pandas==0.14.1
pep8==1.5.6
py==1.4.20
pybtex==0.17
pybtex-docutils==0.2.0
pyparsing==2.0.2
pytest==2.5.2
pytest-cache==1.0
pytest-cov==1.6
pytest-ipdb==0.1-prerelease
pytest-pep8==1.0.6
python-dateutil==2.2
pytz==2014.7
pyzmq==14.2.0
scipy==0.13.3
six==1.8.0
sphinx-bootstrap-theme==0.4.0
sphinxcontrib-bibtex==0.2.9
sphinxcontrib-tikz==0.4.1
tables==3.1.1
tornado==3.2
tox==1.7.1
virtualenv==1.11.6
wsgiref==0.1.2
and pasting the output in your error report, this will show us what packages are installed and used.
The next thing that is useful to do if TARDIS is downloaded in a directory is to run the test suites and find if any of the tests fail:
python setup.py test
running test
running build
running build_py
copying tardis/./macro_atom.c -> build/lib.macosx-10.9-x86_64-2.7/tardis/.
copying tardis/./montecarlo.c -> build/lib.macosx-10.9-x86_64-2.7/tardis/.
running build_ext
skipping 'tardis/montecarlo.c' Cython extension (up-to-date)
skipping 'tardis/macro_atom.c' Cython extension (up-to-date)
running build_scripts
============================= test session starts ==============================
platform darwin -- Python 2.7.8 -- pytest-2.5.1
Running tests in tardis /Users/wkerzend/scripts/python/tardis/docs.
Platform: Darwin-13.4.0-x86_64-i386-64bit
Executable: /Users/wkerzend/.virtualenvs/tardis-devel/bin/python
Full Python Version:
2.7.8 (default, Oct 15 2014, 22:04:42)
[GCC 4.2.1 Compatible Apple LLVM 5.1 (clang-503.0.40)]
encodings: sys: ascii, locale: UTF-8, filesystem: utf-8, unicode bits: 15
byteorder: little
float info: dig: 15, mant_dig: 15
Numpy: 1.9.1
Scipy: 0.13.3
astropy: 1.0.dev10065
yaml: 3.11
cython: 0.20.2
h5py: 2.4.0a0
Matplotlib: 1.4.0
ipython: 1.2.1
plugins: cache, cov, ipdb, pep8
collected 107 items
tardis/io/tests/test_ascii_readers.py .......
tardis/io/tests/test_config_reader.py .............................
tardis/io/tests/test_config_validator.py ........................
tardis/io/tests/test_configuration_namespace.py .........
tardis/tests/test_atomic.py .....s
tardis/tests/test_lte_plasma.py ssssssssss
tardis/tests/test_plasma_nlte.py .....
tardis/tests/test_plasma_simple.py .
tardis/tests/test_tardis_full.py s
tardis/tests/test_util.py ...............
==================== 95 passed, 12 skipped in 8.53 seconds =====================
This will hopefully help us to identify us the problem. We will continue to update this document with other techniques.
Most of the configuration for TARDIS uses the YAML format. If a YAML file is handed to TARDIS it is first checked using a configuration validator to see if the required keywords of the required type exist and produces a dictionary with that has the validated types as well as sensible defaults for missing values built-in. This validator uses a default validation file that already explains a lot about the required configuration items.
tardis_config_version:
property_type: string
default: None
mandatory: True
help: Version of the configuration file
supernova:
luminosity_requested:
property_type: quantity
mandatory: True
default: 1 solLum
help: requested output luminosity for simulation
time_explosion:
property_type: quantity
mandatory: True
default: None
help: time since explosion
distance:
property_type: quantity
mandatory: False
default: None
help: distance to supernova
luminosity_wavelength_start:
property_type: quantity
mandatory: False
default: 0 angstrom
help: start of the integral needed for getting the luminosity right
luminosity_wavelength_end:
property_type: quantity
mandatory: False
default: inf angstrom
help: start of the integral needed for getting the luminosity right
atom_data:
property_type: string
mandatory: True
help: path or filename to the Atomic Data HDF5 file
plasma:
initial_t_inner:
property_type: quantity
mandatory: False
default: -1 K
help: >
initial temperature of the inner boundary black body. If set to -1 K
will result in automatic calculation of boundary
initial_t_rad:
property_type: quantity
mandatory: False
default: -1 K
help: >
initial radiative temperature in all cells. If set to -1 K will result
in automtatic calculation of the initial temperatures
disable_electron_scattering:
property_type: bool
mandatory: False
default: False
help: >
disable electron scattering process in montecarlo part - non-physical only
for tests
ionization:
property_type: string
mandatory: True
default: None
allowed_value: nebular lte
help: ionization treatment mode
excitation:
property_type: string
mandatory: True
default: None
allowed_value: lte dilute-lte
help: excitation treatment mode
radiative_rates_type:
property_type: string
mandatory: True
default: None
allowed_value: dilute-blackbody detailed blackbody
help: radiative rates treatment mode
line_interaction_type:
property_type: string
mandatory: True
default: None
allowed_value: scatter downbranch macroatom
help: line interaction mode
w_epsilon:
property_type: float
mandatory: False
default: 1e-10
help: w to use when j_blues get numerically 0. - avoids numerical complications
delta_treatment:
property_type: float
mandatory: False
default: None
help: >
In the saha calculation set delta equals to the number given in
this configuration item. if set to None (default), normal delta
treatment (as described in Mazzali & Lucy 1993) will be applied
nlte:
species:
property_type: list
mandatory: False
default: []
help: >
Species that are requested to be NLTE treated in the format
['Si 2', 'Ca 1', etc.]
coronal_approximation:
property_type: bool
default: False
mandatory: False
help: set all jblues=0.0
classical_nebular:
property_type: bool
default: False
mandatory: False
help: sets all beta_sobolevs to 1
helium_treatment:
property_type: string
mandatory: False
default: none
allowed_value: none recomb-nlte numerical-nlte
help: none to treat He as the other elements. recomb-nlte to treat with NLTE approximation.
heating_rate_data_file:
property_type: string
mandatory: False
default: none
help: Path to file containing heating rate/light curve data.
model:
structure:
property_type : container-property
type:
property_type: container-declaration
containers: ['file', 'specific']
_file: ['filename','filetype']
+file: ['v_inner_boundary','v_outer_boundary']
_specific: ['velocity', 'density']
filename:
property_type: string
default: None
mandatory: True
help: file name (with path) to structure model
filetype:
property_type: string
default: None
allowed_value: ['simple_ascii', 'artis']
mandatory: True
help: file type
v_inner_boundary:
property_type: quantity
default: 0 km/s
mandatory: False
help: location of the inner boundary chosen from the model
v_outer_boundary:
property_type: quantity
default: inf km/s
mandatory: False
help: location of the inner boundary chosen from the model
velocity:
property_type: quantity_range_sampled
default: None
mandatory: True
help: description of the boundaries of the shells
density:
property_type: container-property
type:
property_type: container-declaration
containers: ['branch85_w7','exponential','power_law', 'uniform']
_uniform: ['value']
_branch85_w7: []
+branch85_w7: ['w7_time_0', 'w7_rho_0', 'w7_v_0']
_power_law: ['time_0', 'rho_0', 'v_0', 'exponent']
_exponential: ['time_0', 'rho_0', 'v_0']
w7_time_0:
property_type: quantity
default: 0.000231481 day
mandatory: False
help: This needs no change - DO NOT TOUCH
w7_rho_0:
property_type: quantity
default: 3e29 g/cm^3
mandatory: False
help: This needs no change - DO NOT TOUCH
w7_v_0:
property_type: quantity
default: 1 km/s
mandatory: False
help: This needs no change - DO NOT TOUCH
time_0:
property_type: quantity
default: None
mandatory: True
help: Time at which the pure model densities are right
rho_0:
property_type: quantity
default: None
mandatory: True
help: density at time_0
v_0:
property_type: quantity
default: None
mandatory: True
help: at what velocity the density rho_0 applies
exponent:
property_type: float
default: None
mandatory: True
help: exponent for exponential density profile
value:
property_type: quantity
default: None
mandatory: True
help: value for uniform density
abundances:
property_type: container-property
type:
property_type: container-declaration
containers: ['file','uniform']
_uniform: []
_file: ['filetype','filename']
filename:
property_type: string
default: None
mandatory: True
help: filename
filetype:
property_type: string
default: None
mandatory: False
help: type of abundance file to read in
montecarlo:
nthreads:
property_type: int
default: 1
mandatory: False
help: The number of OpenMP threads.
seed:
property_type: int
default: 23111963
mandatory: False
help: Seed for the random number generator
no_of_packets:
property_type: int
default: None
mandatory: True
help: Seed for the random number generator
iterations:
property_type: int
default: None
mandatory: True
help: Number of maximum iterations
black_body_sampling:
property_type: quantity_range_sampled
default: [50 angstrom, 200000 angstrom, 1000000]
mandatory: False
help: Sampling of the black-body for energy packet creation (giving maximum and minimum packet frequency)
last_no_of_packets:
property_type: int
default: -1
mandatory: False
help: >
This can set the number of packets for the last run.
If set negative it will remain the same as all other runs.
no_of_virtual_packets:
property_type: int
default: 0
mandatory: False
help: Setting the number of virtual packets for the last iteration.
virtual_spectrum_range:
property_type: quantity_range_sampled
default: [50 angstrom, 250000 angstrom, 1000000]
mandatory: False
help: Limits of virtual packet spectrum (giving maximum and minimum packet frequency)
enable_reflective_inner_boundary:
property_type: bool
default: False
mandatory: False
help: >
experimental feature to enable a reflective boundary.
inner_boundary_albedo:
property_type: float
default: 0.0
mandatory: False
help: albedo of the reflective boundary
convergence_strategy:
property_type : container-property
type:
property_type: container-declaration
containers: ['damped', 'specific']
_damped: []
+damped: ['damping_constant', 't_inner', 't_rad', 'w',
'lock_t_inner_cycles', 't_inner_update_exponent',]
_specific: ['threshold', 'fraction', 'hold_iterations']
+specific: ['t_inner', 't_rad', 'w', 'lock_t_inner_cycles',
'damping_constant', 't_inner_update_exponent']
t_inner_update_exponent:
property_type: float
default: -0.5
mandatory: False
help: L=4*pi*r**2*T^y
lock_t_inner_cycles:
property_type: int
mandatory: False
default: 1
help: >
The number of cycles to lock the update of the inner boundary temperature.
This process helps with convergence. The default is to switch it off (1 cycle)
hold_iterations:
property_type: int
default: 3
mandatory: True
help: >
the number of iterations that the convergence criteria need to be
fulfilled before TARDIS accepts the simulation as converged
fraction:
property_type: float
default: 0.8
mandatory: True
help: >
the fraction of shells that have to converge to the given
convergence threshold. For example, 0.8 means that 80% of shells
have to converge to the threshold that convergence is established
damping_constant:
property_type: float
mandatory: False
default: 0.5
help: damping constant
threshold:
property_type: float
mandatory: True
help: >
specifies the threshold that is taken as convergence
(i.e. 0.05 means that the value does not change more than 5%)
t_inner:
damping_constant:
property_type: float
mandatory: False
default: 0.5
help: damping constant
threshold:
property_type: float
mandatory: False
help: >
specifies the threshold that is taken as convergence
(i.e. 0.05 means that the value does not change more than 5%)
t_rad:
damping_constant:
property_type: float
mandatory: False
default: 0.5
help: damping constant
threshold:
property_type: float
mandatory: True
help: >
specifies the threshold that is taken as convergence
(i.e. 0.05 means that the value does not change more than 5%)
w:
damping_constant:
property_type: float
mandatory: False
default: 0.5
help: damping constant
threshold:
property_type: float
mandatory: True
help: >
specifies the threshold that is taken as convergence
(i.e. 0.05 means that the value does not change more than 5%)
spectrum:
property_type: quantity_range_sampled
default: None
mandatory: True
help: Final spectrum sampling
The default config validator takes a user configuration and a default configuration and creates a consistent and valid configuration for tardis based on the constraints given in the default configuration. Both input data are normally given as a yaml dictionary with a consistent hierarchical structure i.e. for every item in the user configuration there has to be a declaration in the default configuration at the same hierarchical level. This declaration can be either an unspecific empty level declaration like: - Main_level:
- Second_level:
- Third_level:
…
or a declaration of a configuration item like: - item:
- property_type: int
- default: 1
- mandatory: True
- help: ‘This is a doc string.’
This contains always the keywords help, default, mandatory, and property_type. The keyword help is a doc-string which describes the corresponding item. Default specifies the default value which is used in case that no value for this item is specified in the corresponding user configuration item. If the keyword mandatory is True, the item has to be specified in the user configuration. The keyword property_type is used to specify the type of the item. At the moment, the config parser knows the following types: Int: The property type int is for integer like config items. Float: The property type float is for float like config items. String: The property type string is for string like config items. Quantity: The property type quantity is for physical quantities with units given as string. The string contains value and unit separated by a whitespace E.g. 2 cm. Range: The property type range specifies a range via start and end. Note: abs(start - end ) > 0 Quantity_range: Like property type range but with quantities as start and stop. The consistency of the units is checked. Additionally to the four standard keywords the types integer, float, and quantity can have the keywords allowed_value and allowed_type. allowed_value specifies the allowed values in a list, whereas allowed_type specifies a range of allowed values like “x>10”.
For more complex configurations with dependencies, you can use the containers which allow branching in the configuration. A container is declared in the default configuration file by setting the property_type to container property and specifying the properties of the container with keyword type. The property_type of this section is container-declaration which allows you to specify the possible container items with the keyword container. For every specified container item, the code expects the declaration of all sub items. The keywords for this are “_“ + “name of the container item”. If the type declaration for this container is finished you can specify all container items like normal items. Here is an example for a container configuration with two branches
..source: yaml - container_example:
property_type: container-property
- type:
- property_type: container-declaration
- containers: [‘one’, ‘two’, ‘three’]
- _one: [‘one_one’, ‘one_two’]
- _two: [‘two_one’]
- one_one:
- property_type: string
- default: ‘This is a container item’
- mandatory: False
- help: This is a container item from the container one.
- one_two:
- sub_one_two_one:
- property_type: string
- default: ‘This is a container item’
- mandatory: False
- help: This is a container item from the container one.
- sub_one_two_two:
- property_type: string
- default: ‘This is a container item’
- mandatory: False
- help: This is a container item from the container one.
- two_one:
- quantity_range:
- property_type: quantity_range
- default: [1 m,10 cm] #[Start,End]
- mandatory: False
- help: Like property type range but with quantities as start and stop. The consistency of the units is checked.
To use the default parser create a new config object form the class ConfigurationValidator by either from a dictionaries or from yaml files:
My_config = ConfigurationValidator(default configuration dictionary, user configuration dictionary)
or - My_config = ConfigurationValidator.from_yaml(default configuration file, user configuration file) To access the configuration for tardis use the method get_config
Here’s an overview of some of the different modes of operation and models for TARDIS.
TARDIS can be used to compute a synthetic spectrum for a model with a user-specified density and chosen set of abundances. One simple (but flexible) class of model is when the density is parameterised as an exponential function of velocity, as described below.
In this mode, the density profile (function of velocity and time since explosion) is assumed to follow a functional form:
defined by reference density, velocity and time parameters. These parameters are set in the input yaml file, specifically in the “structure” subsection of the “model” section, under the “density” heading (see example below).
For a model with an exponential density profile, a set of uniform abundances can be supplied directly in the input (yaml) file. Elemental abundances are set in the “abundances” subsection of the “model” section, following the “type: uniform” specifier (see example input file below). They are specified as mass fractions. E.g.
Si: 0.6
S: 0.4
will set the mass fraction of silicon (Z=14) to 0.6 and sulphur (Z=16) to 0.4.
Note
The mass fractions must sum to one. If mass fractions are supplied that do not sum to one, TARDIS will renormalise all the supplied abundances and print a “WARNING” message.
Here is an example of an input file that sets up an exponential density profile with a uniform set of abundances:
tardis_config_version: v1.0
supernova:
luminosity_requested: 9.44 log_lsun
time_explosion: 13 day
atom_data: kurucz_atom_pure.h5
model:
structure:
type: specific
velocity:
start : 1.1e4 km/s
stop : 20000 km/s
num: 20
density:
type : exponential
time_0: 2. day
rho_0: 6.e-10 g/cm^3
v_0: 3000. km/s
exponent: 1.0
abundances:
type: uniform
O: 0.19
Mg: 0.03
Si: 0.52
S: 0.19
Ar: 0.04
Ca: 0.03
plasma:
ionization: nebular
excitation: dilute-lte
radiative_rates_type: dilute-blackbody
line_interaction_type: scatter
montecarlo:
seed: 23111963
no_of_packets : 2.0e+5
iterations: 30
last_no_of_packets: 5.0e+5
no_of_virtual_packets: 5
spectrum:
start: 500 angstrom
stop: 20000 angstrom
num: 10000
TARDIS can be used to compute a synthetic spectrum for a model with a user-specified density and chosen set of abundances. One simple (but flexible) class of model is when the density is parameterised as an power law function of velocity, as described below.
In this mode, the density profile (function of velocity and time since explosion) is assumed to follow a functional form:
This form is defined by reference density, velocity and time parameters, and the “exponent”, each of which is set in the input file (“structure” subsection of the “model” section, under the “density” heading; see example below).
For a model with an power law density profile, a set of uniform abundances can be supplied directly in the input (yaml) file. Elemental abundances are set in the “abundances” subsection of the “model” section, following the “type: uniform” specifier (see example input file below). They are specified as mass fractions. E.g.
Si: 0.6
S: 0.4
will set the mass fraction of silicon (Z=14) to 0.6 and sulphur (Z=16) to 0.4.
Note
The mass fractions must sum to one. If mass fractions are supplied that do not sum to one, TARDIS will renormalise all the supplied abundances and print a “WARNING” message.
Here is an example of an input file that sets up a power law density profile with a uniform set of abundances:
tardis_config_version: v1.0
supernova:
luminosity_requested: 9.44 log_lsun
time_explosion: 13 day
atom_data: kurucz_atom_pure.h5
model:
structure:
type: specific
velocity:
start : 1.1e4 km/s
stop : 20000 km/s
num: 20
density:
type : power_law
time_0: 2 day
rho_0: 1e-11 g/cm^3
v_0: 11000 km/s
exponent: -5
abundances:
type: uniform
O: 0.19
Mg: 0.03
Si: 0.52
S: 0.19
Ar: 0.04
Ca: 0.03
plasma:
ionization: nebular
excitation: dilute-lte
radiative_rates_type: dilute-blackbody
line_interaction_type: scatter
montecarlo:
seed: 23111963
no_of_packets : 2.0e+5
iterations: 30
last_no_of_packets: 5.0e+5
no_of_virtual_packets: 5
spectrum:
start: 500 angstrom
stop: 20000 angstrom
num: 10000
TARDIS can be used to compute a synthetic spectrum for a model with a user-specified density and chosen set of abundances - this makes it possible to experiment with 1D density profiles and simple specifications of the composition
The density profile is supplied in the form of a simple ascii file that should look something like this:
1 day
# index velocity (km/s) density (g/cm^3)
0 9000.0000 5.4869692e-10
1 10500.000 2.1759646e-10
2 12000.000 9.7656229e-11
3 13500.000 4.8170911e-11
4 15000.000 2.5600000e-11
5 16500.000 1.4450533e-11
6 18000.000 8.5733893e-12
7 19500.000 5.3037103e-12
8 21000.000 3.3999447e-12
9 22500.000 2.2474623e-12
In this file:
TARDIS will use this table of density versus velocity to specify the density distribution in the ejecta. For the calculation, TARDIS will use the reference time given in the file to scale the mass densities to whatever epoch is requested by assuming homologous expansion:
The values in the example here define a density profile that is dropping off with
Note
The grid of points specified in the input file is interpreted by TARDIS as defining a grid in which the tabulated velocities are taken as the outer boundaries of grid cells and the density is assumed to be uniform with each cell.
Warning
The example given here is to show the format only. It is not a realistic model. In any real calculation better resolution (i.e. more grid points) should be used.
For a model with density profile supplied via a file (see above), uniform abundances can be supplied directly in the input (yaml) file. Elemental abundances are set in the “abundances” subsection of the “model” section, following the “type: uniform” specifier (see example input file below). They are specified as mass fractions. E.g.
Si: 0.6
S: 0.4
will set the mass fraction of silicon (Z=14) to 0.6 and sulphur (Z=16) to 0.4.
Note
The mass fractions must sum to one. If mass fractions are supplied that do not sum to one, TARDIS will renormalise all the supplied abundances and print a “WARNING” message.
If you create a correctly formatted density profile file (called “density.dat”), here is an example of how it can be combined with a uniform set of abundances:
#New configuration for TARDIS based on YAML
#IMPORTANT any pure floats need to have a +/- after the e e.g. 2e+5
#Hopefully pyyaml will fix this soon.
---
#Currently only simple1d is allowed
tardis_config_version: v1.0
supernova:
luminosity_requested: 9.44 log_lsun
time_explosion: 13 day
atom_data: kurucz_atom_pure_simple.h5
model:
structure:
type: file
filename: density.dat
filetype: simple_ascii
v_inner_boundary: 11000 km/s
v_outer_boundary: 20000 km/s
abundances:
type: uniform
O: 0.19
Mg: 0.03
Si: 0.52
S: 0.19
Ar: 0.04
Ca: 0.03
plasma:
disable_electron_scattering: no
ionization: lte
excitation: lte
radiative_rates_type: dilute-blackbody
line_interaction_type: macroatom
montecarlo:
seed: 23111963
no_of_packets : 1.0e+5
iterations: 20
black_body_sampling:
start: 1 angstrom
stop: 1000000 angstrom
num: 1.e+6
last_no_of_packets: 1.e+5
no_of_virtual_packets: 5
convergence_strategy:
type: specific
damping_constant: 1.0
threshold: 0.05
fraction: 0.8
hold_iterations: 3
t_inner:
damping_constant: 1.0
spectrum:
start : 500 angstrom
stop : 20000 angstrom
num: 10000
Note
The specifications for the velocities of the inner and outer boundary values can be neglected (in which case TARDIS will default to using the full velocity range specified in the density.txt file). Values for the boundary velocities that lie outside the range covered by density.txt will not be accepted.
TARDIS can be used to compute a synthetic spectrum for a model with user-specified density and abundance profiles - this makes it possible to experiment with 1D profiles based on explosion models or with any empirical description of a model with stratified abundances or density.
The density profile is supplied in the form of a simple ascii file that should look something like this:
1 day
# index velocity (km/s) density (g/cm^3)
0 9000.0000 5.4869692e-10
1 10500.000 2.1759646e-10
2 12000.000 9.7656229e-11
3 13500.000 4.8170911e-11
4 15000.000 2.5600000e-11
5 16500.000 1.4450533e-11
6 18000.000 8.5733893e-12
7 19500.000 5.3037103e-12
8 21000.000 3.3999447e-12
9 22500.000 2.2474623e-12
In this file:
TARDIS will use this table of density versus velocity to specify the density distribution in the ejecta. For the calculation, TARDIS will use the reference time given in the file to scale the mass densities to whatever epoch is requested by assuming homologous expansion:
The values in the example here define a density profile that is dropping off with
Note
The grid of points specified in the input file is interpreted by TARDIS as defining a grid in which the tabulated velocities are taken as the outer boundaries of grid cells and the density is assumed to be uniform with each cell.
Warning
The example given here is to show the format only. It is not a realistic model. In any real calculation better resolution (i.e. more grid points) should be used.
For a model with density profile supplied via a file (see above), uniform abundances can be supplied as normal. Alternatively, a set of stratified elemental abundances can also be supplied. As with the density, the abundances are specified via an ascii file. An ascii file that could work with the example density file given above should be formatted like this:
# index Z=1 - Z=30
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1 0 0.1 0 0.2 0 0.2 0 0 0 0 0 0.2 0.1 0.1 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1 0 0.1 0 0.2 0 0.2 0 0 0 0 0 0.2 0.1 0.1 0 0
2 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1 0 0.1 0 0.2 0 0.2 0 0 0 0 0 0.2 0.1 0.1 0 0
3 0 0 0 0 0 0 0 0.2 0 0 0 0.1 0 0.2 0 0.2 0 0.2 0 0.1 0 0 0 0 0 0 0 0 0 0
4 0 0 0 0 0 0 0 0.2 0 0 0 0.1 0 0.2 0 0.2 0 0.2 0 0.1 0 0 0 0 0 0 0 0 0 0
5 0 0 0 0 0 0 0 0.2 0 0 0 0.1 0 0.2 0 0.2 0 0.2 0 0.1 0 0 0 0 0 0 0 0 0 0
6 0 0 0 0 0 0 0 0.2 0 0 0 0.1 0 0.2 0 0.2 0 0.2 0 0.1 0 0 0 0 0 0 0 0 0 0
7 0 0 0 0 0 0 0 0.2 0 0 0 0.1 0 0.2 0 0.2 0 0.2 0 0.1 0 0 0 0 0 0 0 0 0 0
8 0 0 0 0 0 0.5 0 0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
9 0 0 0 0 0 0.5 0 0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
In this file:
The abundances are specified as mass fractions (i.e. the sum of columns 1 to 30 in each row should be 1.0). TARDIS does not currently include any elements heavier that Z=30. The mass fractions specified will be adopted directly in the TARDIS calculations - so if your model is e.g. based on an explosion simulation you may need to calculate the state of any radioactive decay chains at the correct epoch.
The example file shown here has three simple layers:
Warning
The example given here is to show the format only. It is not a realistic model. In any real calculation better resolution (i.e. more grid points) should be used.
Warning
The calculation can be no better / more complete than the atomic data set. For further information on the atomic database - including details of how to develop your own dataset to suit your needs, please contact us.
If you create a correctly formatted density profile file (called “density.dat”) and abundance profile file (called “abund.dat”), you can use them in TARDIS by putting the following lines in the model section of the yaml file (and remove all other lines from these sections):
#New configuration for TARDIS based on YAML
#IMPORTANT any pure floats need to have a +/- after the e e.g. 2e+5
#Hopefully pyyaml will fix this soon.
---
#Currently only simple1d is allowed
tardis_config_version: v1.0
supernova:
luminosity_requested: 9.44 log_lsun
time_explosion: 13 day
atom_data: kurucz_atom_pure_simple.h5
model:
structure:
type: file
filename: density.dat
filetype: simple_ascii
v_inner_boundary: 11000 km/s
v_outer_boundary: 20000 km/s
abundances:
type: file
filename: abund.dat
filetype: simple_ascii
plasma:
disable_electron_scattering: no
ionization: lte
excitation: lte
radiative_rates_type: dilute-blackbody
line_interaction_type: macroatom
montecarlo:
seed: 23111963
no_of_packets : 1.0e+5
iterations: 20
black_body_sampling:
start: 1 angstrom
stop: 1000000 angstrom
num: 1.e+6
last_no_of_packets: 1.e+5
no_of_virtual_packets: 5
convergence_strategy:
type: specific
damping_constant: 1.0
threshold: 0.05
fraction: 0.8
hold_iterations: 3
t_inner:
damping_constant: 1.0
spectrum:
start : 500 angstrom
stop : 20000 angstrom
num: 10000
Note
The specifications for the velocities of the inner and outer boundary values can be neglected (in which case TARDIS will default to using the full velocity range specified in the density.txt file). Values for the boundary velocities that lie outside the range covered by density.txt will not be accepted.
Note
Increasing the number of virtual packets will improve the signal-to-noise of TARDIS spectra. You may wish to consider using a filter (e.g. Savitzky–Golay) to suppress the Monte Carlo noise for some applications.
Warning
The usefulness of any TARDIS calculations depends on the quality of the input atomic data. For further information on the atomic data - including details of how to develop your own dataset to suit your needs - please contact us.
Testing a code like TARDIS is an important step to making sure that the output is trustworthy. The TARDIS team has opted to use the py.test testing framework for use with TARDIS (and some astropy extensions of it). To run the tests download the desired version of TARDIS and run it with:
python setup.py test
This will run the currently implemented tests (which are not as many as there should be)
To quickly test TARDIS with any atomic dataset (it will only see if it in general runs):
python setup.py test --args="--atomic-dataset=<path_to_dataset> -s"
For a detailed description of the datasets go to Atomic Data Format Description
| file name | uuid1 | data sources | macroatom | zeta | synpp references | database version |
|---|---|---|---|---|---|---|
| kurucz_cd23_chianti_H_He.zip | 5ca3035ca8b311e3bb684437e69d75d7 | kurucz: Li I-II, Be I-IV, C I-IV, N I-VI O I-VI, F I-VI, Ne I-VI, Na I-VI Mg I-VI, Al I-VI, Si I-VI P I-VI, S I-VI, Cl I-V, Ar I-V K I-V, Ca I-IX, Sc I-IX, Ti I-IX V I-IX, Cr I-IX, Mn I-IX, Fe I-IX Co I-IX, Ni I-IX, Cu I-II Zn I-III chianti: H I,He I-II tardis_artificial_missing_ion: Li III, B V, C V-VI, N VII, O VII-VIII F VII-IX, Ne VII-X, Na VII-XI Mg VII-XII, Al VII-XIII Si VII-XIV, P VII-XV, S VII-XVI Cl VI-XVII, Ar VI-XVIII K VI-XIX, Ca X-XX, Sc X-XXI Ti X-XXII, V X-XXIII, Cr X-XXIV Mn X-XXV, Fe X-XXVI, Co X-XXVII Ni X-XXVIII, Cu III-XXIX Zn IV-XXX |
True | True | True | v0.9 |
Many of the Development workflow is taken from Astropy and credit belongs to the Astropy team for designing it.
The first step is to setup up a python environment. Python Environment for TARDIS explains how to do that.
In TARDIS we aim to stick to a test driven development. This uses the testing framework extensively starting with a test that shows this feature lacking via the implementation of the feature until the merging of the code to the main repository.
In most cases we try to break down big features into small, quantifiable goals which are then acted upon.
- Document feature to be added in an issue and maybe ask the mailing list if this feature exists
- Write a test that demonstrates what feature will be added.
- Run the test to verify that it fails in the way you think it should.
- If it fails in an unexpected way, your test may be wrong. This is a great time to ask the group for guidance
- If it passes, you are done! You just added test coverage to an already existing feature, and that is great! (unlikely)
- Add the feature (also known as “a simple matter of programming”).
- Run the test to verify that it passes.
- Write documentation about your feature.
- close issue/partial PR and add to changelog.
This section describes some of the Physics going on in TARDIS and will be gradually expanded
In order to compute the synthetic spectrum, Tardis must either be told or must calculate many physical properties of the model. To understand and test the code it can be important to look at these values. One easy way to do this is to run Tardis in an interactive mode and then inspect the model properties.
With iPython installed launch a session using
ipython --pylab
and then run your calculation, which is based on your “my_config.yml” file
from tardis import run_tardis
model = run_tardis('myconfig.yml')
If all goes well, the simulation should run as usual. Afterwards, the information from the simulation will all exist in “radial1d” and can be examined. Some examples for useful/interesting quantities are given below (but much more information is available: contact us via tardis-sn-users if you need further help).
For example, two of our important quantities are the parameters of the radiation field model, \(T_{\rm rad}\) and \(W\). These exist as Astropy Quantities. Thus
model.t_rads.cgs
will give you a list of the \(T_{\rm rad}\)-values for the model zones in cgs units. To obtain an array of the values (without units) use
model.t_rads.cgs.value
Similarly, the \(W\)-values can be accessed using
model.ws.cgs
Several important quantities that were setup when the model was defined by the configuration file are located in the “tardis_config” section. For example, the inner and outer velocity boundaries of the zones in the model
model.tardis_config.structure.v_inner.cgs
model.tardis_config.structure.v_outer.cgs
and the average density in the zones
model.tardis_config.structure.mean_densities.cgs
Many other interesting quantities are stored in the “plasma_array”. For example the calculated ion populations or level populations:
model.plasma_array.ion_populations
model.plasma_array.level_populations
These are stored as Pandas DataFrames. An index can be supplied to obtain the population in a particular zone. E.g., for the ion populations of the innermost zone (index = 0)
model.plasma_array.ion_populations[0]
Ion populations for a particular ionization stage of a particular element can be accessed by specifying an appropriate tuple \((Z,C)\), which identifies the element (via atomic number \(Z\) ) and the charge (via the ion charge \(C\) ). Thus,
model.plasma_array.ion_populations.ix[(14,1)]
will identify the ion popuations for Si II (\(Z=14, C=1\)) in all the zones. The above examples can be combined to obtain e.g. the Si II population in the innermost zone
model.plasma_array.ion_populations[0].ix[(14,1)]
The level populations are stored (and can be accessed) in a similar way - a third label can be used to pick out a particular atomic level. E.g., to pull out the population of the ground state (index 0) of Si II
model.plasma_array.level_populations.ix[(14,1,0)]
Note
If you prefer to work in SI units, all the astropy Quantities may instead by accessed with “xxx.si”.
Note
Information that is not stored as astropy Quantities (e.g. the ion an level populations used in the example above) are usually stored in cgs units (i.e. \({\rm cm}^{-3}\) for the populations).
We assume the supernova is in homologous expansion and we assume that there is photon injection from an inner boundary and we assume an outer boundary where these photons can leave the simulation again.
Montecarlo packets have a frequency, angle (\(\mu=\cos{\theta}\) ) and an energy: \(P(\nu, \mu, \epsilon)\). A large number \(n\) are generated at the start of each Montecarlo run in TARDIS by drawing \(\nu\) from a black-body distribution distribution with a \(\mu\) being drawn from a \(\sqrt[2]{z}\) distribution and an energy that is \(1/n\).
These packets are then launched into the simulation and there are two possible outcomes for each packet.
- packet leaves the simulation through the outer bound and count towards the
spectrum
- packet leaves the simulation through the inner boundary and are lost
(reabsorbed)
The packets can gain and loose energy throughout the simulation. If these packets scatter (through the various mechanisms), we transform into the comoving frame at that position \(v(r) = r / t_\textrm{explosion}\), where \(t_explosion\) is the time since explosion with the doppler factor \(d_\textrm{rest \rightarrow comoving}\) then change direction from \(\mu_\textrm{in}\) to \(mu_\textrm{out}\) and transform back and change the energy accordingly:
This way the montecarlo packets can gain or loose energy in the simulation:
(Source code, png, hires.png, pdf)
The spectrum is then generated as a weighted histogram. For each bin with edges \(\nu_{n}, \nu_{n+1}\) we get all the packets that left the simulation through the outer boundary and add up their remaining energies.
The calculations for the distance to the outer boundary:
The calculations for the distance to the inner boundary:
During the monte-carlo run we collect two estimators for the radiation field:
where \(\epsilon, \nu\) are comoving energy and comoving frequency of a packet respectively.
To calculate the temperature and dilution factor we first calculate the mean intensity in each cell ( \(J = \frac{1}{4\pi\, \Delta t\, V} J_\textrm{estimator}\) )., [Lucy03].
The weighted mean frequency is used to obtain the radiation temperature. Specifically, the radiation temperature is chosen as the temperature of a black body that has the same weighted mean frequency as has been computed in the simulation. Accordingly,
where the evaluation comes from the mean value of
and so
With the radiation temperature known, we can then obtain our estimate for for the dilution factor. Our radiation field model in the nebular approximation is
i.e. a dilute blackbody. Therefore we use our value of the mean intensity derrived from the estimator (above) to obtain the dilution factor
There endeth the lesson.
![digraph g{
a -> b -> c
c -> d [label="d_inner or \nd_outer"]
c -> e [label="d_line"]
d -> f [label="yes"]
d -> g [label="no"]
g -> a
e -> a [label="no"]
e -> h [label="yes"]
h -> a
a [label="We have a packet.",shape=box,fillcolor="white",style="filled,rounded"];
b [label="Calculate\nd_line, d_electron,\nd_inner and d_outer.",shape=box,fillcolor="white",style="filled,rounded"];
c [label="Which distance\nis smallest?", shape="diamond", fillcolor="white", style="filled"]
d [label="Are we leaving\nsimulation area?", shape="diamond", fillcolor="white", style="filled"]
e [label="Does the\npacket interact?", shape="diamond", fillcolor="white", style="filled"]
f [label="Packet is re-absorbed\nor emitted.\nThis ends the loop.", shape="box", fillcolor="white", style="filled,rounded"]
g [label="Update line\nprobabilities.", shape="box", fillcolor="white", style="filled,rounded"]
h [label="New random direction,\nupdated energy,\nmoving packet to current position,\nupdating event random number.", shape="box", fillcolor="white", style="filled,rounded"]
}](_images/graphviz-18fc53bc7bdea0590006e3769e4129f53732a0cf.png)
The LTEPlasma plasma class is the child of BasePlasma but is the first class that actually calculates plasma conditions. After running exactley through the same steps as BasePlasma, LTEPlasma will start calculating the partition functions.
, where Z is the partition function, g is the degeneracy factor, E the energy of the level and T the temperature of the radiation field.
The next step is to calculate the ionization balance using the Saha ionization equation. and then calculating the Number density of the ions (and an electron number density) in a second step. First \(g_e=\left(\frac{2 \pi m_e k_\textrm{B}T_\textrm{rad}}{h^2}\right)^{3/2}\) is calculated (in LTEPlasma.update_t_rad), followed by calculating the ion fractions (LTEPlasma.calculate_saha).
In the second step (LTEPlasma.calculate_ionization_balance), we calculate in an iterative process the electron density and the number density for each ion species.
Initially, we set the electron density (\(N_e\)) to the sum of all atom number densities. After having calculated the ion species number densities we recalculated the electron density by weighting the ion species number densities with their ion number (e.g. neutral ion number densities don’t contribute at all to the electron number density, once ionized contribute with a factor of 1, twice ionized contribute with a factor of two, ....).
Finally we calculate the level populations (LTEPlasma.calculate_level_populations), by using the calculated ion species number densities:
This concludes the calculation of the plasma. In the code, the next step is calculating the \(\tau_\textrm{Sobolev}\) using the quantities calculated here.
import os
from matplotlib import pyplot as plt
from matplotlib import colors
from tardis import atomic, plasma_array, util
import numpy as np
import pandas as pd
from astropy import units as u
#Making 2 Figures for ionization balance and level populations
plt.figure(1).clf()
ax1 = plt.figure(1).add_subplot(111)
plt.figure(2).clf()
ax2 = plt.figure(2).add_subplot(111)
# expanding the tilde to the users directory
atom_fname = os.path.join(os.path.dirname(atomic.__file__), 'data', 'atom_data.h5')
# reading in the HDF5 File
atom_data = atomic.AtomData.from_hdf5(atom_fname)
#The atom_data needs to be prepared to create indices. The Class needs to know which atomic numbers are needed for the
#calculation and what line interaction is needed (for "downbranch" and "macroatom" the code creates special tables)
atom_data.prepare_atom_data([14], 'scatter')
#Initializing the NebularPlasma class using the from_abundance class method.
#This classmethod is normally only needed to test individual plasma classes
#Usually the plasma class just gets the number densities from the model class
lte_plasma = plasma_array.BasePlasmaArray.from_abundance({'Si':1.0}, 1e-14*u.g/u.cm**3, atom_data, 10*u.day)
lte_plasma.update_radiationfield([10000], [1.0])
#Initializing a dataframe to store the ion populations and level populations for the different temperatures
ion_number_densities = pd.DataFrame(index=lte_plasma.ion_populations.index)
level_populations = pd.DataFrame(index=lte_plasma.level_populations.ix[14, 1].index)
t_rads = np.linspace(2000, 20000, 100)
#Calculating the different ion populations and level populuatios for the given temperatures
for t_rad in t_rads:
lte_plasma.update_radiationfield([t_rad], ws=[1.0])
#getting total si number density
si_number_density = lte_plasma.number_densities.get_value(14, 0)
#Normalizing the ion populations
ion_density = lte_plasma.ion_populations / si_number_density
ion_number_densities[t_rad] = ion_density
#normalizing the level_populations for Si II
current_level_population = lte_plasma.level_populations[0].ix[14, 1] / lte_plasma.ion_populations.get_value((14, 1), 0)
#normalizing with statistical weight
current_level_population /= atom_data.levels.ix[14, 1].g
level_populations[t_rad] = current_level_population
ion_colors = ['b', 'g', 'r', 'k']
for ion_number in [0, 1, 2, 3]:
current_ion_density = ion_number_densities.ix[14, ion_number]
ax1.plot(current_ion_density.index, current_ion_density.values, '%s-' % ion_colors[ion_number],
label='Si %s W=1.0' % util.int_to_roman(ion_number + 1).upper())
#only plotting every 5th radiation temperature
t_rad_normalizer = colors.Normalize(vmin=2000, vmax=20000)
t_rad_color_map = plt.cm.ScalarMappable(norm=t_rad_normalizer, cmap=plt.cm.jet)
for t_rad in t_rads[::5]:
ax2.plot(level_populations[t_rad].index, level_populations[t_rad].values, color=t_rad_color_map.to_rgba(t_rad))
ax2.semilogy()
t_rad_color_map.set_array(t_rads)
cb = plt.figure(2).colorbar(t_rad_color_map)
ax1.set_xlabel('T [K]')
ax1.set_ylabel('Number Density Fraction')
ax1.legend()
ax2.set_xlabel('Level Number for Si II')
ax2.set_ylabel('Number Density Fraction')
cb.set_label('T [K]')
plt.show()
The NebularPlasma class is a more complex description of the Plasma state than the LTEPlasma. It takes a dilution factor (W) into account, which deals with the dilution of the radiation field due to geometric, line-blocking and other effects.
The calculations follow the same steps as LTEPlasma, however the calculations are different and often take into account if a particular level is meta-stable or not. NebularPlasma will start calculating the partition functions.
, where Z is the partition function, g is the degeneracy factor, E the energy of the level, T the temperature of the radiation field and W the dilution factor.
The next step is to calculate the ionization balance using the Saha ionization equation. and then calculating the Number density of the ions (and an electron number density) in a second step. In the first step, we calculate the ionization balance using the LTE approximation (\(\Phi_{i, j}(\textrm{LTE})\)). Then we adjust the ionization balance using two factors \(\zeta\) and \(\delta\).
\(\zeta\) is read in for specific temperatures and then interpolated for the target temperature.
\(\delta\) is a radiation field correction factors which is calculated according to Mazzali & Lucy 1993 ([MazzaliLucy93]; henceforth ML93)
In ML93 the radiation field correction factor is denoted as \(\delta\) and is calculated in Formula 15 & 20
The radiation correction factor changes according to a ionization energy threshold \(\chi_\textrm{T}\) and the species ionization threshold (from the ground state) \(\chi_0\).
For \(\chi_\textrm{T} \ge \chi_0\)
For \(\chi_\textrm{T} < \chi_0\)
where \(T_\textrm{R}\) is the radiation field Temperature, \(T_\textrm{e}\) is the electron temperature and W is the dilution factor.
Now we can calculate the ionization balance using equation 14 in [MazzaliLucy93]:
In the last step, we calculate the ion number densities according using the methods in LTEPlasma
Finally we calculate the level populations (NebularPlasma.calculate_level_populations()),
by using the calculated ion species number densities:
This concludes the calculation of the nebular plasma. In the code, the next step is calculating the \(\tau_\textrm{Sobolev}\) using the quantities calculated here.
import os
from matplotlib import colors
from tardis import atomic, plasma_array, util
from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
#Making 2 Figures for ionization balance and level populations
plt.figure(1).clf()
ax1 = plt.figure(1).add_subplot(111)
plt.figure(2).clf()
ax2 = plt.figure(2).add_subplot(111)
# expanding the tilde to the users directory
atom_fname = os.path.expanduser('~/.tardis/si_kurucz.h5')
# reading in the HDF5 File
atom_data = atomic.AtomData.from_hdf5(atom_fname)
#The atom_data needs to be prepared to create indices. The Class needs to know which atomic numbers are needed for the
#calculation and what line interaction is needed (for "downbranch" and "macroatom" the code creates special tables)
atom_data.prepare_atom_data([14], 'scatter')
#Initializing the NebularPlasma class using the from_abundance class method.
#This classmethod is normally only needed to test individual plasma classes
#Usually the plasma class just gets the number densities from the model class
nebular_plasma = plasma.NebularPlasma.from_abundance(10000, 0.5, {'Si': 1}, 1e-13, atom_data, 10.)
#Initializing a dataframe to store the ion populations and level populations for the different temperatures
ion_number_densities = pd.DataFrame(index=nebular_plasma.ion_populations.index)
level_populations = pd.DataFrame(index=nebular_plasma.level_populations.ix[14, 1].index)
t_rads = np.linspace(2000, 20000, 100)
#Calculating the different ion populations and level populuatios for the given temperatures
for t_rad in t_rads:
nebular_plasma.update_radiationfield(t_rad, w=1.0)
#getting total si number density
si_number_density = nebular_plasma.number_density.get_value(14)
#Normalizing the ion populations
ion_density = nebular_plasma.ion_populations / si_number_density
ion_number_densities[t_rad] = ion_density
#normalizing the level_populations for Si II
current_level_population = nebular_plasma.level_populations.ix[14, 1] / nebular_plasma.ion_populations.ix[14, 1]
#normalizing with statistical weight
current_level_population /= atom_data.levels.ix[14, 1].g
level_populations[t_rad] = current_level_population
ion_colors = ['b', 'g', 'r', 'k']
for ion_number in [0, 1, 2, 3]:
current_ion_density = ion_number_densities.ix[14, ion_number]
ax1.plot(current_ion_density.index, current_ion_density.values, '%s-' % ion_colors[ion_number],
label='Si %s W=1.0' % util.int_to_roman(ion_number + 1).upper())
#only plotting every 5th radiation temperature
t_rad_normalizer = colors.Normalize(vmin=2000, vmax=20000)
t_rad_color_map = plt.cm.ScalarMappable(norm=t_rad_normalizer, cmap=plt.cm.jet)
for t_rad in t_rads[::5]:
ax2.plot(level_populations[t_rad].index, level_populations[t_rad].values, color=t_rad_color_map.to_rgba(t_rad))
ax2.semilogy()
#Calculating the different ion populations for the given temperatures with W=0.5
ion_number_densities = pd.DataFrame(index=nebular_plasma.ion_populations.index)
for t_rad in t_rads:
nebular_plasma.update_radiationfield(t_rad, w=0.5)
#getting total si number density
si_number_density = nebular_plasma.number_density.get_value(14)
#Normalizing the ion populations
ion_density = nebular_plasma.ion_populations / si_number_density
ion_number_densities[t_rad] = ion_density
#normalizing the level_populations for Si II
current_level_population = nebular_plasma.level_populations.ix[14, 1] / nebular_plasma.ion_populations.ix[14, 1]
#normalizing with statistical weight
current_level_population /= atom_data.levels.ix[14, 1].g
level_populations[t_rad] = current_level_population
#Plotting the ion fractions
for ion_number in [0, 1, 2, 3]:
print "w=0.5"
current_ion_density = ion_number_densities.ix[14, ion_number]
ax1.plot(current_ion_density.index, current_ion_density.values, '%s--' % ion_colors[ion_number],
label='Si %s W=0.5' % util.int_to_roman(ion_number + 1).upper())
for t_rad in t_rads[::5]:
ax2.plot(level_populations[t_rad].index, level_populations[t_rad].values, color=t_rad_color_map.to_rgba(t_rad),
linestyle='--')
ax2.semilogy()
t_rad_color_map.set_array(t_rads)
cb = plt.figure(2).colorbar(t_rad_color_map)
ax1.set_xlabel('T [K]')
ax1.set_ylabel('Number Density Fraction')
ax1.legend()
ax2.set_xlabel('Level Number for Si II')
ax2.set_ylabel('Number Density Fraction')
cb.set_label('T [K]')
plt.show()
The macro atom is described in detail in [Lucy02]. The basic principal is that when an energy packet is absorbed that the macro atom is on a certain level. Three probabilities govern the next step the Pup, Pdown and Pdown emission being the probability for going to a higher level, a lower level and a lower level and emitting a photon while doing this respectively (see Figure 1 in [Lucy02] ).
The macro atom is the most complex idea to implement as a data structure. The setup is done in ~tardisatomic, but we will nonetheless discuss it here (as ~tardisatomic is even less documented than this one).
For each level we look at the line list to see what transitions (upwards or downwards are possible). We create a two arrays, the first is a long one-dimensional array containing the probabilities. Each level contains a set of probabilities, The first part of each set contains the upwards probabilities (internal upward), the second part the downwards probabilities (internal downward), and the last part is the downward and emission probability.
each set is stacked after the other one to make one long one dimensional ~numpy.ndarray.
The second array is for book-keeping it has exactly the length as levels (with an example for the Si II level 15):
| Level ID | Probability index | Nup | Ndown | Ntotal |
|---|---|---|---|---|
| 14001015 | ??? | 17 | 5 | 17 + 5*2 = 27 |
We now will calculate the transition probabilites, using the radiative rates in Equation 20, 21, and 22 in [Lucy02]. Then we calculate the downward emission probability from Equation 5, the downward and upward internal transition probabilities in [Lucy03].
We ignore the probability to emit a k-packet as TARDIS only works with photon packets. Next we calculate the radidative rates using equation 10 in [Lucy03].
with \(\beta_\textrm{Sobolev} = \frac{1}{\tau_\textrm{Sobolev}}(1-e^{-\tau_\textrm{Sobolev}})\) .
using the Einstein coefficients
we get
This results in the transition probabilities:
There are two parts for each of the probabilities, one that is pre-computed by ~tardisatomic and is in the HDF5 File, and one that is computed during the plasma calculations:
NLTE treatment of lines is available both in ~LTEPlasma and the ~NebularPlasma class. This can be enabled by specifying which species should be treated as NLTE with a simple list of tuples (e.g. [(20,1)] for Ca II).
First let’s dive into the basics:
There are two rates to consider from a given level.
where \(\bar{J}_\nu\) (in LTE this is \(B(\nu, T)\)) denotes the mean intensity at the frequency of the line and \(n_e\) the number density of electrons.
Next, we calculate the rate of change of a level by adding up all outgoing and all incoming transitions from level \(j\).
In a statistical equilibrium all incoming rates and outgoing rates add up to 0 (\(\frac{dn_j}{dt}=0\)). We use this to calculate the level populations using the rate coefficients (\(r_ij, r_ji\)).
with the additional constrained that all the level number populations need to add up to the current ion population $N$ we change this to
For a three level atom we have:
which can be written in matrix from:
To solve for the level populations we need an additional constraint: \(n_1 + n_2 + n_3 = N\). By setting \(N = 1\): we can get the relative rates:
Now we go back and look at the rate coefficients used for a level population - as an example \(\frac{dn_2}{dt}\):
Next we will group the stimulated emission and stimulated absorption terms as we can assume \(\bar{J_{12}} = \bar{J_{21}}\):
Paper version:
This section will summarise the equations used in the calculation of the helium state for the recomb-NLTE approximation in Tardis. A full physical justification for these equations will be provided in an upcoming paper. All of the level populations are given as a function of the He II ground state population (\(n_{2,1,0}\)), and the values are then normalised using the helium number density to give the correct level number densities.
- i: atomic number
- j: ion number
- k: level number
(Note: An extra factor of \(\frac{1}{W}\) is included for the metastable states of He I.)
Code Version:
In the Tardis plasma, some of these equations are re-written slightly to make use of existing property methods (e.g. PhiSahaLTE, PhiSahaNebular) often using the relation:
Another helium_treatment option offered by Tardis is numerical-nlte. The use of this plasma property requires an additional code that is the property of Stephan Hachinger (see arXiv:1201.1506) and is not distributed with Tardis. Tardis also requires a specific atomic datafile to use this module. This plasma option is included so that people who have access to and permission to use the necessary module may use it. Otherwise, the recomb-NLTE option provides a very accurate alternative approximation.
with a 4x speed increase.
| [Lucy02] | L. B. Lucy. Monte Carlo transition probabilities. \aap , 384:725–735, March 2002. arXiv:arXiv:astro-ph/0107377, doi:10.1051/0004-6361:20011756. |
| [Lucy03] | L. B. Lucy. Monte Carlo transition probabilities. II. \aap , 403:261–275, May 2003. arXiv:arXiv:astro-ph/0303202, doi:10.1051/0004-6361:20030357. |
| [MazzaliLucy93] | P. A. Mazzali and L. B. Lucy. The application of Monte Carlo methods to the synthesis of early-time supernovae spectra. \aap , 279:447–456, November 1993. |
- Wolfgang Kerzendorf
- Stuart Sim
- Michael Klauser
- Markus Kromer
- Maryam Patel (documentation, test environment)
- Adam Suban-Loewen (GUI, profiling and parallelization)
- Chris Sasaki and Mike Reid (Logo)
- Knox Long & Christian Knigge (comparisons with “Python” and ionization data)
- Erik Bray (help with incorporating Astropy’s setuphelpers)
- Robert Ryans (testing installation process)
TARDIS-SN is licensed under a 3-clause BSD style license:
Copyright (c) 2014, TARDIS-SN Developers All rights reserved.
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
Neither the name of the TARDIS-SN nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
The code is built on a few principles:
- open - the code is fully open source and we invite usage and contributions from the community
- modular - the code has different microphysics modules and can be easily extended
- fast - the code is geared towards rapid spectral synthesis to fit supernovae and other transients
- easy - the code is designed to be easily installed and run as well as a detailed documentation
We encourage you to subscribe to tardis-sn-user to ask questions about TARDIS.
If you use this code for any publications or presentations please acknowledge it accordingly. For this first version please mention the website and cite Kerzendorf & Sim 2014.
User modifications and additions that lead to publications need to be handed back to the community by incorporating them into this publicly available version of TARDIS.
The current stable version of TARDIS is 0.9.2 and can be downloaded here, further installation instructions are available here Installation.
A file containing an example configuration file and an atomic database can be found in the section Running TARDIS
If you’re interested in contributing to the code, either contact us or you can contribute directly via github. We are using Astropy’s excellent workflow - more details can be found at http://astropy.readthedocs.org/en/latest/development/workflow/maintainer_workflow.html.
Warning
Currently TARDIS only works on 64-bit python installations. We’re working on making it work on 32-bit python distributions.