[yt-svn] commit/yt: brittonsmith: Merged in jzuhone/yt-3.x/yt-3.0 (pull request #963)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Jul 2 23:57:38 PDT 2014
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/c3d17397f8c1/
Changeset: c3d17397f8c1
Branch: yt-3.0
User: brittonsmith
Date: 2014-07-03 08:57:30
Summary: Merged in jzuhone/yt-3.x/yt-3.0 (pull request #963)
Porting spectral integrator to 3.0
Affected #: 10 files
diff -r 57ecc4cbdf24f36a5d522c983612468861ffde35 -r c3d17397f8c1f1f82df9d7d255e3e8f8964df57d doc/install_script.sh
--- a/doc/install_script.sh
+++ b/doc/install_script.sh
@@ -567,8 +567,10 @@
mkdir -p ${DEST_DIR}/data
cd ${DEST_DIR}/data
-echo 'de6d8c6ea849f0206d219303329a0276b3cce7c051eec34377d42aacbe0a4f47ac5145eb08966a338ecddd2b83c8f787ca9956508ad5c39ee2088ad875166410 xray_emissivity.h5' > xray_emissivity.h5.sha512
-[ ! -e xray_emissivity.h5 ] && get_ytdata xray_emissivity.h5
+echo 'de6d8c6ea849f0206d219303329a0276b3cce7c051eec34377d42aacbe0a4f47ac5145eb08966a338ecddd2b83c8f787ca9956508ad5c39ee2088ad875166410 cloudy_emissivity.h5' > cloudy_emissivity.h5.sha512
+[ ! -e cloudy_emissivity.h5 ] && get_ytdata cloudy_emissivity.h5
+echo '0f714ae2eace0141b1381abf1160dc8f8a521335e886f99919caf3beb31df1fe271d67c7b2a804b1467949eb16b0ef87a3d53abad0e8160fccac1e90d8d9e85f apec_emissivity.h5' > apec_emissivity.h5.sha512
+[ ! -e apec_emissivity.h5 ] && get_ytdata apec_emissivity.h5
# Set paths to what they should be when yt is activated.
export PATH=${DEST_DIR}/bin:$PATH
diff -r 57ecc4cbdf24f36a5d522c983612468861ffde35 -r c3d17397f8c1f1f82df9d7d255e3e8f8964df57d doc/source/analyzing/analysis_modules/photon_simulator.rst
--- a/doc/source/analyzing/analysis_modules/photon_simulator.rst
+++ b/doc/source/analyzing/analysis_modules/photon_simulator.rst
@@ -1,6 +1,11 @@
Constructing Mock X-ray Observations
------------------------------------
+.. note::
+
+ If you just want to create derived fields for X-ray emission,
+ you should go `here <xray_emission_fields.html>`_ instead.
+
The ``photon_simulator`` analysis module enables the creation of
simulated X-ray photon lists of events from datasets that ``yt`` is able
to read. The simulated events then can be exported to X-ray telescope
diff -r 57ecc4cbdf24f36a5d522c983612468861ffde35 -r c3d17397f8c1f1f82df9d7d255e3e8f8964df57d doc/source/analyzing/analysis_modules/xray_emission_fields.rst
--- a/doc/source/analyzing/analysis_modules/xray_emission_fields.rst
+++ b/doc/source/analyzing/analysis_modules/xray_emission_fields.rst
@@ -2,41 +2,46 @@
X-ray Emission Fields
=====================
-.. sectionauthor:: Britton Smith <brittonsmith at gmail.com>
+.. sectionauthor:: Britton Smith <brittonsmith at gmail.com>, John ZuHone <jzuhone at gmail.com>
+
+.. note::
+
+ If you came here trying to figure out how to create simulated X-ray photons and observations,
+ you should go `here <photon_simulator.html>`_ instead.
This functionality provides the ability to create metallicity-dependent
-X-ray luminosity, emissivity, and photo emissivity fields for a given
+X-ray luminosity, emissivity, and photon emissivity fields for a given
photon energy range. This works by interpolating from emission tables
-created with the photoionization code, `Cloudy <http://nublado.org/>`_.
-If you installed yt with the install script, the data should be located in
-the *data* directory inside the installation directory. Emission fields can
-be made for any interval between 0.1 keV and 100 keV.
+created from the photoionization code `Cloudy <http://nublado.org/>`_ or
+the collisional ionization database `AtomDB <http://www.atomdb.org>`_. If
+you installed yt with the install script, these data files should be located in
+the *data* directory inside the installation directory, or can be downloaded
+from `<http://yt-project.org/data>`_. Emission fields can be made for any
+interval between 0.1 keV and 100 keV.
Adding Emission Fields
----------------------
-Fields can be created for luminosity (erg/s), emissivity (erg/s/cm^3),
-and photon emissivity (photons/s/cm^3). The only required arguments are
-the minimum and maximum energies.
+Fields will be created for luminosity :math:`{\rm (erg~s^{-1})}`, emissivity :math:`{\rm (erg~s^{-1}~cm^{-3})}`,
+and photon emissivity :math:`{\rm (photons~s^{-1}~cm^{-3})}`. The only required arguments are the
+dataset object, and the minimum and maximum energies of the energy band.
.. code-block:: python
- from yt.mods import *
+ import yt
from yt.analysis_modules.spectral_integrator.api import \
- add_xray_luminosity_field, \
- add_xray_emissivity_field, \
- add_xray_photon_emissivity_field
+ add_xray_emissivity_field
- add_xray_luminosity_field(0.5, 7)
- add_xray_emissivity_field(0.5, 7)
- add_xray_photon_emissivity_field(0.5, 7)
+ xray_fields = add_xray_emissivity_field(0.5, 7.0)
Additional keyword arguments are:
- * **filename** (*string*): Path to data file containing emissivity
- values. If None, a file called xray_emissivity.h5 is used. This file
- contains emissivity tables for primordial elements and for metals at
- solar metallicity for the energy range 0.1 to 100 keV. Default: None.
+ * **filename** (*string*): Path to data file containing emissivity values. If None,
+ a file called "cloudy_emissivity.h5" is used, for photoionized plasmas. A second
+ option, for collisionally ionized plasmas, is in the file "apec_emissivity.h5",
+ available at http://yt-project.org/data. These files contain emissivity tables
+ for primordial elements and for metals at solar metallicity for the energy range
+ 0.1 to 100 keV. Default: None.
* **with_metals** (*bool*): If True, use the metallicity field to add the
contribution from metals. If False, only the emission from H/He is
@@ -46,24 +51,27 @@
metallicity for the emission from metals. The *with_metals* keyword
must be set to False to use this. Default: None.
-The resulting fields can be used like all normal fields.
+The resulting fields can be used like all normal fields. The function will return the names of
+the created fields in a Python list.
-.. python-script::
+.. code-block:: python
- from yt.mods import *
+ import yt
from yt.analysis_modules.spectral_integrator.api import \
- add_xray_luminosity_field, \
- add_xray_emissivity_field, \
- add_xray_photon_emissivity_field
+ add_xray_emissivity_field
- add_xray_luminosity_field(0.5, 7)
- add_xray_emissivity_field(0.5, 7)
- add_xray_photon_emissivity_field(0.5, 7)
+ xray_fields = add_xray_emissivity_field(0.5, 7.0, filename="apec_emissivity.h5")
- pf = load("enzo_tiny_cosmology/DD0046/DD0046")
- plot = SlicePlot(pf, 'x', 'Xray_Luminosity_0.5_7keV')
+ ds = yt.load("enzo_tiny_cosmology/DD0046/DD0046")
+ plot = yt.SlicePlot(ds, 'x', 'xray_luminosity_0.5_7.0_keV')
plot.save()
- plot = ProjectionPlot(pf, 'x', 'Xray_Emissivity_0.5_7keV')
+ plot = yt.ProjectionPlot(ds, 'x', 'xray_emissivity_0.5_7.0_keV')
plot.save()
- plot = ProjectionPlot(pf, 'x', 'Xray_Photon_Emissivity_0.5_7keV')
+ plot = yt.ProjectionPlot(ds, 'x', 'xray_photon_emissivity_0.5_7.0_keV')
plot.save()
+
+.. warning::
+
+ The X-ray fields depend on the number density of hydrogen atoms, in the yt field
+ ``H_number_density``. If this field is not defined (either in the dataset or by the user),
+ the primordial hydrogen mass fraction (X = 0.76) will be used to construct it.
\ No newline at end of file
diff -r 57ecc4cbdf24f36a5d522c983612468861ffde35 -r c3d17397f8c1f1f82df9d7d255e3e8f8964df57d yt/analysis_modules/photon_simulator/photon_simulator.py
--- a/yt/analysis_modules/photon_simulator/photon_simulator.py
+++ b/yt/analysis_modules/photon_simulator/photon_simulator.py
@@ -956,7 +956,7 @@
if isinstance(self.parameters["Area"], basestring):
mylog.error("Writing SIMPUT files is only supported if you didn't convolve with an ARF.")
- raise TypeError
+ raise TypeError("Writing SIMPUT files is only supported if you didn't convolve with an ARF.")
if emin is None:
emin = self.events["eobs"].min().value*0.95
@@ -1032,7 +1032,10 @@
f = h5py.File(h5file, "w")
f.create_dataset("/exp_time", data=float(self.parameters["ExposureTime"]))
- f.create_dataset("/area", data=float(self.parameters["Area"]))
+ area = self.parameters["Area"]
+ if not isinstance(area, basestring):
+ area = float(area)
+ f.create_dataset("/area", data=area)
f.create_dataset("/redshift", data=self.parameters["Redshift"])
f.create_dataset("/d_a", data=float(self.parameters["AngularDiameterDistance"]))
if "ARF" in self.parameters:
diff -r 57ecc4cbdf24f36a5d522c983612468861ffde35 -r c3d17397f8c1f1f82df9d7d255e3e8f8964df57d yt/analysis_modules/photon_simulator/spectral_models.py
--- a/yt/analysis_modules/photon_simulator/spectral_models.py
+++ b/yt/analysis_modules/photon_simulator/spectral_models.py
@@ -181,7 +181,7 @@
Examples
--------
>>> apec_model = TableApecModel("/Users/jzuhone/Data/atomdb_v2.0.2/", 0.05, 50.0,
- 1000, thermal_broad=True)
+ ... 1000, thermal_broad=True)
"""
def __init__(self, apec_root, emin, emax, nchan,
apec_vers="2.0.2", thermal_broad=False):
@@ -226,7 +226,7 @@
tmpspec = np.zeros((self.nchan))
- i = np.where((self.line_handle[tindex].data.field('element')==element) &
+ i = np.where((self.line_handle[tindex].data.field('element') == element) &
(self.line_handle[tindex].data.field('lambda') > self.minlam) &
(self.line_handle[tindex].data.field('lambda') < self.maxlam))[0]
@@ -242,24 +242,24 @@
vec += np.diff(cdf)*a
else:
ie = np.searchsorted(ebins, E0, side='right')-1
- for i,a in zip(ie,amp): vec[i] += a
+ for i, a in zip(ie, amp): vec[i] += a
tmpspec += vec
- ind = np.where((self.coco_handle[tindex].data.field('Z')==element) &
- (self.coco_handle[tindex].data.field('rmJ')==0))[0]
- if len(ind)==0:
+ ind = np.where((self.coco_handle[tindex].data.field('Z') == element) &
+ (self.coco_handle[tindex].data.field('rmJ') == 0))[0]
+ if len(ind) == 0:
return tmpspec
else:
- ind=ind[0]
+ ind = ind[0]
- n_cont=self.coco_handle[tindex].data.field('N_Cont')[ind]
- e_cont=self.coco_handle[tindex].data.field('E_Cont')[ind][:n_cont]
+ n_cont = self.coco_handle[tindex].data.field('N_Cont')[ind]
+ e_cont = self.coco_handle[tindex].data.field('E_Cont')[ind][:n_cont]
continuum = self.coco_handle[tindex].data.field('Continuum')[ind][:n_cont]
tmpspec += np.interp(self.emid.ndarray_view(), e_cont, continuum)*self.de.ndarray_view()
- n_pseudo=self.coco_handle[tindex].data.field('N_Pseudo')[ind]
- e_pseudo=self.coco_handle[tindex].data.field('E_Pseudo')[ind][:n_pseudo]
+ n_pseudo = self.coco_handle[tindex].data.field('N_Pseudo')[ind]
+ e_pseudo = self.coco_handle[tindex].data.field('E_Pseudo')[ind][:n_pseudo]
pseudo = self.coco_handle[tindex].data.field('Pseudo')[ind][:n_pseudo]
tmpspec += np.interp(self.emid.ndarray_view(), e_pseudo, pseudo)*self.de.ndarray_view()
diff -r 57ecc4cbdf24f36a5d522c983612468861ffde35 -r c3d17397f8c1f1f82df9d7d255e3e8f8964df57d yt/analysis_modules/spectral_integrator/api.py
--- a/yt/analysis_modules/spectral_integrator/api.py
+++ b/yt/analysis_modules/spectral_integrator/api.py
@@ -15,6 +15,4 @@
from .spectral_frequency_integrator import \
EmissivityIntegrator, \
- add_xray_emissivity_field, \
- add_xray_luminosity_field, \
- add_xray_photon_emissivity_field
+ add_xray_emissivity_field
\ No newline at end of file
diff -r 57ecc4cbdf24f36a5d522c983612468861ffde35 -r c3d17397f8c1f1f82df9d7d255e3e8f8964df57d yt/analysis_modules/spectral_integrator/spectral_frequency_integrator.py
--- a/yt/analysis_modules/spectral_integrator/spectral_frequency_integrator.py
+++ b/yt/analysis_modules/spectral_integrator/spectral_frequency_integrator.py
@@ -1,6 +1,6 @@
"""
Integrator classes to deal with interpolation and integration of input spectral
-bins. Currently only supports Cloudy-style data.
+bins. Currently only supports Cloudy and APEC-style data.
@@ -23,19 +23,21 @@
mylog, \
only_on_root
-from yt.fields.local_fields import add_field
+from yt.utilities.exceptions import YTFieldNotFound
from yt.utilities.exceptions import YTException
from yt.utilities.linear_interpolators import \
- BilinearFieldInterpolator
+ UnilinearFieldInterpolator, BilinearFieldInterpolator
from yt.utilities.physical_constants import \
- erg_per_eV, hcgs
-from yt.units import keV, Hz
-keV_per_Hz = keV/Hz/hcgs
+ hcgs, mp
+from yt.units.yt_array import YTArray, YTQuantity
+from yt.utilities.physical_ratios import \
+ primordial_H_mass_fraction, erg_per_keV
xray_data_version = 1
-def _get_data_file():
- data_file = "xray_emissivity.h5"
+def _get_data_file(data_file=None):
+ if data_file is None:
+ data_file = "cloudy_emissivity.h5"
data_url = "http://yt-project.org/data"
if "YT_DEST" in os.environ and \
os.path.isdir(os.path.join(os.environ["YT_DEST"], "data")):
@@ -62,8 +64,9 @@
class ObsoleteDataException(YTException):
def __str__(self):
- return "X-ray emissivity data is out of data.\nDownload the latest data from http://yt-project.org/data/xray_emissivity.h5 and move it to %s." % \
- os.path.join(os.environ["YT_DEST"], "data", "xray_emissivity.h5")
+ return "X-ray emissivity data is out of date.\n" + \
+ "Download the latest data from http://yt-project.org/data/cloudy_emissivity.h5 and move it to %s." % \
+ os.path.join(os.environ["YT_DEST"], "data", "cloudy_emissivity.h5")
class EmissivityIntegrator(object):
r"""Class for making X-ray emissivity fields with hdf5 data tables
@@ -75,9 +78,11 @@
----------
filename: string, default None
Path to data file containing emissivity values. If None,
- a file called xray_emissivity.h5 is used. This file contains
- emissivity tables for primordial elements and for metals at
- solar metallicity for the energy range 0.1 to 100 keV.
+ a file called "cloudy_emissivity.h5" is used, for photoionized
+ plasmas. A second option, for collisionally ionized plasmas, is
+ in the file "apec_emissivity.h5", available at http://yt-project.org/data.
+ These files contain emissivity tables for primordial elements and
+ for metals at solar metallicity for the energy range 0.1 to 100 keV.
Default: None.
"""
@@ -89,7 +94,8 @@
default_filename = True
if not os.path.exists(filename):
- raise IOError("File does not exist: %s." % filename)
+ mylog.warning("File %s does not exist, will attempt to find it." % filename)
+ filename = _get_data_file(data_file=filename)
only_on_root(mylog.info, "Loading emissivity data from %s." % filename)
in_file = h5py.File(filename, "r")
if "info" in in_file.attrs:
@@ -103,51 +109,51 @@
for field in ["emissivity_primordial", "emissivity_metals",
"log_nH", "log_T", "log_E"]:
- setattr(self, field, in_file[field][:])
+ if field in in_file:
+ setattr(self, field, in_file[field][:])
in_file.close()
E_diff = np.diff(self.log_E)
self.E_bins = \
- np.power(10, np.concatenate([self.log_E[:-1] - 0.5 * E_diff,
- [self.log_E[-1] - 0.5 * E_diff[-1],
- self.log_E[-1] + 0.5 * E_diff[-1]]]))
- self.dnu = keV_per_Hz * np.diff(self.E_bins)
+ YTArray(np.power(10, np.concatenate([self.log_E[:-1] - 0.5 * E_diff,
+ [self.log_E[-1] - 0.5 * E_diff[-1],
+ self.log_E[-1] + 0.5 * E_diff[-1]]])),
+ "keV")
+ self.dnu = (np.diff(self.E_bins)/hcgs).in_units("Hz")
- def _get_interpolator(self, data, e_min, e_max):
- r"""Create an interpolator for total emissivity in a
- given energy range.
-
- Parameters
- ----------
- e_min: float
- the minimum energy in keV for the energy band.
- e_min: float
- the maximum energy in keV for the energy band.
-
- """
+ def get_interpolator(self, data, e_min, e_max):
+ e_min = YTQuantity(e_min, "keV")
+ e_max = YTQuantity(e_max, "keV")
if (e_min - self.E_bins[0]) / e_min < -1e-3 or \
(e_max - self.E_bins[-1]) / e_max > 1e-3:
- raise EnergyBoundsException(np.power(10, self.E_bins[0]),
- np.power(10, self.E_bins[-1]))
+ raise EnergyBoundsException(self.E_bins[0], self.E_bins[-1])
e_is, e_ie = np.digitize([e_min, e_max], self.E_bins)
e_is = np.clip(e_is - 1, 0, self.E_bins.size - 1)
e_ie = np.clip(e_ie, 0, self.E_bins.size - 1)
- my_dnu = np.copy(self.dnu[e_is: e_ie])
+ my_dnu = self.dnu[e_is: e_ie].copy()
# clip edge bins if the requested range is smaller
- my_dnu[0] -= e_min - self.E_bins[e_is]
- my_dnu[-1] -= self.E_bins[e_ie] - e_max
+ my_dnu[0] -= ((e_min - self.E_bins[e_is])/hcgs).in_units("Hz")
+ my_dnu[-1] -= ((self.E_bins[e_ie] - e_max)/hcgs).in_units("Hz")
interp_data = (data[..., e_is:e_ie] * my_dnu).sum(axis=-1)
- return BilinearFieldInterpolator(np.log10(interp_data),
- [self.log_nH[0], self.log_nH[-1],
- self.log_T[0], self.log_T[-1]],
- ["log_nH", "log_T"], truncate=True)
+ if len(data.shape) == 2:
+ emiss = UnilinearFieldInterpolator(np.log10(interp_data),
+ [self.log_T[0], self.log_T[-1]],
+ "log_T", truncate=True)
+ else:
+ emiss = BilinearFieldInterpolator(np.log10(interp_data),
+ [self.log_nH[0], self.log_nH[-1],
+ self.log_T[0], self.log_T[-1]],
+ ["log_nH", "log_T"], truncate=True)
-def add_xray_emissivity_field(e_min, e_max, filename=None,
+ return emiss
+
+def add_xray_emissivity_field(ds, e_min, e_max,
+ filename=None,
with_metals=True,
constant_metallicity=None):
- r"""Create an X-ray emissivity field for a given energy range.
+ r"""Create X-ray emissivity fields for a given energy range.
Parameters
----------
@@ -155,197 +161,117 @@
the minimum energy in keV for the energy band.
e_min: float
the maximum energy in keV for the energy band.
-
- Other Parameters
- ----------------
- filename: string
+ filename: string, optional
Path to data file containing emissivity values. If None,
- a file called xray_emissivity.h5 is used. This file contains
- emissivity tables for primordial elements and for metals at
- solar metallicity for the energy range 0.1 to 100 keV.
+ a file called "cloudy_emissivity.h5" is used, for photoionized
+ plasmas. A second option, for collisionally ionized plasmas, is
+ in the file "apec_emissivity.h5", available at http://yt-project.org/data.
+ These files contain emissivity tables for primordial elements and
+ for metals at solar metallicity for the energy range 0.1 to 100 keV.
Default: None.
- with_metals: bool
+ with_metals: bool, optional
If True, use the metallicity field to add the contribution from
metals. If False, only the emission from H/He is considered.
Default: True.
- constant_metallicity: float
+ constant_metallicity: float, optional
If specified, assume a constant metallicity for the emission
from metals. The *with_metals* keyword must be set to False
to use this.
Default: None.
- This will create a field named "Xray_Emissivity_{e_min}_{e_max}keV".
- The units of the field are erg s^-1 cm^-3.
+ This will create three fields:
+
+ "xray_emissivity_{e_min}_{e_max}_keV" (erg s^-1 cm^-3)
+ "xray_luminosity_{e_min}_{e_max}_keV" (erg s^-1)
+ "xray_photon_emissivity_{e_min}_{e_max}_keV" (photons s^-1 cm^-3)
Examples
--------
>>> from yt.mods import *
>>> from yt.analysis_modules.spectral_integrator.api import *
- >>> add_xray_emissivity_field(0.5, 2)
- >>> pf = load(dataset)
- >>> p = ProjectionPlot(pf, 'x', "Xray_Emissivity_0.5_2keV")
+ >>> ds = load(dataset)
+ >>> add_xray_emissivity_field(ds, 0.5, 2)
+ >>> p = ProjectionPlot(ds, 'x', "xray_emissivity_0.5_2_keV")
>>> p.save()
"""
+ if with_metals:
+ try:
+ ds._get_field_info("metal_density")
+ except YTFieldNotFound:
+ raise RuntimeError("Your dataset does not have a \"metal_density\" field! " +
+ "Perhaps you should specify a constant metallicity?")
+
my_si = EmissivityIntegrator(filename=filename)
- em_0 = my_si._get_interpolator(my_si.emissivity_primordial, e_min, e_max)
+ em_0 = my_si.get_interpolator(my_si.emissivity_primordial, e_min, e_max)
em_Z = None
if with_metals or constant_metallicity is not None:
- em_Z = my_si._get_interpolator(my_si.emissivity_metals, e_min, e_max)
+ em_Z = my_si.get_interpolator(my_si.emissivity_metals, e_min, e_max)
+
+ energy_erg = np.power(10, my_si.log_E) * erg_per_keV
+ emp_0 = my_si.get_interpolator((my_si.emissivity_primordial[..., :] / energy_erg),
+ e_min, e_max)
+ emp_Z = None
+ if with_metals or constant_metallicity is not None:
+ emp_Z = my_si.get_interpolator((my_si.emissivity_metals[..., :] / energy_erg),
+ e_min, e_max)
+
+ try:
+ ds._get_field_info("H_number_density")
+ except YTFieldNotFound:
+ mylog.warning("Could not find a field for \"H_number_density\". Assuming primordial H " +
+ "mass fraction.")
+ def _nh(field, data):
+ return primordial_H_mass_fraction*data["gas","density"]/mp
+ ds.add_field(("gas","H_number_density"), function=_nh, units="cm**-3")
def _emissivity_field(field, data):
- dd = {"log_nH" : np.log10(data["H_NumberDensity"]),
- "log_T" : np.log10(data["Temperature"])}
+ dd = {"log_nH" : np.log10(data["gas","H_number_density"]),
+ "log_T" : np.log10(data["gas","temperature"])}
my_emissivity = np.power(10, em_0(dd))
if em_Z is not None:
if with_metals:
- my_Z = data["Metallicity"]
+ my_Z = data["gas","metallicity"]
elif constant_metallicity is not None:
my_Z = constant_metallicity
my_emissivity += my_Z * np.power(10, em_Z(dd))
- return data["H_NumberDensity"]**2 * my_emissivity
+ return data["gas","H_number_density"]**2 * YTArray(my_emissivity, "erg*cm**3/s")
- field_name = "Xray_Emissivity_%s_%skeV" % (e_min, e_max)
- add_field(field_name, function=_emissivity_field,
- projection_conversion="cm",
- display_name=r"\epsilon_{X}\/(%s-%s\/keV)" % (e_min, e_max),
- units=r"\rm{erg}\/\rm{cm}^{-3}\/\rm{s}^{-1}")
- return field_name
-
-def add_xray_luminosity_field(e_min, e_max, filename=None,
- with_metals=True,
- constant_metallicity=None):
- r"""Create an X-ray luminosity field for a given energy range.
-
- Parameters
- ----------
- e_min: float
- the minimum energy in keV for the energy band.
- e_min: float
- the maximum energy in keV for the energy band.
-
- Other Parameters
- ----------------
- filename: string
- Path to data file containing emissivity values. If None,
- a file called xray_emissivity.h5 is used. This file contains
- emissivity tables for primordial elements and for metals at
- solar metallicity for the energy range 0.1 to 100 keV.
- Default: None.
- with_metals: bool
- If True, use the metallicity field to add the contribution from
- metals. If False, only the emission from H/He is considered.
- Default: True.
- constant_metallicity: float
- If specified, assume a constant metallicity for the emission
- from metals. The *with_metals* keyword must be set to False
- to use this.
- Default: None.
-
- This will create a field named "Xray_Luminosity_{e_min}_{e_max}keV".
- The units of the field are erg s^-1.
-
- Examples
- --------
-
- >>> from yt.mods import *
- >>> from yt.analysis_modules.spectral_integrator.api import *
- >>> add_xray_luminosity_field(0.5, 2)
- >>> pf = load(dataset)
- >>> sp = pf.sphere('max', (2., 'mpc'))
- >>> print sp.quantities['TotalQuantity']('Xray_Luminosity_0.5_2keV')
-
- """
-
- em_field = add_xray_emissivity_field(e_min, e_max, filename=filename,
- with_metals=with_metals,
- constant_metallicity=constant_metallicity)
+ emiss_name = "xray_emissivity_%s_%s_keV" % (e_min, e_max)
+ ds.add_field(emiss_name, function=_emissivity_field,
+ display_name=r"\epsilon_{X}\/(%s-%s\/keV)" % (e_min, e_max),
+ units="erg/cm**3/s")
def _luminosity_field(field, data):
- return data[em_field] * data["CellVolume"]
- field_name = "Xray_Luminosity_%s_%skeV" % (e_min, e_max)
- add_field(field_name, function=_luminosity_field,
- display_name=r"\rm{L}_{X}\/(%s-%s\/keV)" % (e_min, e_max),
- units=r"\rm{erg}\/\rm{s}^{-1}")
- return field_name
+ return data[emiss_name] * data["cell_volume"]
-def add_xray_photon_emissivity_field(e_min, e_max, filename=None,
- with_metals=True,
- constant_metallicity=None):
- r"""Create an X-ray photon emissivity field for a given energy range.
+ lum_name = "xray_luminosity_%s_%s_keV" % (e_min, e_max)
+ ds.add_field(lum_name, function=_luminosity_field,
+ display_name=r"\rm{L}_{X}\/(%s-%s\/keV)" % (e_min, e_max),
+ units="erg/s")
- Parameters
- ----------
- e_min: float
- the minimum energy in keV for the energy band.
- e_min: float
- the maximum energy in keV for the energy band.
+ def _photon_emissivity_field(field, data):
+ dd = {"log_nH" : np.log10(data["gas","H_number_density"]),
+ "log_T" : np.log10(data["gas","temperature"])}
- Other Parameters
- ----------------
- filename: string
- Path to data file containing emissivity values. If None,
- a file called xray_emissivity.h5 is used. This file contains
- emissivity tables for primordial elements and for metals at
- solar metallicity for the energy range 0.1 to 100 keV.
- Default: None.
- with_metals: bool
- If True, use the metallicity field to add the contribution from
- metals. If False, only the emission from H/He is considered.
- Default: True.
- constant_metallicity: float
- If specified, assume a constant metallicity for the emission
- from metals. The *with_metals* keyword must be set to False
- to use this.
- Default: None.
-
- This will create a field named "Xray_Photon_Emissivity_{e_min}_{e_max}keV".
- The units of the field are photons s^-1 cm^-3.
-
- Examples
- --------
-
- >>> from yt.mods import *
- >>> from yt.analysis_modules.spectral_integrator.api import *
- >>> add_xray_emissivity_field(0.5, 2)
- >>> pf = load(dataset)
- >>> p = ProjectionPlot(pf, 'x', "Xray_Emissivity_0.5_2keV")
- >>> p.save()
-
- """
-
- my_si = EmissivityIntegrator(filename=filename)
- energy_erg = np.power(10, my_si.log_E) * erg_per_eV
-
- em_0 = my_si._get_interpolator((my_si.emissivity_primordial[..., :] / energy_erg),
- e_min, e_max)
- em_Z = None
- if with_metals or constant_metallicity is not None:
- em_Z = my_si._get_interpolator((my_si.emissivity_metals[..., :] / energy_erg),
- e_min, e_max)
-
- def _emissivity_field(field, data):
- dd = {"log_nH" : np.log10(data["H_NumberDensity"]),
- "log_T" : np.log10(data["Temperature"])}
-
- my_emissivity = np.power(10, em_0(dd))
- if em_Z is not None:
+ my_emissivity = np.power(10, emp_0(dd))
+ if emp_Z is not None:
if with_metals:
- my_Z = data["Metallicity"]
+ my_Z = data["gas","metallicity"]
elif constant_metallicity is not None:
my_Z = constant_metallicity
- my_emissivity += my_Z * np.power(10, em_Z(dd))
+ my_emissivity += my_Z * np.power(10, emp_Z(dd))
- return data["H_NumberDensity"]**2 * my_emissivity
+ return data["gas","H_number_density"]**2 * YTArray(my_emissivity, "photons*cm**3/s")
- field_name = "Xray_Photon_Emissivity_%s_%skeV" % (e_min, e_max)
- add_field(field_name, function=_emissivity_field,
- projection_conversion="cm",
- display_name=r"\epsilon_{X}\/(%s-%s\/keV)" % (e_min, e_max),
- units=r"\rm{photons}\/\rm{cm}^{-3}\/\rm{s}^{-1}")
- return field_name
+ phot_name = "xray_photon_emissivity_%s_%s_keV" % (e_min, e_max)
+ ds.add_field(phot_name, function=_photon_emissivity_field,
+ display_name=r"\epsilon_{X}\/(%s-%s\/keV)" % (e_min, e_max),
+ units="photons/cm**3/s")
+
+ return emiss_name, lum_name, phot_name
\ No newline at end of file
diff -r 57ecc4cbdf24f36a5d522c983612468861ffde35 -r c3d17397f8c1f1f82df9d7d255e3e8f8964df57d yt/units/unit_lookup_table.py
--- a/yt/units/unit_lookup_table.py
+++ b/yt/units/unit_lookup_table.py
@@ -89,6 +89,7 @@
"angstrom": (cm_per_ang, dimensions.length),
"Jy": (jansky_cgs, dimensions.specific_flux),
"counts": (1.0, dimensions.dimensionless),
+ "photons": (1.0, dimensions.dimensionless),
# for AstroPy compatibility
"solMass": (mass_sun_grams, dimensions.mass),
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list