[yt-svn] commit/yt: ngoldbaum: Merged in jzuhone/yt (pull request #2465)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Thu Jan 19 09:10:56 PST 2017
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/f0fa273ba464/
Changeset: f0fa273ba464
Branch: yt
User: ngoldbaum
Date: 2017-01-19 17:10:29+00:00
Summary: Merged in jzuhone/yt (pull request #2465)
Port the spectral_integrator analysis module into yt.fields
Affected #: 11 files
diff -r f00742795e0970ba78ded9ea01644a6a40821b13 -r f0fa273ba464df9040cf988857236495897424bb doc/source/analyzing/XrayEmissionFields.ipynb
--- /dev/null
+++ b/doc/source/analyzing/XrayEmissionFields.ipynb
@@ -0,0 +1,220 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "> Note: If you came here trying to figure out how to create simulated X-ray photons and observations,\n",
+ " you should go [here](analysis_modules/photon_simulator.html) instead."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This functionality provides the ability to create metallicity-dependent X-ray luminosity, emissivity, and photon emissivity fields for a given photon energy range. This works by interpolating from emission tables created from the photoionization code [Cloudy](http://nublado.org/) or the collisional ionization database [AtomDB](http://www.atomdb.org). These can be downloaded from http://yt-project.org/data from the command line like so:\n",
+ "\n",
+ "`# Put the data in a directory you specify` \n",
+ "`yt download cloudy_emissivity_v2.h5 /path/to/data`\n",
+ "\n",
+ "`# Put the data in the location set by \"supp_data_dir\"` \n",
+ "`yt download apec_emissivity_v2.h5 supp_data_dir`\n",
+ "\n",
+ "The data path can be a directory on disk, or it can be \"supp_data_dir\", which will download the data to the directory specified by the `\"supp_data_dir\"` yt configuration entry. It is easiest to put these files in the directory from which you will be running yt or `\"supp_data_dir\"`, but see the note below about putting them in alternate locations."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Emission fields can be made for any energy interval between 0.1 keV and 100 keV, and will always be created for luminosity $(\\rm{erg~s^{-1}})$, emissivity $\\rm{(erg~s^{-1}~cm^{-3})}$, and photon emissivity $\\rm{(photons~s^{-1}~cm^{-3})}$. The only required arguments are the\n",
+ "dataset object, and the minimum and maximum energies of the energy band. However, typically one needs to decide what will be used for the metallicity. This can either be a floating-point value representing a spatially constant metallicity, or a prescription for a metallicity field, e.g. `(\"gas\", \"metallicity\")`. For this first example, where the dataset has no metallicity field, we'll just assume $Z = 0.3~Z_\\odot$ everywhere:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "import yt\n",
+ "\n",
+ "ds = yt.load(\"GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150\")\n",
+ "\n",
+ "xray_fields = yt.add_xray_emissivity_field(ds, 0.5, 7.0, table_type='apec', metallicity=0.3)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "> Note: If you place the HDF5 emissivity tables in a location other than the current working directory or the location \n",
+ " specified by the \"supp_data_dir\" configuration value, you will need to specify it in the call to \n",
+ " `add_xray_emissivity_field`: \n",
+ " `xray_fields = yt.add_xray_emissivity_field(ds, 0.5, 7.0, data_dir=\"/path/to/data\", table_type='apec', metallicity=0.3)`"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Having made the fields, one can see which fields were made:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "print (xray_fields)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The luminosity field is useful for summing up in regions like this:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "sp = ds.sphere(\"c\", (2.0, \"Mpc\"))\n",
+ "print (sp.quantities.total_quantity(\"xray_luminosity_0.5_7.0_keV\"))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Whereas the emissivity fields may be useful in derived fields or for plotting:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "slc = yt.SlicePlot(ds, 'z', ['xray_emissivity_0.5_7.0_keV','xray_photon_emissivity_0.5_7.0_keV'],\n",
+ " width=(0.75, \"Mpc\"))\n",
+ "slc.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The emissivity and the luminosity fields take the values one would see in the frame of the source. However, if one wishes to make projections of the X-ray emission from a cosmologically distant object, the energy band will be redshifted. For this case, one can supply a `redshift` parameter and a `Cosmology` object (either from the dataset or one made on your own) to compute X-ray intensity fields along with the emissivity and luminosity fields.\n",
+ "\n",
+ "This example shows how to do that, Where we also use a spatially dependent metallicity field and the Cloudy tables instead of the APEC tables we used previously:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "ds2 = yt.load(\"D9p_500/10MpcBox_HartGal_csf_a0.500.d\")\n",
+ "\n",
+ "# In this case, use the redshift and cosmology from the dataset, \n",
+ "# but in theory you could put in something different\n",
+ "xray_fields2 = yt.add_xray_emissivity_field(ds2, 0.5, 2.0, redshift=ds2.current_redshift, cosmology=ds2.cosmology,\n",
+ " metallicity=(\"gas\", \"metallicity\"), table_type='cloudy')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now, one can see that two new fields have been added, corresponding to X-ray intensity / surface brightness when projected:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "print (xray_fields2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Note also that the energy range now corresponds to the *observer* frame, whereas in the source frame the energy range is between `emin*(1+redshift)` and `emax*(1+redshift)`. Let's zoom in on a galaxy and make a projection of the intensity fields:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "prj = yt.ProjectionPlot(ds2, \"x\", [\"xray_intensity_0.5_2.0_keV\", \"xray_photon_intensity_0.5_2.0_keV\"],\n",
+ " center=\"max\", width=(40, \"kpc\"))\n",
+ "prj.set_zlim(\"xray_intensity_0.5_2.0_keV\", 1.0e-32, 5.0e-24)\n",
+ "prj.set_zlim(\"xray_photon_intensity_0.5_2.0_keV\", 1.0e-24, 5.0e-16)\n",
+ "prj.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "> Warning: The X-ray fields depend on the number density of hydrogen atoms, given by the yt field\n",
+ " `H_nuclei_density`. In the case of the APEC model, this assumes that all of the hydrogen in your\n",
+ " dataset is ionized, whereas in the Cloudy model the ionization level is taken into account. If \n",
+ " this field is not defined (either in the dataset or by the user), it will be constructed using\n",
+ " abundance information from your dataset. Finally, if your dataset contains no abundance information,\n",
+ " a primordial hydrogen mass fraction (X = 0.76) will be assumed."
+ ]
+ }
+ ],
+ "metadata": {
+ "anaconda-cloud": {},
+ "kernelspec": {
+ "display_name": "Python [default]",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.5.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff -r f00742795e0970ba78ded9ea01644a6a40821b13 -r f0fa273ba464df9040cf988857236495897424bb 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
@@ -3,6 +3,14 @@
Constructing Mock X-ray Observations
------------------------------------
+.. warning::
+
+ The ``photon_simulator`` analysis module has been deprecated; it is
+ no longer being updated, and it will be removed in a future version
+ of yt. Users are encouraged to download and use the
+ `pyXSIM <http://hea-www.cfa.harvard.edu/~jzuhone/pyxsim>`_ package
+ instead.
+
.. note::
If you just want to create derived fields for X-ray emission,
diff -r f00742795e0970ba78ded9ea01644a6a40821b13 -r f0fa273ba464df9040cf988857236495897424bb doc/source/analyzing/analysis_modules/xray_emission_fields.rst
--- a/doc/source/analyzing/analysis_modules/xray_emission_fields.rst
+++ /dev/null
@@ -1,78 +0,0 @@
-.. _xray_emission_fields:
-
-X-ray Emission Fields
-=====================
-.. 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 photon emissivity fields for a given
-photon energy range. This works by interpolating from emission tables
-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 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
-
- import yt
- from yt.analysis_modules.spectral_integrator.api import \
- add_xray_emissivity_field
-
- xray_fields = add_xray_emissivity_field(ds, 0.5, 7.0)
-
-Additional keyword arguments are:
-
- * **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
- 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. It should be given in unit of solar metallicity.
- Default: None.
-
-The resulting fields can be used like all normal fields. The function will return the names of
-the created fields in a Python list.
-
-.. code-block:: python
-
- import yt
- from yt.analysis_modules.spectral_integrator.api import \
- add_xray_emissivity_field
-
- xray_fields = add_xray_emissivity_field(ds, 0.5, 7.0, filename="apec_emissivity.h5")
-
- ds = yt.load("enzo_tiny_cosmology/DD0046/DD0046")
- plot = yt.SlicePlot(ds, 'x', 'xray_luminosity_0.5_7.0_keV')
- plot.save()
- plot = yt.ProjectionPlot(ds, 'x', 'xray_emissivity_0.5_7.0_keV')
- plot.save()
- 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.
diff -r f00742795e0970ba78ded9ea01644a6a40821b13 -r f0fa273ba464df9040cf988857236495897424bb doc/source/analyzing/xray_emission_fields.rst
--- /dev/null
+++ b/doc/source/analyzing/xray_emission_fields.rst
@@ -0,0 +1,3 @@
+.. _xray_emission_fields:
+
+.. notebook:: XrayEmissionFields.ipynb
diff -r f00742795e0970ba78ded9ea01644a6a40821b13 -r f0fa273ba464df9040cf988857236495897424bb tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -53,7 +53,7 @@
local_tipsy_002:
- yt/frontends/tipsy/tests/test_outputs.py
- local_varia_007:
+ local_varia_008:
- yt/analysis_modules/radmc3d_export
- yt/frontends/moab/tests/test_c5.py
- yt/analysis_modules/photon_simulator/tests/test_spectra.py
@@ -61,6 +61,7 @@
- yt/visualization/volume_rendering/tests/test_vr_orientation.py
- yt/visualization/volume_rendering/tests/test_mesh_render.py
- yt/visualization/tests/test_mesh_slices.py:test_tri2
+ - yt/fields/tests/test_xray_fields.py
local_orion_001:
- yt/frontends/boxlib/tests/test_orion.py
diff -r f00742795e0970ba78ded9ea01644a6a40821b13 -r f0fa273ba464df9040cf988857236495897424bb yt/__init__.py
--- a/yt/__init__.py
+++ b/yt/__init__.py
@@ -125,7 +125,8 @@
ValidateSpatial, \
ValidateGridType, \
add_field, \
- derived_field
+ derived_field, \
+ add_xray_emissivity_field
from yt.data_objects.api import \
DatasetSeries, ImageArray, \
diff -r f00742795e0970ba78ded9ea01644a6a40821b13 -r f0fa273ba464df9040cf988857236495897424bb yt/analysis_modules/spectral_integrator/api.py
--- a/yt/analysis_modules/spectral_integrator/api.py
+++ b/yt/analysis_modules/spectral_integrator/api.py
@@ -1,18 +1,8 @@
-"""
-API for spectral_integrator
-
-
-
-"""
+from yt.funcs import issue_deprecation_warning
-#-----------------------------------------------------------------------------
-# Copyright (c) 2013, yt Development Team.
-#
-# Distributed under the terms of the Modified BSD License.
-#
-# The full license is in the file COPYING.txt, distributed with this software.
-#-----------------------------------------------------------------------------
+issue_deprecation_warning("The spectral_integrator module is deprecated. "
+ "'add_xray_emissivity_field' can now be imported "
+ "from the yt module.")
-from .spectral_frequency_integrator import \
- EmissivityIntegrator, \
+from yt.fields.xray_emission_fields import \
add_xray_emissivity_field
\ No newline at end of file
diff -r f00742795e0970ba78ded9ea01644a6a40821b13 -r f0fa273ba464df9040cf988857236495897424bb yt/analysis_modules/spectral_integrator/spectral_frequency_integrator.py
--- a/yt/analysis_modules/spectral_integrator/spectral_frequency_integrator.py
+++ /dev/null
@@ -1,279 +0,0 @@
-"""
-Integrator classes to deal with interpolation and integration of input spectral
-bins. Currently only supports Cloudy and APEC-style data.
-
-
-
-"""
-
-#-----------------------------------------------------------------------------
-# Copyright (c) 2013, yt Development Team.
-#
-# Distributed under the terms of the Modified BSD License.
-#
-# The full license is in the file COPYING.txt, distributed with this software.
-#-----------------------------------------------------------------------------
-
-from yt.utilities.on_demand_imports import _h5py as h5py
-import numpy as np
-import os
-
-from yt.funcs import \
- download_file, \
- mylog, \
- only_on_root
-
-from yt.utilities.exceptions import YTFieldNotFound
-from yt.utilities.exceptions import YTException
-from yt.utilities.linear_interpolators import \
- UnilinearFieldInterpolator, BilinearFieldInterpolator
-from yt.utilities.physical_constants import \
- 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=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")):
- data_dir = os.path.join(os.environ["YT_DEST"], "data")
- else:
- data_dir = "."
- data_path = os.path.join(data_dir, data_file)
- if not os.path.exists(data_path):
- mylog.info("Attempting to download supplementary data from %s to %s." %
- (data_url, data_dir))
- fn = download_file(os.path.join(data_url, data_file), data_path)
- if fn != data_path:
- raise RuntimeError("Failed to download supplementary data.")
- return data_path
-
-class EnergyBoundsException(YTException):
- def __init__(self, lower, upper):
- self.lower = lower
- self.upper = upper
-
- def __str__(self):
- return "Energy bounds are %e to %e keV." % \
- (self.lower, self.upper)
-
-class ObsoleteDataException(YTException):
- def __str__(self):
- 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
- from Cloudy.
-
- Initialize an EmissivityIntegrator object.
-
- Parameters
- ----------
- filename: string, default None
- 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.
-
- """
- def __init__(self, filename=None):
-
- default_filename = False
- if filename is None:
- filename = _get_data_file()
- default_filename = True
-
- if not os.path.exists(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:
- only_on_root(mylog.info, in_file.attrs["info"])
- if default_filename and \
- in_file.attrs["version"] < xray_data_version:
- raise ObsoleteDataException()
- else:
- only_on_root(mylog.info, "X-ray emissivity data version: %s." % \
- in_file.attrs["version"])
-
- for field in ["emissivity_primordial", "emissivity_metals",
- "log_nH", "log_T", "log_E"]:
- if field in in_file:
- setattr(self, field, in_file[field][:])
- in_file.close()
-
- E_diff = np.diff(self.log_E)
- 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):
- 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(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 = 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])/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)
- 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)
-
- return emiss
-
-def add_xray_emissivity_field(ds, e_min, e_max,
- filename=None,
- with_metals=True,
- constant_metallicity=None):
- r"""Create X-ray emissivity fields 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.
- filename: string, optional
- 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, 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, optional
- If specified, assume a constant metallicity for the emission
- from metals. The *with_metals* keyword must be set to False
- to use this. It should be given in unit of solar metallicity.
- Default: None.
-
- 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 *
- >>> 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_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)
-
- 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["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["gas","metallicity"]
- elif constant_metallicity is not None:
- my_Z = constant_metallicity
- my_emissivity += my_Z * np.power(10, em_Z(dd))
-
- return data["gas","H_number_density"]**2 * \
- YTArray(my_emissivity, "erg*cm**3/s")
-
- emiss_name = "xray_emissivity_%s_%s_keV" % (e_min, e_max)
- ds.add_field(("gas", 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[emiss_name] * data["cell_volume"]
-
- lum_name = "xray_luminosity_%s_%s_keV" % (e_min, e_max)
- ds.add_field(("gas", lum_name), function=_luminosity_field,
- display_name=r"\rm{L}_{X} (%s-%s keV)" % (e_min, e_max),
- units="erg/s")
-
- def _photon_emissivity_field(field, data):
- dd = {"log_nH" : np.log10(data["gas","H_number_density"]),
- "log_T" : np.log10(data["gas","temperature"])}
-
- my_emissivity = np.power(10, emp_0(dd))
- if emp_Z is not None:
- if with_metals:
- my_Z = data["gas","metallicity"]
- elif constant_metallicity is not None:
- my_Z = constant_metallicity
- my_emissivity += my_Z * np.power(10, emp_Z(dd))
-
- return data["gas","H_number_density"]**2 * \
- YTArray(my_emissivity, "photons*cm**3/s")
-
- phot_name = "xray_photon_emissivity_%s_%s_keV" % (e_min, e_max)
- ds.add_field(("gas", 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
diff -r f00742795e0970ba78ded9ea01644a6a40821b13 -r f0fa273ba464df9040cf988857236495897424bb yt/fields/api.py
--- a/yt/fields/api.py
+++ b/yt/fields/api.py
@@ -43,3 +43,5 @@
FieldDetector
from .field_info_container import \
FieldInfoContainer
+from .xray_emission_fields import \
+ add_xray_emissivity_field
\ No newline at end of file
diff -r f00742795e0970ba78ded9ea01644a6a40821b13 -r f0fa273ba464df9040cf988857236495897424bb yt/fields/tests/test_xray_fields.py
--- /dev/null
+++ b/yt/fields/tests/test_xray_fields.py
@@ -0,0 +1,40 @@
+from yt.fields.xray_emission_fields import \
+ add_xray_emissivity_field
+from yt.utilities.answer_testing.framework import \
+ requires_ds, can_run_ds, data_dir_load, \
+ ProjectionValuesTest, FieldValuesTest
+
+def setup():
+ from yt.config import ytcfg
+ ytcfg["yt","__withintesting"] = "True"
+
+def check_xray_fields(ds_fn, fields):
+ if not can_run_ds(ds_fn): return
+ dso = [ None, ("sphere", ("m", (0.1, 'unitary')))]
+ for field in fields:
+ for axis in [0, 1, 2]:
+ for dobj_name in dso:
+ yield ProjectionValuesTest(ds_fn, axis, field,
+ None, dobj_name)
+ yield FieldValuesTest(ds_fn, field, dobj_name)
+
+sloshing = "GasSloshingLowRes/sloshing_low_res_hdf5_plt_cnt_0300"
+ at requires_ds(sloshing, big_data=True)
+def test_sloshing_apec():
+ ds = data_dir_load(sloshing)
+ fields = add_xray_emissivity_field(ds, 0.5, 7.0, table_type="apec",
+ metallicity=0.3)
+ for test in check_xray_fields(ds, fields):
+ test_sloshing_apec.__name__ = test.description
+ yield test
+
+d9p = "D9p_500/10MpcBox_HartGal_csf_a0.500.d"
+ at requires_ds(d9p, big_data=True)
+def test_d9p_cloudy():
+ ds = data_dir_load(d9p)
+ fields = add_xray_emissivity_field(ds, 0.5, 2.0, redshift=ds.current_redshift,
+ table_type="cloudy", cosmology=ds.cosmology,
+ metallicity=("gas", "metallicity"))
+ for test in check_xray_fields(ds, fields):
+ test_d9p_cloudy.__name__ = test.description
+ yield test
diff -r f00742795e0970ba78ded9ea01644a6a40821b13 -r f0fa273ba464df9040cf988857236495897424bb yt/fields/xray_emission_fields.py
--- /dev/null
+++ b/yt/fields/xray_emission_fields.py
@@ -0,0 +1,304 @@
+"""
+Integrator classes to deal with interpolation and integration of input spectral
+bins. Currently only supports Cloudy and APEC-style data.
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+from yt.utilities.on_demand_imports import _h5py as h5py
+import numpy as np
+import os
+
+from yt.config import ytcfg
+from yt.fields.derived_field import DerivedField
+from yt.funcs import mylog, only_on_root, issue_deprecation_warning
+from yt.utilities.exceptions import YTFieldNotFound
+from yt.utilities.exceptions import YTException
+from yt.utilities.linear_interpolators import \
+ UnilinearFieldInterpolator, BilinearFieldInterpolator
+from yt.units.yt_array import YTArray, YTQuantity
+from yt.utilities.cosmology import Cosmology
+
+data_version = {"cloudy": 2,
+ "apec": 2}
+
+data_url = "http://yt-project.org/data"
+
+def _get_data_file(table_type, data_dir=None):
+ data_file = "%s_emissivity_v%d.h5" % (table_type, data_version[table_type])
+ if data_dir is None:
+ supp_data_dir = ytcfg.get("yt", "supp_data_dir")
+ data_dir = supp_data_dir if os.path.exists(supp_data_dir) else "."
+ data_path = os.path.join(data_dir, data_file)
+ if not os.path.exists(data_path):
+ msg = "Failed to find emissivity data file %s! " % data_file + \
+ "Please download from http://yt-project.org/data!"
+ mylog.error(msg)
+ raise IOError(msg)
+ return data_path
+
+class EnergyBoundsException(YTException):
+ def __init__(self, lower, upper):
+ self.lower = lower
+ self.upper = upper
+
+ def __str__(self):
+ return "Energy bounds are %e to %e keV." % \
+ (self.lower, self.upper)
+
+class ObsoleteDataException(YTException):
+ def __init__(self, table_type):
+ data_file = "%s_emissivity_v%d.h5" % (table_type, data_version[table_type])
+ self.msg = "X-ray emissivity data is out of date.\n"
+ self.msg += "Download the latest data from %s/%s." % (data_url, data_file)
+
+ def __str__(self):
+ return self.msg
+
+class XrayEmissivityIntegrator(object):
+ r"""Class for making X-ray emissivity fields. Uses hdf5 data tables
+ generated from Cloudy and AtomDB/APEC.
+
+ Initialize an XrayEmissivityIntegrator object.
+
+ Parameters
+ ----------
+ table_type: string
+ The type of data to use when computing the emissivity values. If "cloudy",
+ a file called "cloudy_emissivity.h5" is used, for photoionized
+ plasmas. If, "apec", a file called "apec_emissivity.h5" is used for
+ collisionally ionized plasmas. These files contain emissivity tables
+ for primordial elements and for metals at solar metallicity for the
+ energy range 0.1 to 100 keV.
+ redshift : float, optional
+ The cosmological redshift of the source of the field. Default: 0.0.
+ data_dir : string, optional
+ The location to look for the data table in. If not supplied, the file
+ will be looked for in the location of the YT_DEST environment variable
+ or in the current working directory.
+ use_metals : boolean, optional
+ If set to True, the emissivity will include contributions from metals.
+ Default: True
+ """
+ def __init__(self, table_type, redshift=0.0, data_dir=None, use_metals=True):
+
+ filename = _get_data_file(table_type, data_dir=data_dir)
+ only_on_root(mylog.info, "Loading emissivity data from %s." % filename)
+ in_file = h5py.File(filename, "r")
+ if "info" in in_file.attrs:
+ only_on_root(mylog.info, in_file.attrs["info"].decode('utf8'))
+ if in_file.attrs["version"] != data_version[table_type]:
+ raise ObsoleteDataException(table_type)
+ else:
+ only_on_root(mylog.info, "X-ray '%s' emissivity data version: %s." % \
+ (table_type, in_file.attrs["version"]))
+
+ self.log_T = in_file["log_T"][:]
+ self.emissivity_primordial = in_file["emissivity_primordial"][:]
+ if "log_nH" in in_file:
+ self.log_nH = in_file["log_nH"][:]
+ if use_metals:
+ self.emissivity_metals = in_file["emissivity_metals"][:]
+ self.ebin = YTArray(in_file["E"], "keV")
+ in_file.close()
+ self.dE = np.diff(self.ebin)
+ self.emid = 0.5*(self.ebin[1:]+self.ebin[:-1]).to("erg")
+ self.redshift = redshift
+
+ def get_interpolator(self, data_type, e_min, e_max, energy=True):
+ data = getattr(self, "emissivity_%s" % data_type)
+ if not energy:
+ data = data[..., :] / self.emid.v
+ e_min = YTQuantity(e_min, "keV")*(1.0+self.redshift)
+ e_max = YTQuantity(e_max, "keV")*(1.0+self.redshift)
+ if (e_min - self.ebin[0]) / e_min < -1e-3 or \
+ (e_max - self.ebin[-1]) / e_max > 1e-3:
+ raise EnergyBoundsException(self.ebin[0], self.ebin[-1])
+ e_is, e_ie = np.digitize([e_min, e_max], self.ebin)
+ e_is = np.clip(e_is - 1, 0, self.ebin.size - 1)
+ e_ie = np.clip(e_ie, 0, self.ebin.size - 1)
+
+ my_dE = self.dE[e_is: e_ie].copy()
+ # clip edge bins if the requested range is smaller
+ my_dE[0] -= e_min - self.ebin[e_is]
+ my_dE[-1] -= self.ebin[e_ie] - e_max
+
+ interp_data = (data[..., e_is:e_ie]*my_dE).sum(axis=-1)
+ if data.ndim == 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)
+
+ return emiss
+
+def add_xray_emissivity_field(ds, e_min, e_max, redshift=0.0,
+ metallicity=("gas", "metallicity"),
+ table_type="cloudy", data_dir=None,
+ cosmology=None, **kwargs):
+ r"""Create X-ray emissivity fields 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.
+ redshift : float, optional
+ The cosmological redshift of the source of the field. Default: 0.0.
+ metallicity : field or float, optional
+ Either the name of a metallicity field or a single floating-point
+ number specifying a spatially constant metallicity. Must be in
+ solar units. If set to None, no metals will be assumed. Default:
+ ("gas", "metallicity")
+ table_type : string, optional
+ The type of emissivity table to be used when creating the fields.
+ Options are "cloudy" or "apec". Default: "cloudy"
+ data_dir : string, optional
+ The location to look for the data table in. If not supplied, the file
+ will be looked for in the location of the YT_DEST environment variable
+ or in the current working directory.
+ cosmology : :class:`~yt.utilities.cosmology.Cosmology`, optional
+ If set and redshift > 0.0, this cosmology will be used when computing the
+ cosmological dependence of the emission fields. If not set, yt's default
+ LCDM cosmology will be used.
+
+ 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
+ --------
+
+ >>> import yt
+ >>> ds = yt.load("sloshing_nomag2_hdf5_plt_cnt_0100")
+ >>> yt.add_xray_emissivity_field(ds, 0.5, 2)
+ >>> p = yt.ProjectionPlot(ds, 'x', "xray_emissivity_0.5_2_keV")
+ >>> p.save()
+ """
+ # The next several if constructs are for backwards-compatibility
+ if "constant_metallicity" in kwargs:
+ issue_deprecation_warning("The \"constant_metallicity\" parameter is deprecated. Set "
+ "the \"metallicity\" parameter to a constant float value instead.")
+ metallicity = kwargs["constant_metallicity"]
+ if "with_metals" in kwargs:
+ issue_deprecation_warning("The \"with_metals\" parameter is deprecated. Use the "
+ "\"metallicity\" parameter to choose a constant or "
+ "spatially varying metallicity.")
+ if kwargs["with_metals"] and isinstance(metallicity, float):
+ raise RuntimeError("\"with_metals=True\", but you specified a constant metallicity!")
+ if not kwargs["with_metals"] and not isinstance(metallicity, float):
+ raise RuntimeError("\"with_metals=False\", but you didn't specify a constant metallicity!")
+ if not isinstance(metallicity, float) and metallicity is not None:
+ try:
+ metallicity = ds._get_field_info(*metallicity)
+ except YTFieldNotFound:
+ raise RuntimeError("Your dataset does not have a {} field! ".format(metallicity) +
+ "Perhaps you should specify a constant metallicity instead?")
+
+ my_si = XrayEmissivityIntegrator(table_type, data_dir=data_dir, redshift=redshift)
+
+ em_0 = my_si.get_interpolator("primordial", e_min, e_max)
+ emp_0 = my_si.get_interpolator("primordial", e_min, e_max, energy=False)
+ if metallicity is not None:
+ em_Z = my_si.get_interpolator("metals", e_min, e_max)
+ emp_Z = my_si.get_interpolator("metals", e_min, e_max, energy=False)
+
+ def _emissivity_field(field, data):
+ dd = {"log_nH": np.log10(data["gas", "H_nuclei_density"]),
+ "log_T": np.log10(data["gas", "temperature"])}
+
+ my_emissivity = np.power(10, em_0(dd))
+ if metallicity is not None:
+ if isinstance(metallicity, DerivedField):
+ my_Z = data[metallicity.name]
+ else:
+ my_Z = metallicity
+ my_emissivity += my_Z * np.power(10, em_Z(dd))
+
+ return data["gas","H_nuclei_density"]**2 * \
+ YTArray(my_emissivity, "erg*cm**3/s")
+
+ emiss_name = "xray_emissivity_%s_%s_keV" % (e_min, e_max)
+ ds.add_field(("gas", emiss_name), function=_emissivity_field,
+ display_name=r"\epsilon_{X} (%s-%s keV)" % (e_min, e_max),
+ sampling_type="cell", units="erg/cm**3/s")
+
+ def _luminosity_field(field, data):
+ return data[emiss_name] * data["cell_volume"]
+
+ lum_name = "xray_luminosity_%s_%s_keV" % (e_min, e_max)
+ ds.add_field(("gas", lum_name), function=_luminosity_field,
+ display_name=r"\rm{L}_{X} (%s-%s keV)" % (e_min, e_max),
+ sampling_type="cell", units="erg/s")
+
+ def _photon_emissivity_field(field, data):
+ dd = {"log_nH": np.log10(data["gas", "H_nuclei_density"]),
+ "log_T": np.log10(data["gas", "temperature"])}
+
+ my_emissivity = np.power(10, emp_0(dd))
+ if metallicity is not None:
+ if isinstance(metallicity, DerivedField):
+ my_Z = data[metallicity.name]
+ else:
+ my_Z = metallicity
+ my_emissivity += my_Z * np.power(10, emp_Z(dd))
+
+ return data["gas", "H_nuclei_density"]**2 * \
+ YTArray(my_emissivity, "photons*cm**3/s")
+
+ phot_name = "xray_photon_emissivity_%s_%s_keV" % (e_min, e_max)
+ ds.add_field(("gas", phot_name), function=_photon_emissivity_field,
+ display_name=r"\epsilon_{X} (%s-%s keV)" % (e_min, e_max),
+ sampling_type="cell", units="photons/cm**3/s")
+
+ fields = [emiss_name, lum_name, phot_name]
+
+ if redshift > 0.0:
+
+ if cosmology is None:
+ if hasattr(ds, "cosmology"):
+ cosmology = ds.cosmology
+ else:
+ cosmology = Cosmology()
+
+ D_L = cosmology.luminosity_distance(0.0, redshift)
+ angular_scale = 1.0/cosmology.angular_scale(0.0, redshift)
+ dist_fac = 1.0/(4.0*np.pi*D_L*D_L*angular_scale*angular_scale)
+
+ ei_name = "xray_intensity_%s_%s_keV" % (e_min, e_max)
+ def _intensity_field(field, data):
+ I = dist_fac*data[emiss_name]
+ return I.in_units("erg/cm**3/s/arcsec**2")
+ ds.add_field(("gas", ei_name), function=_intensity_field,
+ display_name=r"I_{X} (%s-%s keV)" % (e_min, e_max),
+ sampling_type="cell", units="erg/cm**3/s/arcsec**2")
+
+ i_name = "xray_photon_intensity_%s_%s_keV" % (e_min, e_max)
+ def _photon_intensity_field(field, data):
+ I = (1.0+redshift)*dist_fac*data[phot_name]
+ return I.in_units("photons/cm**3/s/arcsec**2")
+ ds.add_field(("gas", i_name), function=_photon_intensity_field,
+ display_name=r"I_{X} (%s-%s keV)" % (e_min, e_max),
+ sampling_type="cell", units="photons/cm**3/s/arcsec**2")
+
+ fields += [ei_name, i_name]
+
+ [mylog.info("Adding %s field." % field) for field in fields]
+
+ return fields
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