In [1]:
import numpy as np
import matplotlib.pyplot as plt
import pyneb as pn
This is the class characterizing emission lines. An emission line is identified by an element and a spectrum (which identify the emitting ion), a wavelength in Angstrom, a blend flag, a label in the standard PyNeb format, an observed intensity, a reddening-corrected intensity, an expression describing how the intensity depends on the included wavelengths, an observational error and an error on the corrected intensity. Other programs determine one or more of these values.
To instantiate an Emission Line object, use the following:
In [2]:
line = pn.EmissionLine('O', 3, 5007, obsIntens=[1.4, 1.3])
line2 = pn.EmissionLine(label = 'O3_5007A', obsIntens=320, corrected=True)
obsIntens is a value, a list or a numpy array of values corresponding to the observed intensity of the given emission line.
In [3]:
print(line)
To know how the label of a given line is exactly spelled, you can print the dictionary pn.LINE_LABEL_LIST
In [4]:
print(pn.LINE_LABEL_LIST['O3'])
It is possible to instantiate a line not contained in the pn.LINE_LABEL_LIST. In this case a warning is issued, but the code doesn't stop.
The observed intensity is stored in line.obsIntens and the extinction-corrected intensity is stored in line.corrIntens. line.corrIntens is set to 0.0 when the line is instantiated, unless the parameter corrected is set to True, in which case the observed value obsIntens is copied to the corrIntens slot (the same applies for corrError, which is set to obsError).
The corrIntens value can also be computed using an instantiation of the pn.RedCorr class:
In [5]:
redcorr = pn.RedCorr(E_BV = 0.87, law = 'F99')
In [6]:
line.correctIntens(redcorr) #redcorr is used to compute line.corrIntens
The line information is printed using:
In [7]:
line.printLine()
Most of the times, users will not need to define or manipulate EmissionLine objects, since most of the work on the EmissionLine objects will be performed from the Observation class (read data, extinction correction); see next section.
WARNING: Note that the wavelengths assigned to EmissionLine objects are simply the numerical part of the label:
In [8]:
Hb1 = pn.EmissionLine(label='H1r_4861A').wave
print(Hb1)
whereas pn.Atom computes them as the difference from energy levels:
In [9]:
Hb2 = pn.RecAtom('H', 1).getWave(4, 2)
print(Hb2)
This can cause small errors when both methods are used simultaneously. For instance, the extinction correction at Hb1 is slightly different from the expected value of 1:
In [10]:
rc = pn.RedCorr()
rc.E_BV = 1.34
rc.law = 'F99'
rc.getCorrHb(Hb1)
Out[10]:
This happens because the ExtCorr uses the precise H$\beta$ value computed from energy levels. If this roundoff error exceeds your tolerance, a possible workaround is forcing the emission line to have exactly the wavelength computed from the energy levels:
In [11]:
O3w = pn.EmissionLine('O', 3, wave=5007).wave
print(rc.getCorrHb(O3w))
In [12]:
Hb11 = pn.EmissionLine('H', 1, wave=Hb2).wave
print(rc.getCorrHb(Hb11))
# This will generate a warning (as the transition is not included in the inventory for the specified atom),
# but the code won't stop.
pn.Observation is the class characterizing observation records. An observation record is composed of an object identifier, the observed intensity of one or more emission lines, and, optionally, the dereddened line intensities and the identifier of the extinction law used, and the value of c(Hbeta).
Observations can be initialized by reading data files, which can be organized with different emission lines either in rows or columns (usually, in a survey of many objects with few emission lines emission lines change across columns; and in a high-resolution observation of a small sample of objects lines change across rows).
The following is an example of how to define an observation:
In [13]:
obs = pn.Observation()
In [14]:
%%writefile observations1.dat
LINE SMC_24
S4_10.5m 7.00000
Ne2_12.8m 8.3000
Ne3_15.6m 34.10
S3_18.7m 10.
O2_3726A 39.700
O2_3729A 18.600
Ne3_3869A 18.90
Ne3_3968A 6.4
S2_4069A 0.85
S2_4076A 0.450
O3_4363A 4.36
H1r_4861A 100.00
O3_5007A 435.09
N2_5755A 0.510000
S3_6312A 0.76
O1_6300A 1.69
O1_6364A 0.54
N2_6548A 6.840000
H1r_6563A 345.00
N2_6584A 19.00
S2_6716A 1.220000
S2_6731A 2.180000
Ar3_7136A 4.91
O2_7319A+ 6.540000
O2_7330A+ 5.17
In [15]:
obs.readData('observations1.dat', fileFormat='lines_in_rows', err_default=0.05) # fill obs with data read from smc_24.dat
In [16]:
obs.printIntens(returnObs=True)
In [17]:
obs.extinction.law = 'CCM89' # define the extinction law from Cardelli et al.
obs.correctData() # the dereddened data are computed
The data can be read by the readData method as above or directly while instantiating the object:
In [18]:
obs = pn.Observation('observations1.dat', fileFormat='lines_in_rows', corrected=True)
The format of the data file from which the emission line intensities are read can be one of three kinds: "lines_in_rows" as above, or "lines_in_cols" like this one:
In [19]:
%%writefile observations2.dat
NAME O2_3726A O2_3726Ae O2_3729A O2_3729Ae
NGC3132 0.93000 0.05000 0.17224200 0.10
IC418 1.28000 0.05000 0.09920000 0.05
M33 0.03100 0.080 0.03100 0.10
In [20]:
obs2 = pn.Observation('observations2.dat', fileFormat='lines_in_cols', corrected=True)
or fileFormat='lines_in_rows_err_cols' (errors labeled “err”. Don’t name an observation “err”!) like this one:
In [21]:
%%writefile observations3.dat
LINE TT err TT2 err TT3 err
cHbeta 1.2 0.0 1.5 0.2 1.1 0.2
O3_5007A 1.5 0.15 1.3 .2 1.1 0.1
H1_6563A 2.89 0.05 1.6 0.3 1.3 0.1
N2_6584A 1. 0.20 0.3 0.5 1.5 0.1
In [22]:
#obs3 = pn.Observation('observations3.dat', fileFormat='lines_in_rows_err_cols', corrected=False)
The delimiter between the columns is any sequence of spaces or TAB, but it can be changed using the delimiter parameter. The line names are defined by a label starting with the name of the atom (‘O2’), followed by an underscore, followed by a wavelength and ending with a unit (‘A’ or ‘m’). The list of all the lines managed by PyNeb, ordered by atoms, is obtained by entering:
In [23]:
for atom in pn.LINE_LABEL_LIST:
print(atom, pn.LINE_LABEL_LIST[atom])
The presence of a trailing “e” at the end of the label points to the error associated to the line. The error is considered to be relative to the intensity (i.e., 0.05 means 5% of the intensity), unless the parameter errIsRelative is set to False. A common value for all the errors can be defined by the parameter err_default (0.10 is the default value).
Once the data have been read, they have to be corrected from extinction. An instantiation of RedCorr() is available inside the Observation object as obs.extinction.
If the data file contains cHbeta or E(B-V) alongside of line labels, the corresponding information on extinction is transmitted to the extinction correction object. Otherwise, the extinction parameters must be set manually; for example:
In [24]:
obs = pn.Observation('observations1.dat', fileFormat='lines_in_rows', corrected=True)
obs.extinction.cHbeta = 1.2
obs.extinction.E_BV = 0.34
An extinction law has to be specified in either case:
In [25]:
obs.extinction.law = 'F99'
To correct all the lines at once:
In [26]:
obs.correctData()
In [27]:
obs.printIntens(returnObs=True)
In [28]:
obs.printIntens()
If you want the corrected line intensities to be normalized to a given wavelength, use the following:
In [29]:
obs.correctData(normWave=4861.)
The extinction correction can be determined by comparing the observed values to a theoretical ratio, as in the following:
In [30]:
obs.printIntens()
In [31]:
obs.def_EBV(label1="H1r_6563A", label2="H1r_4861A", r_theo=2.85)
print(obs.extinction.E_BV)
obs.correctData(normWave=4861.)
In [32]:
obs.printIntens()
By default, this method prints out the corrected intensities. To print the observed intensities, use the returnObs=True parameter.
The method getSortedLines returns the lines sorted in alphabetical order according to either the emitting atoms (default) or the wavelength (using the crit='wave' parameter):
In [33]:
for line in obs.getSortedLines():
print(line.label, line.corrIntens[0])
The following method, which gives the list of all the atoms implied in the observed emission lines, will be useful later:
In [34]:
atomList = obs.getUniqueAtoms()
In [35]:
atomList
Out[35]:
Once an Observation object is instantiated, you can add a new observation (corresponding, e.g., to a new object or a new fiber) by using:
In [36]:
obs.addObs('test', np.random.rand(25))
where ‘test’ is the name of the new observation. The new observation must have the same size of obs, that is, it must contain obs.n_lines lines.
In [37]:
obs.printIntens()
In [38]:
line = pn.EmissionLine(label='Cl3_5518A', obsIntens=[3.5, 2.5])
obs.addLine(line)
In [39]:
obs.printIntens()
You can extract the line intensities from an Observation object by, for example:
In [40]:
obs.names
Out[40]:
In [41]:
obs.getIntens(obsName='SMC_24')
Out[41]:
In [42]:
obs.getIntens()['O2_7330A+']
Out[42]:
Once the electron temperature and density are determined, it is easy to obtain the ionic abundances from a set of emission lines included in an Observation object:
In [43]:
obs = pn.Observation()
obs.readData('observations1.dat', fileFormat='lines_in_rows', err_default=0.05) # fill obs with data read from observations1.dat
obs.def_EBV(label1="H1r_6563A", label2="H1r_4861A", r_theo=2.85)
obs.correctData(normWave=4861.)
Te = [10000.]
Ne = [1e3]
# Define a dictionary to hold all the Atom objects needed
all_atoms = pn.getAtomDict(atom_list=obs.getUniqueAtoms())
# define a dictionary to store the abundances
ab_dict = {}
# we use the following lines to determine the ionic abundances
ab_labels = ['N2_6584A', 'O2_3726A', 'O3_5007A', 'S2_6716A',
'S3_6312A', 'Ar3_7136A', 'Ne3_3869A']
for line in obs.getSortedLines():
if line.label in ab_labels:
ab = all_atoms[line.atom].getIonAbundance(line.corrIntens, Te, Ne,
to_eval=line.to_eval, Hbeta=100)
ab_dict[line.atom] = ab
In [44]:
ab_dict
Out[44]: