[yt-svn] commit/yt: 5 new changesets
Bitbucket
commits-noreply at bitbucket.org
Thu Dec 13 10:29:11 PST 2012
5 new commits in yt:
https://bitbucket.org/yt_analysis/yt/changeset/ea797f1009ac/
changeset: ea797f1009ac
branch: yt
user: brittonsmith
date: 2012-12-12 21:48:20
summary: Adding new metallicity dependent x-ray emissivity and luminosity
fields created from tabulated Cloudy data.
affected #: 2 files
diff -r 04e4a6a7d10836f01864993c2507ffc6f89a00aa -r ea797f1009acf08ffda1b5b6e824fcacc66b1a50 yt/analysis_modules/spectral_integrator/api.py
--- a/yt/analysis_modules/spectral_integrator/api.py
+++ b/yt/analysis_modules/spectral_integrator/api.py
@@ -30,4 +30,7 @@
from .spectral_frequency_integrator import \
SpectralFrequencyIntegrator, \
- create_table_from_textfiles
+ create_table_from_textfiles, \
+ EmissivityIntegrator, \
+ add_xray_emissivity_field, \
+ add_xray_luminosity_field
diff -r 04e4a6a7d10836f01864993c2507ffc6f89a00aa -r ea797f1009acf08ffda1b5b6e824fcacc66b1a50 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
@@ -4,9 +4,11 @@
Author: Matthew Turk <matthewturk at gmail.com>
Affiliation: KIPAC/SLAC/Stanford
+Author: Britton Smith <brittons at origins.colorado.edu>
+Affiliation: Michigan State University
Homepage: http://yt-project.org/
License:
- Copyright (C) 2007-2011 Matthew Turk. All Rights Reserved.
+ Copyright (C) 2007-2012 Matthew Turk. All Rights Reserved.
This file is part of yt.
@@ -24,15 +26,17 @@
along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
+from exceptions import IOError
+import h5py
import numpy as np
+import os
from yt.funcs import *
from yt.data_objects.field_info_container import add_field
+from yt.utilities.exceptions import YTException
from yt.utilities.linear_interpolators import \
- UnilinearFieldInterpolator, \
- BilinearFieldInterpolator, \
- TrilinearFieldInterpolator
+ BilinearFieldInterpolator
class SpectralFrequencyIntegrator(object):
def __init__(self, table, field_names,
@@ -80,8 +84,8 @@
return 10**interp(dd)
add_field(name, function=frequency_bin_field,
projection_conversion="cm",
- units=r"\rm{ergs}\/\rm{cm}^{-3}\/\rm{s}^{-1}",
- projected_units=r"\rm{ergs}\/\rm{cm}^{-2}\/\rm{s}^{-1}")
+ units=r"\rm{ergs}\ \rm{cm}^{-3}\ \rm{s}^{-1}",
+ projected_units=r"\rm{ergs}\ \rm{cm}^{-2}\ \rm{s}^{-1}")
return name
def create_table_from_textfiles(pattern, rho_spec, e_spec, T_spec):
@@ -100,3 +104,212 @@
table[ri,:,ei] = [float(l.split()[-1]) for l in open(pattern%(i+1)) if l[0] != "#"]
return table
+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 EmissivityIntegrator(object):
+ r"""Class for making X-ray emissivity fields with hdf5 data tables
+ from Cloudy.
+ """
+ def __init__(self, filename=None):
+ r"""Initialize an EmissivityIntegrator object.
+
+ Keyword 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.
+
+ """
+ if filename is None:
+ filename = os.path.join(os.path.dirname(__file__),
+ "xray_emissivity.h5")
+ if not os.path.exists(filename):
+ raise IOError("File does not exist: %s." % 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"])
+ for field in ["emissivity_primordial", "emissivity_metals",
+ "log_nH", "log_T", "log_E"]:
+ 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 = 2.41799e17 * np.diff(self.E_bins)
+ del self.log_E
+
+ 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.
+
+ """
+ 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]))
+ 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])
+ # 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
+
+ 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)
+
+def add_xray_emissivity_field(e_min, e_max, filename=None,
+ with_metals=True,
+ constant_metallicity=None):
+ r"""Create an X-ray emissivity 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.
+
+ Keyword 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_Emissivity_{e_min}_{e_max}keV".
+ The units of the field are erg 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)
+
+ 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)
+
+ 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:
+ if with_metals:
+ my_Z = data["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
+
+ 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}",
+ projected_units=r"\rm{erg}\ \rm{cm}^{-2}\ \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.
+
+ Keyword 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.h.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)
+
+ 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
https://bitbucket.org/yt_analysis/yt/changeset/cdeb85f88fee/
changeset: cdeb85f88fee
branch: yt
user: brittonsmith
date: 2012-12-13 16:00:59
summary: Changing default filename for xray emissivity data to a directory
within the yt install.
affected #: 1 file
diff -r ea797f1009acf08ffda1b5b6e824fcacc66b1a50 -r cdeb85f88feeadd4f639f9e0f181f790212053be 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
@@ -131,8 +131,9 @@
"""
if filename is None:
- filename = os.path.join(os.path.dirname(__file__),
- "xray_emissivity.h5")
+ filename = os.path.join(os.environ["YT_DEST"],
+ "data", "xray_emissivity.h5")
+
if not os.path.exists(filename):
raise IOError("File does not exist: %s." % filename)
only_on_root(mylog.info, "Loading emissivity data from %s." % filename)
https://bitbucket.org/yt_analysis/yt/changeset/3d1d4b49a7ad/
changeset: 3d1d4b49a7ad
branch: yt
user: brittonsmith
date: 2012-12-13 17:48:56
summary: Adding photon emissivity field [photons s^{-1} cm^{-3}] and adding a
version check for the emissivity data. If the version does not match,
instructions are printed for downloading the latest data.
affected #: 2 files
diff -r cdeb85f88feeadd4f639f9e0f181f790212053be -r 3d1d4b49a7adedde4e21eb38acfff4ea4aeaf67d yt/analysis_modules/spectral_integrator/api.py
--- a/yt/analysis_modules/spectral_integrator/api.py
+++ b/yt/analysis_modules/spectral_integrator/api.py
@@ -33,4 +33,5 @@
create_table_from_textfiles, \
EmissivityIntegrator, \
add_xray_emissivity_field, \
- add_xray_luminosity_field
+ add_xray_luminosity_field, \
+ add_xray_photon_emissivity_field
diff -r cdeb85f88feeadd4f639f9e0f181f790212053be -r 3d1d4b49a7adedde4e21eb38acfff4ea4aeaf67d 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
@@ -38,6 +38,8 @@
from yt.utilities.linear_interpolators import \
BilinearFieldInterpolator
+xray_data_version = 1
+
class SpectralFrequencyIntegrator(object):
def __init__(self, table, field_names,
bounds, ev_bounds):
@@ -113,6 +115,11 @@
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 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")
+
class EmissivityIntegrator(object):
r"""Class for making X-ray emissivity fields with hdf5 data tables
from Cloudy.
@@ -130,9 +137,12 @@
Default: None.
"""
+
+ default_filename = False
if filename is None:
filename = os.path.join(os.environ["YT_DEST"],
"data", "xray_emissivity.h5")
+ default_filename = True
if not os.path.exists(filename):
raise IOError("File does not exist: %s." % filename)
@@ -140,6 +150,13 @@
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"]:
setattr(self, field, in_file[field][:])
@@ -151,7 +168,6 @@
[self.log_E[-1] - 0.5 * E_diff[-1],
self.log_E[-1] + 0.5 * E_diff[-1]]]))
self.dnu = 2.41799e17 * np.diff(self.E_bins)
- del self.log_E
def _get_interpolator(self, data, e_min, e_max):
r"""Create an interpolator for total emissivity in a
@@ -314,3 +330,80 @@
display_name=r"\rm{L}_{X}\ (%s-%s\ keV)" % (e_min, e_max),
units=r"\rm{erg}\ \rm{s}^{-1}")
return field_name
+
+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.
+
+ 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.
+
+ Keyword 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) * 1.60217646e-9
+
+ 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:
+ if with_metals:
+ my_Z = data["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
+
+ 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}",
+ projected_units=r"\rm{photons}\ \rm{cm}^{-2}\ \rm{s}^{-1}")
+ return field_name
https://bitbucket.org/yt_analysis/yt/changeset/8b3cf91894e3/
changeset: 8b3cf91894e3
branch: yt
user: brittonsmith
date: 2012-12-13 17:49:48
summary: Adding xray emissivity data to install script.
affected #: 1 file
diff -r 3d1d4b49a7adedde4e21eb38acfff4ea4aeaf67d -r 8b3cf91894e3b279ce3b78eb63daee265770c305 doc/install_script.sh
--- a/doc/install_script.sh
+++ b/doc/install_script.sh
@@ -399,6 +399,14 @@
( ${SHASUM} -c $1.sha512 2>&1 ) 1>> ${LOG_FILE} || do_exit
}
+function get_ytdata
+{
+ echo "Downloading $1 from yt-project.org"
+ [ -e $1 ] && return
+ ${GETFILE} "http://yt-project.org/data/$1" || do_exit
+ ( ${SHASUM} -c $1.sha512 2>&1 ) 1>> ${LOG_FILE} || do_exit
+}
+
ORIG_PWD=`pwd`
if [ -z "${DEST_DIR}" ]
@@ -407,6 +415,13 @@
exit 1
fi
+# Get supplemental data.
+
+mkdir -p ${DEST_DIR}/data
+cd ${DEST_DIR}/data
+echo 'de6d8c6ea849f0206d219303329a0276b3cce7c051eec34377d42aacbe0a4f47ac5145eb08966a338ecddd2b83c8f787ca9956508ad5c39ee2088ad875166410 xray_emissivity.h5' > xray_emissivity.h5.sha512
+get_ytdata xray_emissivity.h5
+
mkdir -p ${DEST_DIR}/src
cd ${DEST_DIR}/src
https://bitbucket.org/yt_analysis/yt/changeset/dbba2cffd0c3/
changeset: dbba2cffd0c3
branch: yt
user: jzuhone
date: 2012-12-13 19:29:10
summary: Merged in brittonsmith/yt (pull request #370)
affected #: 3 files
diff -r 7dff8c5bd37734e3b24c1c3e033b30d5cfb95f2d -r dbba2cffd0c3a6efe73e78109751618b82c01863 doc/install_script.sh
--- a/doc/install_script.sh
+++ b/doc/install_script.sh
@@ -399,6 +399,14 @@
( ${SHASUM} -c $1.sha512 2>&1 ) 1>> ${LOG_FILE} || do_exit
}
+function get_ytdata
+{
+ echo "Downloading $1 from yt-project.org"
+ [ -e $1 ] && return
+ ${GETFILE} "http://yt-project.org/data/$1" || do_exit
+ ( ${SHASUM} -c $1.sha512 2>&1 ) 1>> ${LOG_FILE} || do_exit
+}
+
ORIG_PWD=`pwd`
if [ -z "${DEST_DIR}" ]
@@ -407,6 +415,13 @@
exit 1
fi
+# Get supplemental data.
+
+mkdir -p ${DEST_DIR}/data
+cd ${DEST_DIR}/data
+echo 'de6d8c6ea849f0206d219303329a0276b3cce7c051eec34377d42aacbe0a4f47ac5145eb08966a338ecddd2b83c8f787ca9956508ad5c39ee2088ad875166410 xray_emissivity.h5' > xray_emissivity.h5.sha512
+get_ytdata xray_emissivity.h5
+
mkdir -p ${DEST_DIR}/src
cd ${DEST_DIR}/src
diff -r 7dff8c5bd37734e3b24c1c3e033b30d5cfb95f2d -r dbba2cffd0c3a6efe73e78109751618b82c01863 yt/analysis_modules/spectral_integrator/api.py
--- a/yt/analysis_modules/spectral_integrator/api.py
+++ b/yt/analysis_modules/spectral_integrator/api.py
@@ -30,4 +30,8 @@
from .spectral_frequency_integrator import \
SpectralFrequencyIntegrator, \
- create_table_from_textfiles
+ create_table_from_textfiles, \
+ EmissivityIntegrator, \
+ add_xray_emissivity_field, \
+ add_xray_luminosity_field, \
+ add_xray_photon_emissivity_field
diff -r 7dff8c5bd37734e3b24c1c3e033b30d5cfb95f2d -r dbba2cffd0c3a6efe73e78109751618b82c01863 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
@@ -4,9 +4,11 @@
Author: Matthew Turk <matthewturk at gmail.com>
Affiliation: KIPAC/SLAC/Stanford
+Author: Britton Smith <brittons at origins.colorado.edu>
+Affiliation: Michigan State University
Homepage: http://yt-project.org/
License:
- Copyright (C) 2007-2011 Matthew Turk. All Rights Reserved.
+ Copyright (C) 2007-2012 Matthew Turk. All Rights Reserved.
This file is part of yt.
@@ -24,16 +26,20 @@
along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
+from exceptions import IOError
+import h5py
import numpy as np
+import os
from yt.funcs import *
from yt.data_objects.field_info_container import add_field
+from yt.utilities.exceptions import YTException
from yt.utilities.linear_interpolators import \
- UnilinearFieldInterpolator, \
- BilinearFieldInterpolator, \
- TrilinearFieldInterpolator
+ BilinearFieldInterpolator
+xray_data_version = 1
+
class SpectralFrequencyIntegrator(object):
def __init__(self, table, field_names,
bounds, ev_bounds):
@@ -80,8 +86,8 @@
return 10**interp(dd)
add_field(name, function=frequency_bin_field,
projection_conversion="cm",
- units=r"\rm{ergs}\/\rm{cm}^{-3}\/\rm{s}^{-1}",
- projected_units=r"\rm{ergs}\/\rm{cm}^{-2}\/\rm{s}^{-1}")
+ units=r"\rm{ergs}\ \rm{cm}^{-3}\ \rm{s}^{-1}",
+ projected_units=r"\rm{ergs}\ \rm{cm}^{-2}\ \rm{s}^{-1}")
return name
def create_table_from_textfiles(pattern, rho_spec, e_spec, T_spec):
@@ -100,3 +106,304 @@
table[ri,:,ei] = [float(l.split()[-1]) for l in open(pattern%(i+1)) if l[0] != "#"]
return table
+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 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")
+
+class EmissivityIntegrator(object):
+ r"""Class for making X-ray emissivity fields with hdf5 data tables
+ from Cloudy.
+ """
+ def __init__(self, filename=None):
+ r"""Initialize an EmissivityIntegrator object.
+
+ Keyword 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.
+
+ """
+
+ default_filename = False
+ if filename is None:
+ filename = os.path.join(os.environ["YT_DEST"],
+ "data", "xray_emissivity.h5")
+ default_filename = True
+
+ if not os.path.exists(filename):
+ raise IOError("File does not exist: %s." % 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"]:
+ 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 = 2.41799e17 * np.diff(self.E_bins)
+
+ 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.
+
+ """
+ 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]))
+ 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])
+ # 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
+
+ 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)
+
+def add_xray_emissivity_field(e_min, e_max, filename=None,
+ with_metals=True,
+ constant_metallicity=None):
+ r"""Create an X-ray emissivity 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.
+
+ Keyword 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_Emissivity_{e_min}_{e_max}keV".
+ The units of the field are erg 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)
+
+ 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)
+
+ 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:
+ if with_metals:
+ my_Z = data["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
+
+ 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}",
+ projected_units=r"\rm{erg}\ \rm{cm}^{-2}\ \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.
+
+ Keyword 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.h.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)
+
+ 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
+
+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.
+
+ 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.
+
+ Keyword 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) * 1.60217646e-9
+
+ 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:
+ if with_metals:
+ my_Z = data["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
+
+ 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}",
+ projected_units=r"\rm{photons}\ \rm{cm}^{-2}\ \rm{s}^{-1}")
+ return field_name
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