First example code for use of densities

All densities are seperated in atomic hydrogen (HI), ioninised hydrogen (HII), and molecular hydrogen (H2).

This example contains * constant densities * superposition of different types * superposition of same types

Further on, there is another example for all implemented density models for the Milky Way.

from crpropa import *

First, we start with the use of a constant density. The first argument sets the HI, the second the HII, and the third the H2 density numbers, in SI units.

CD = ConstantDensity(1, 2, 0.5)

To see the output option we check the density at a random position.

The output options are: * getDensity: for the sum of all densities * getHIDensity: for the HI part * getHIIDensity: for the HII part * getH2Density: for the H2 part * getNucleonDensity: for the sum of nuclei (\(n_{HI} + n_{HII} + 2\cdot n_{H2}\))

position = Vector3d(2, 1, 3)

n_tot = CD.getDensity(position)
n_HI = CD.getHIDensity(position)
n_HII = CD.getHIIDensity(position)
n_H2 = CD.getH2Density(position)
n_nucl = CD.getNucleonDensity(position)


print('total density n = %f' % n_tot)
print('HI density n_HI = %f' % n_HI)
print('HII density n_HII = %f' % n_HII)
print('H2 density n_H2 = %f' % n_H2)
print('nucleon density n_nucl = %f' % n_nucl)

total density n = 3.500000
HI density n_HI = 1.000000
HII density n_HII = 2.000000
H2 density n_H2 = 0.500000
nucleon density n_nucl = 4.000000

The ConstantDensity can be adjusted to new values and the usage of several parts can be chosen.

There are methods to change and activate (setter functions) and methods to see actual configuration (getisfor-functions).

#see the actual configuration
print('HI: %s, HII: %s, H2: %s, \n \n' % (CD.getIsForHI(), CD.getIsForHII(), CD.getIsForH2()))

# change activity
CD.setHI(False)

# change activity and density number
CD.setHII(False, 1.5)

# change density number
CD.setH2(1.3)

# see the changes in the Description
print(CD.getDescription())

HI: True, HII: True, H2: True,


ConstantDensity:HI component is not activ and has a density of 1e+06 cm^-3      HII component is not activ and  has a density of 1.5e+06 cm^-3      H2 component is activ and  has a density of 1.3e+06 cm^-3

The output of the getDensity and getNucleonDensity depends on the activity of the types. Only activated types are used for summing up

n_tot = CD.getDensity(position)
n_HI = CD.getHIDensity(position)
n_HII = CD.getHIIDensity(position)
n_H2 = CD.getH2Density(position)
n_nucl = CD.getNucleonDensity(position)


print('total density n = %f' % n_tot)
print('HI density n_HI = %f' % n_HI)
print('HII density n_HII = %f' % n_HII)
print('H2 density n_H2 = %f' % n_H2)
print('nucleon density n_nucl = %f' % n_nucl)

total density n = 1.300000
HI density n_HI = 1.000000
HII density n_HII = 1.500000
H2 density n_H2 = 1.300000
nucleon density n_nucl = 2.600000

Customize a density with different density models

To customize a density, use the DensityList. In a superposition of global models of density distributions of the Milky Way, normalisation is automatically taken care of. Therefore, you can just add components by deactivating the others.

CD1 = ConstantDensity(0, 2, 0)     # for use HII
CD2 = ConstantDensity(3, 1, 2.5)   # for use HI, H2

CUS = DensityList()

# first deactivate not wanted parts

CD1.setHI(False)
CD1.setH2(False)
CD2.setHII(False)

# add density
CUS.addDensity(CD1)
CUS.addDensity(CD2)

# get output
n_tot = CUS.getDensity(position)
n_HI = CUS.getHIDensity(position)
n_HII = CUS.getHIIDensity(position)
n_H2 = CUS.getH2Density(position)
n_nucl = CUS.getNucleonDensity(position)

print('total density n = %f' % n_tot)
print('HI density n_HI = %f' % n_HI)
print('HII density n_HII = %f' % n_HII)
print('H2 density n_H2 = %f' % n_H2)
print('nucleon density n_nucl = %f' % n_nucl)

total density n = 7.500000
HI density n_HI = 3.000000
HII density n_HII = 3.000000
H2 density n_H2 = 2.500000
nucleon density n_nucl = 10.000000

DensityList

you can also superposition total models without deactivating several components

#set wanted density
CD1 = ConstantDensity(1, 3, 4)
CD2 = ConstantDensity(1.5, 2, 3.3)

DL = DensityList()

# add density to list

DL.addDensity(CD1)
DL.addDensity(CD2)

# see output
n_tot = DL.getDensity(position)
n_HI = DL.getHIDensity(position)
n_HII = DL.getHIIDensity(position)
n_H2 = DL.getH2Density(position)
n_nucl = DL.getNucleonDensity(position)

print('total density n = %f' % n_tot)
print('HI density n_HI = %f' % n_HI)
print('HII density n_HII = %f' % n_HII)
print('H2 density n_H2 = %f' % n_H2)
print('nucleon density n_nucl = %f' % n_nucl)

total density n = 14.800000
HI density n_HI = 2.500000
HII density n_HII = 5.000000
H2 density n_H2 = 7.300000
nucleon density n_nucl = 22.100000