GETTING STARTED To get started, move into directory where PyNeb resides and enter python
In [1]:
%matplotlib inline
# import code and modules
import pyneb as pn
import numpy as np
import matplotlib.pyplot as plt
In [2]:
#######################################################################
# DEFINING ATOMS
# define an OII atom
O2 = pn.Atom("O", 2)
# alternate syntax to define an atom (spec is a string)
N2 = pn.Atom("N", "2")
In [3]:
# check atom definition
print O2.elem
print O2.spec
print O2.atom
print O2.name
In [4]:
# explore the atom: builtin data
print O2.gs # ground-state configuration
In [5]:
# array of stat weights
print O2.getStatWeight()
In [6]:
# stat weight of a given level
lev_i = 2
print O2.getStatWeight(lev_i)
In [7]:
# explore the atom: adopted atomic data
pn.atomicData.getPredefinedDataFileDict() # we suggest using the tab for this command...
pn.atomicData.getDirForFile('o_ii_atom_WFD96.fits') # wanna know where the file lies?
O2.printSources() # print bibliographic references of data used to build O2
In [8]:
print O2.NLevels # number of levels in the selected data
print O2.getEnergy(2) # energy of first excited level (ground = 1) in Angstrom^-1
print O2.getA(2,1) # transition probability of 2->1
In [9]:
N2.plotGrotrian()
In [10]:
# set temperature and density
tem = 15000.
den = 1000.
print O2.getPopulations(tem, den) # compute populations
In [11]:
print O2.getCritDensity(tem, level=2) # critical density of level 2 at tem
In [12]:
print O2.getOmega(tem, 2, 1) # effective collision strength of transition 2->1 at T=10000K
In [13]:
print O2.getOmegaArray(2, 1) # array of effective collision strengths for 2->1 as a function of T
In [14]:
print O2.getTemArray() # print array of temperatures of tabulated Omegas
In [15]:
print O2.getCollRates(tem) # print collisional Rates at T=tem
In [16]:
pn.atomicData.getAllAvailableFiles('O3')
Out[16]:
In [17]:
# This bit calls the script DataPlot.py to plot atomic data.
dataplot = pn.DataPlot('O', 3)
dataplot.plotAllA(figsize=(14, 10)) # transition probabilities plot
In [18]:
pn.config.setNlevels(6, 'O3', 'coll')
dataplot.plotOmega() # collision strength plot
In [19]:
# customize atomic data
# First step: check which directories are searched for atomic data files
pn.atomicData.getAllDataFilePaths()
Out[19]:
In [20]:
# Add your selected directory to the list
pn.atomicData.addDataFilePath('/tmp')
In [21]:
# Check if it's been added
pn.atomicData.getAllDataFilePaths()
Out[21]:
In [22]:
# Remove it if you gave the wrong dir
pn.atomicData.removeDataFilePath('/tmp')
In [23]:
# Set 'o_iii_.fits' to be the OIII atom file
pn.atomicData.setDataFile('o_iii_coll_Pal12-AK99.dat')
In [24]:
# define an atom with the new data
O2test = pn.Atom("O", 2)
In [25]:
# define all atoms at once and put them in a dictionary
# (all of them defined with the latest dataset selected)
atoms = pn.getAtomDict() #this generates a lot of warnings as not all element-spectrum combination exist
In [26]:
# see what atoms have been built
atoms
Out[26]:
In [27]:
# build only some atoms
atoms = pn.getAtomDict(atom_list=['O1', 'O2', 'O3', 'N2', 'N3'])
In [28]:
# explore some specific atom in the atoms collection
atoms['N2'].elem
Out[28]:
In [29]:
# if you want to be able to access them directly rather than through a dictionary:
for key in atoms.keys():
vars()[key]=atoms[key]
In [30]:
# for example
O2.elem
Out[30]:
In [31]:
# list all atom features
dir(O2)
Out[31]:
In [32]:
#######################################################################
# MAKING CALCULATIONS
# set temperature and density
tem = 15000.
den = 1000.
# compute populations
O2.getPopulations(tem, den)
Out[32]:
In [33]:
# compute emissivity of transition (lev_i, lev_j)
O2.getEmissivity(tem, den, 3, 2)
Out[33]:
In [34]:
# also works if tem is an array
tem = np.array([10000, 20000, 30000])
print O2.getOmega(tem, 2, 1)
print O2.getCollRates(tem)
print O2.getPopulations(tem, den)
In [35]:
# tem and den can be arrays as well as single numbers
tem = np.array([10000, 12000, 13000])
den = np.array([100, 200, 300])
print O2.getPopulations(tem, den) # returns the n_tem x n_den x n_levels array of populations
print O2.getPopulations(tem, den, product=False) # element-by-element multiplication of tem and den (no scalar product: returns [pop(tem_1, den_1), pop(tem_2, den_2), ... pop(tem_n, den_n)]
In [36]:
# find transition corresponding to given wavelength
N2.printTransition(6584)
In [37]:
# temperature determination from an intensity ratio
N2.getTemDen(0.01, den=1000., wave1=5755, wave2=6584)
Out[37]:
In [38]:
# same as above, by specifying the levels involved
N2.getTemDen(0.01, den=1000., lev_i1=5, lev_j1=4, lev_i2=4, lev_j2=3)
Out[38]:
In [39]:
# same as above, by specifying the levels involved
N2.getTemDen(0.01, den=1000., to_eval = 'L(5755) / L(6584)')
Out[39]:
In [40]:
# no formal difference between temperature and density diagnostics, so beware of what you do
print N2.getTemDen(0.01, tem=8782., wave1=5755, wave2=6584)
print N2.getTemDen(0.01, tem=8882., wave1=5755, wave2=6584)
In [41]:
# ionic abundance (intensity, temperature, density, transition)
O2.getIonAbundance(100, 1.5e4, 100., wave=3727)
Out[41]:
In [42]:
# printout as in old nebular
O2.printIonic() # only prints transitions and corresponding wavelengths. Useful for a quick glance at the atom.
In [43]:
O2.printIonic(tem=10000, den=100) # also prints line emissivities
In [44]:
O2.printIonic(tem=10000, den=100, printA=True, printPop=True, printCrit=True) # also prints populations and critical densities
In [45]:
# Compute Hb emissivity at T=10000K
H1 = pn.RecAtom('H', 1)
H1.getEmissivity(10000, 1e2, 4, 2)
Out[45]:
In [46]:
# simultaneously compute temperature and density from pairs of line ratios
# First of all, a Diagnostics object must be created and initialized with the relevant diagnostics.
diags = pn.Diagnostics() # this creates the object
diags.getAllDiags() # see what Diagnostics exist
tem, den = diags.getCrossTemDen('[NII] 5755/6548', '[SII] 6731/6716', 0.050, 1.0)
In [47]:
#TO BE CONTINUED FROM HERE
print tem, den
In [47]: