In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import pyneb as pn
Most plots are obtained by operating on emission maps, which are grids of emissivities as a function of temperature and density generated by the EmisGrid class.
EmisGrid instantiates an atom and computes the emissivities of all its lines for the (tem, den) values of a regularly spaced grid (may be log or linear in the case of the density). Each line is represented in a 2D array (a grid), and there are as many arrays transitions in the atom. The results can be operated on, saved for a later use in a cPickle file, or restored.
The following command instantiates an [O III] atom and computes the emissivity of all its lines in a 30x30 grid:
In [2]:
O3_EG = pn.EmisGrid('O', 3, n_tem=30, n_den=30)
The arguments are described in more details in the Reference Manual. Here is the list:
In [3]:
O3_EG = pn.EmisGrid(elem='O', spec=3, n_tem=100, n_den=100,
tem_min=5000., tem_max=20000., den_min=10.,
den_max=1.e8, restore_file=None, atomObj=None)
The emissivity grid of a specific line can be obtained by means of:
In [4]:
O3_5007 = O3_EG.getGrid(wave=5007)
In [5]:
O3_5007.shape
Out[5]:
The emissivity grid of a combination of lines can also be computed:
In [6]:
O3_Te = O3_EG.getGrid(to_eval = 'L(4363)/L(5007)')
There are two plotting tools integrated in the EmisGrid object:
In [7]:
O3_EG.plotImage(to_eval = 'L(4363)/L(5007)')
In [8]:
f, ax = plt.subplots(figsize=(7,5))
O3_EG.plotContours(to_eval = 'L(4363)/L(5007)', ax=ax)
f.savefig('OIII_diag.pdf')
See the Reference Manual for more option and the Diagnostic class for producing plots combining different atoms.
It is quite common to have to instantiate various EmisGrid objects, especially if you want to make a diagnostic diagram. This can easily be done using the getEmisGridDict method, used for example as follows:
In [9]:
emisgrids = pn.getEmisGridDict(atom_list=['O2', 'O3', 'N2'])
In [10]:
emisgrids
Out[10]:
This command generates a dictionary of emission grids for [O II], [O III] and [N II]. The resulting maps are saved in a directory defined by default when PyNeb is started, in the pn.config.pypic_path variable. It first tries to use the $HOME/.pypics directory; if it fails, it tries to use /tmp/pypic; if it fails too, the value is set to None and a user-defined value has to be provided by changing pn.config.pypic_path or using the pypic_path parameter when calling getEmisGridDict.
If a Diagnostic object is already available (see next Section), it can be used to determine the relevant atoms for which a grid must be computed or restored:
In [11]:
diags = pn.Diagnostics() # See next section
emisgrids = pn.getEmisGridDict(atomDict=diags.atomDict, den_max=1e6)
This EmisGrid dictionary will be very useful to plot diagnostic diagrams with the Diagnostic object, as is described in the next section.
Diagnostics is the class used to evaluate temperatures and densities from line ratios. It is also the class that plots the diagnostic Te-Ne diagrams. The object is instantiated like this:
In [12]:
diags = pn.Diagnostics() # instantiate the Diagnostic class
An optional parameter addAll=True (default is False) lets the object load all the available diagnostics. Most of the time this option is not used and the diagnostics are added manually as they are needed:
In [13]:
diags.addDiag(['[NI] 5198/5200','[NII] 5755/6548','[OII] 3726/3729'])
Each diagnostic is defined by a label and is associated to a tuple containing 3 elements: the atom corresponding to the diagnostic lines, the algebraic definition of the line ratios and the algebraic definition of the error of the diagnostic, which depends on the error of each line involved. In the present case, the diagnostic is the ratios of two [O III] lines, 4363/5007, and the error is the quadratic sum of the relative error of each line (E(lambda)): RMS(a, b) = sqrt(a2 + b2).
Users can also define their own diagnostics, for example using:
In [14]:
diags.addDiag('[OIII] 4363/4959', ('O3', 'L(4363)/L(4959)', 'RMS([E(4363),E(4959)])'))
Notice that the diagnostics are defined so that they tend to increase with the main parameter they trace: [OIII] 4363/5007 increases with the electron temperature.
The diagnostics contained in a Diagnostics object are listed by means of the diags.getDiagLabels() and diags.getDiags() methods. Once added to the Diagnostics object, they can be used either to compute Te and Ne via getCrossTemDen or to plot diagrams (see below). A diagnostic can be removed from the list with the delDiag method.
The getCrossTemDen method cross-converges the temperature and density derived from two sensitive line ratios (not necessarily from the same atom), by inputting the quantity derived with one line ratio into the other and then iterating. When the iteration process ends, the two diagnostics give self-consistent results. The first line ratio must be a temperature-sensitive one and the second a density-sensitive one. The temperature and density can be individual numbers as well as arrays (provided they are equal in shape).
In [15]:
diags.getCrossTemDen('[NII] 5755/6548', '[SII] 6731/6716', 0.02, 1.0)
Out[15]:
In [16]:
for diag in sorted(pn.diags_dict.keys()):
print('"{0}" : {1}'.format(diag, pn.diags_dict[diag]))
The plotting tool included in the Diagnostics class requires an EmisGrid dictionary (as returned by pn.getEmisGridDict; see the previous section) and an Observation object. The plot is obtained by:
In [17]:
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.)
In [18]:
diags = pn.Diagnostics()
diags.addDiagsFromObs(obs)
diags.diags
Out[18]:
In [19]:
obs.getIntens()
Out[19]:
In [20]:
import matplotlib as mpl
%config InlineBackend.figure_format = 'png'
mpl.rc("savefig", dpi=150)
emisgrids = pn.getEmisGridDict(atomDict=diags.atomDict)
diags.plot(emisgrids, obs)
If there is more than one spectrum in the Observation object, the index of the observation to be used is given by i_obs:
In [21]:
diags.plot(emisgrids, obs, i_obs = 0)
In [22]:
diags.diags
Out[22]:
The label used to identify emission lines can be changed in the diags object, using for example:
In [23]:
diags.addClabel('[OIII] 4363/5007', '[OIII]na')
In [24]:
diags.plot(emisgrids, obs, i_obs = 0)