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



In [2]:

    
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



In [3]:

    
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}$)



In [4]:

    
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)



In [7]:

    
#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



In [8]:

    
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 globel models of density distribution of the Milkyway keep care of normalisation. Therfore you can just add components by deactivating the other



In [10]:

    
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



In [11]:

    
#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



In [12]:



In [ ]: