First example code for the use of densities

All densities are separated in atomic hydrogen (HI), ionised hydrogen (HII) and molecular hydrogen (H2).

This example contains:

• constant densities
• superposition of different types
• superposition of the same types

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

from crpropa import *

First we start with the use of a constant density.

The first double sets the HI, the second the HII and the third the H2 density-number 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.

For this purpose methods are available to change and activate (set-function) and methodes to see the 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 active and has a density of 1e+06 cm^-3
HII component is not active and has a density of 1.5e+06 cm^-3
H2 component is active 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

To customize a density use the DensityList. In a superposition of global models of the density distribution of the Milky Way, the normalisation of the different components has to be considered. The relative normalisations are conserved when you just add components after deactivating the unwanted parts.

CD1 = ConstantDensity(0,2,0)     # for use of HII
CD2 = ConstantDensity(3,1,2.5)   # for use of 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