first example code for use of densitys

all densitys are seperated in atomic hydrogyn (HI), ioninised hydrogyn (HII) and molecular hydrogn (H2)

this example contains

• constant densitys
• superposition of different typses
• superposition of 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 set 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 densitys * 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 useage of several parts can be chosen

therefore are methodes to change and activate (set-function) and methodes 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 acitvated 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 globel models of density distribution of the Milkyway keep care of normalisation. Therfore you can just add components by deactivating the other

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