[yt-svn] commit/yt: chummels: Merged in jzuhone/yt (pull request #1616)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Thu Jul 9 16:38:13 PDT 2015
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/98c8bdfbc1e0/
Changeset: 98c8bdfbc1e0
Branch: yt
User: chummels
Date: 2015-07-09 23:38:02+00:00
Summary: Merged in jzuhone/yt (pull request #1616)
Speedups and improvements for photon_simulator
Affected #: 4 files
diff -r c66deb7d69cab3114e9563d6cf9ed57521675ff4 -r 98c8bdfbc1e00d4ee2388663a7f4072bbf98fc57 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
@@ -10,6 +10,10 @@
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
simulators to produce realistic observations or can be analyzed in-line.
+
+For detailed information about the design of the algorithm in yt, check
+out `the SciPy 2014 Proceedings. <http://conference.scipy.org/proceedings/scipy2014/zuhone.html>`_.
+
The algorithm is based off of that implemented in
`PHOX <http://www.mpa-garching.mpg.de/~kdolag/Phox/>`_ for SPH datasets
by Veronica Biffi and Klaus Dolag. There are two relevant papers:
@@ -139,6 +143,12 @@
the optional keyword ``thermal_broad`` is set to ``True``, the spectral
lines will be thermally broadened.
+.. note::
+
+ ``SpectralModel`` objects based on XSPEC models (both the thermal
+ emission and Galactic absorption models mentioned below) only work
+ in Python 2.7, since currently PyXspec only works with Python 2.x.
+
Now that we have our ``SpectralModel`` that gives us a spectrum, we need
to connect this model to a ``PhotonModel`` that will connect the field
data in the ``data_source`` to the spectral model to actually generate
@@ -148,7 +158,8 @@
.. code:: python
thermal_model = ThermalPhotonModel(apec_model, X_H=0.75, Zmet=0.3,
- photons_per_chunk=100000000)
+ photons_per_chunk=100000000,
+ method="invert_cdf")
Where we pass in the ``SpectralModel``, and can optionally set values for
the hydrogen mass fraction ``X_H`` and metallicity ``Z_met``. If
@@ -165,6 +176,18 @@
this parameter needs to be set higher, or if you are looking to decrease memory
usage, you might set this parameter lower.
+The ``method`` keyword argument is also optional, and determines how the individual
+photon energies are generated from the spectrum. It may be set to one of two values:
+
+* ``method="invert_cdf"``: Construct the cumulative distribution function of the spectrum and invert
+ it, using uniformly drawn random numbers to determine the photon energies (fast, but relies
+ on construction of the CDF and interpolation between the points, so for some spectra it
+ may not be accurate enough).
+* ``method="accept_reject"``: Generate the photon energies from the spectrum using an acceptance-rejection
+ technique (accurate, but likely to be slow).
+
+``method="invert_cdf"`` (the default) should be sufficient for most cases.
+
Next, we need to specify "fiducial" values for the telescope collecting
area, exposure time, and cosmological redshift. Remember, the initial
photon generation will act as a source for Monte-Carlo sampling for more
@@ -191,12 +214,29 @@
By default, the angular diameter distance to the object is determined
from the ``cosmology`` and the cosmological ``redshift``. If a
``Cosmology`` instance is not provided, one will be made from the
-default cosmological parameters. If your source is local to the galaxy,
-you can set its distance directly, using a tuple, e.g.
-``dist=(30, "kpc")``. In this case, the ``redshift`` and ``cosmology``
-will be ignored. Finally, if the photon generating function accepts any
-parameters, they can be passed to ``from_scratch`` via a ``parameters``
-dictionary.
+default cosmological parameters. The ``center`` keyword argument specifies
+the center of the photon distribution, and the photon positions will be
+rescaled with this value as the origin. This argument accepts the following
+values:
+
+* A NumPy array or list corresponding to the coordinates of the center in
+ units of code length.
+* A ``YTArray`` corresponding to the coordinates of the center in some
+ length units.
+* ``"center"`` or ``"c"`` corresponds to the domain center.
+* ``"max"`` or ``"m"`` corresponds to the location of the maximum gas density.
+* A two-element tuple specifying the max or min of a specific field, e.g.,
+ ``("min","gravitational_potential")``, ``("max","dark_matter_density")``
+
+If ``center`` is not specified, ``from_scratch`` will attempt to use the
+``"center"`` field parameter of the ``data_source``.
+
+``from_scratch`` takes a few other optional keyword arguments. If your
+source is local to the galaxy, you can set its distance directly, using
+a tuple, e.g. ``dist=(30, "kpc")``. In this case, the ``redshift`` and
+``cosmology`` will be ignored. Finally, if the photon generating
+function accepts any parameters, they can be passed to ``from_scratch``
+via a ``parameters`` dictionary.
At this point, the ``photons`` are distributed in the three-dimensional
space of the ``data_source``, with energies in the rest frame of the
@@ -265,7 +305,7 @@
abs_model = TableAbsorbModel("tbabs_table.h5", 0.1)
Now we're ready to project the photons. First, we choose a line-of-sight
-vector ``L``. Second, we'll adjust the exposure time and the redshift.
+vector ``normal``. Second, we'll adjust the exposure time and the redshift.
Third, we'll pass in the absorption ``SpectrumModel``. Fourth, we'll
specify a ``sky_center`` in RA,DEC on the sky in degrees.
@@ -274,26 +314,40 @@
course far short of a full simulation of a telescope ray-trace, but it's
a quick-and-dirty way to get something close to the real thing. We'll
discuss how to get your simulated events into a format suitable for
-reading by telescope simulation codes later.
+reading by telescope simulation codes later. If you just want to convolve
+the photons with an ARF, you may specify that as the only response, but some
+ARFs are unnormalized and still require the RMF for normalization. Check with
+the documentation associated with these files for details. If we are using the
+RMF to convolve energies, we must set ``convolve_energies=True``.
.. code:: python
ARF = "chandra_ACIS-S3_onaxis_arf.fits"
RMF = "chandra_ACIS-S3_onaxis_rmf.fits"
- L = [0.0,0.0,1.0]
- events = photons.project_photons(L, exp_time_new=2.0e5, redshift_new=0.07, absorb_model=abs_model,
- sky_center=(187.5,12.333), responses=[ARF,RMF])
+ normal = [0.0,0.0,1.0]
+ events = photons.project_photons(normal, exp_time_new=2.0e5, redshift_new=0.07, dist_new=None,
+ absorb_model=abs_model, sky_center=(187.5,12.333), responses=[ARF,RMF],
+ convolve_energies=True, no_shifting=False, north_vector=None,
+ psf_sigma=None)
-Also, the optional keyword ``psf_sigma`` specifies a Gaussian standard
-deviation to scatter the photon sky positions around with, providing a
-crude representation of a PSF.
+In this case, we chose a three-vector ``normal`` to specify an arbitrary
+line-of-sight, but ``"x"``, ``"y"``, or ``"z"`` could also be chosen to
+project along one of those axes.
-.. warning::
+``project_photons`` takes several other optional keyword arguments.
- The binned images that result, even if you convolve with responses,
- are still of the same resolution as the finest cell size of the
- simulation dataset. If you want a more accurate simulation of a
- particular X-ray telescope, you should check out `Storing events for future use and for reading-in by telescope simulators`_.
+* ``no_shifting`` (default ``False``) controls whether or not Doppler
+ shifting of photon energies is turned on.
+* ``dist_new`` is a (value, unit) tuple that is used to set a new
+ angular diameter distance by hand instead of having it determined
+ by the cosmology and the value of the redshift. Should only be used
+ for simulations of nearby objects.
+* For off-axis ``normal`` vectors, the ``north_vector`` argument can
+ be used to control what vector corresponds to the "up" direction in
+ the resulting event list.
+* ``psf_sigma`` may be specified to provide a crude representation of
+ a PSF, and corresponds to the standard deviation (in degress) of a
+ Gaussian PSF model.
Let's just take a quick look at the raw events object:
@@ -343,19 +397,27 @@
Which is starting to look like a real observation!
+.. warning::
+
+ The binned images that result, even if you convolve with responses,
+ are still of the same resolution as the finest cell size of the
+ simulation dataset. If you want a more accurate simulation of a
+ particular X-ray telescope, you should check out `Storing events for future use and for reading-in by telescope simulators`_.
+
We can also bin up the spectrum into energy bins, and write it to a FITS
table file. This is an example where we've binned up the spectrum
according to the unconvolved photon energy:
.. code:: python
- events.write_spectrum("virgo_spec.fits", energy_bins=True, emin=0.1, emax=10.0, nchan=2000, clobber=True)
+ events.write_spectrum("virgo_spec.fits", bin_type="energy", emin=0.1, emax=10.0, nchan=2000, clobber=True)
-If we don't set ``energy_bins=True``, and we have convolved our events
+We can also set ``bin_type="channel"``. If we have convolved our events
with response files, then any other keywords will be ignored and it will
try to make a spectrum from the channel information that is contained
-within the RMF, suitable for analyzing in XSPEC. For now, we'll stick
-with the energy spectrum, and plot it up:
+within the RMF. Otherwise, the channels will be determined from the ``emin``,
+``emax``, and ``nchan`` keywords, and will be numbered from 1 to ``nchan``.
+For now, we'll stick with the energy spectrum, and plot it up:
.. code:: python
diff -r c66deb7d69cab3114e9563d6cf9ed57521675ff4 -r 98c8bdfbc1e00d4ee2388663a7f4072bbf98fc57 yt/analysis_modules/photon_simulator/photon_models.py
--- a/yt/analysis_modules/photon_simulator/photon_models.py
+++ b/yt/analysis_modules/photon_simulator/photon_models.py
@@ -25,7 +25,7 @@
from yt.extern.six import string_types
import numpy as np
from yt.funcs import *
-from yt.utilities.physical_constants import mp, kboltz
+from yt.utilities.physical_constants import mp
from yt.utilities.parallel_tools.parallel_analysis_interface import \
parallel_objects
from yt.units.yt_array import uconcatenate
@@ -34,6 +34,12 @@
kT_min = 8.08e-2
kT_max = 50.
+photon_units = {"Energy":"keV",
+ "dx":"kpc"}
+for ax in "xyz":
+ photon_units[ax] = "kpc"
+ photon_units["v"+ax] = "km/s"
+
class PhotonModel(object):
def __init__(self):
@@ -46,7 +52,7 @@
class ThermalPhotonModel(PhotonModel):
r"""
Initialize a ThermalPhotonModel from a thermal spectrum.
-
+
Parameters
----------
@@ -61,30 +67,37 @@
photons_per_chunk : integer
The maximum number of photons that are allocated per chunk. Increase or decrease
as needed.
+ method : string, optional
+ The method used to generate the photon energies from the spectrum:
+ "invert_cdf": Invert the cumulative distribution function of the spectrum.
+ "accept_reject": Acceptance-rejection method using the spectrum.
+ The first method should be sufficient for most cases.
"""
- def __init__(self, spectral_model, X_H=0.75, Zmet=0.3, photons_per_chunk=10000000):
+ def __init__(self, spectral_model, X_H=0.75, Zmet=0.3,
+ photons_per_chunk=10000000, method="invert_cdf"):
self.X_H = X_H
self.Zmet = Zmet
self.spectral_model = spectral_model
self.photons_per_chunk = photons_per_chunk
+ self.method = method
def __call__(self, data_source, parameters):
-
+
ds = data_source.ds
exp_time = parameters["FiducialExposureTime"]
area = parameters["FiducialArea"]
redshift = parameters["FiducialRedshift"]
D_A = parameters["FiducialAngularDiameterDistance"].in_cgs()
- dist_fac = 1.0/(4.*np.pi*D_A.value*D_A.value*(1.+redshift)**3)
+ dist_fac = 1.0/(4.*np.pi*D_A.value*D_A.value*(1.+redshift)**2)
src_ctr = parameters["center"]
- vol_scale = 1.0/np.prod(ds.domain_width.in_cgs().to_ndarray())
-
my_kT_min, my_kT_max = data_source.quantities.extrema("kT")
- self.spectral_model.prepare()
- energy = self.spectral_model.ebins
+ self.spectral_model.prepare_spectrum(redshift)
+ emid = self.spectral_model.emid
+ ebins = self.spectral_model.ebins
+ nchan = len(emid)
citer = data_source.chunks([], "io")
@@ -99,7 +112,13 @@
photons["Energy"] = []
photons["NumberOfPhotons"] = []
- spectral_norm = area.v*exp_time.v*dist_fac/vol_scale
+ spectral_norm = area.v*exp_time.v*dist_fac
+
+ tot_num_cells = data_source.ires.shape[0]
+
+ pbar = get_pbar("Generating photons ", tot_num_cells)
+
+ cell_counter = 0
for chunk in parallel_objects(citer):
@@ -114,7 +133,7 @@
if isinstance(self.Zmet, string_types):
metalZ = chunk[self.Zmet].v
else:
- metalZ = self.Zmet
+ metalZ = self.Zmet*np.ones(num_cells)
idxs = np.argsort(kT)
@@ -133,7 +152,7 @@
n += bcount
kT_idxs = np.unique(kT_idxs)
- cell_em = EM[idxs]*vol_scale
+ cell_em = EM[idxs]*spectral_norm
number_of_photons = np.zeros(num_cells, dtype="uint64")
energies = np.zeros(self.photons_per_chunk)
@@ -141,8 +160,6 @@
start_e = 0
end_e = 0
- pbar = get_pbar("Generating photons for chunk ", num_cells)
-
for ibegin, iend, ikT in zip(bcell, ecell, kT_idxs):
kT = kT_bins[ikT] + 0.5*dkT
@@ -151,58 +168,57 @@
cem = cell_em[ibegin:iend]
- em_sum_c = cem.sum()
- if isinstance(self.Zmet, string_types):
- em_sum_m = (metalZ[ibegin:iend]*cem).sum()
- else:
- em_sum_m = metalZ*em_sum_c
-
cspec, mspec = self.spectral_model.get_spectrum(kT)
- cumspec_c = np.cumsum(cspec.d)
- tot_ph_c = cumspec_c[-1]*spectral_norm*em_sum_c
- cumspec_c /= cumspec_c[-1]
- cumspec_c = np.insert(cumspec_c, 0, 0.0)
-
- cumspec_m = np.cumsum(mspec.d)
- tot_ph_m = cumspec_m[-1]*spectral_norm*em_sum_m
- cumspec_m /= cumspec_m[-1]
- cumspec_m = np.insert(cumspec_m, 0, 0.0)
+ tot_ph_c = cspec.d.sum()
+ tot_ph_m = mspec.d.sum()
u = np.random.random(size=n_current)
- cell_norm_c = tot_ph_c*cem/em_sum_c
- cell_n_c = np.uint64(cell_norm_c) + np.uint64(np.modf(cell_norm_c)[0] >= u)
-
- if isinstance(self.Zmet, string_types):
- cell_norm_m = tot_ph_m*metalZ[ibegin:iend]*cem/em_sum_m
- else:
- cell_norm_m = tot_ph_m*metalZ*cem/em_sum_m
- cell_n_m = np.uint64(cell_norm_m) + np.uint64(np.modf(cell_norm_m)[0] >= u)
-
- number_of_photons[ibegin:iend] = cell_n_c + cell_n_m
+ cell_norm_c = tot_ph_c*cem
+ cell_norm_m = tot_ph_m*metalZ[ibegin:iend]*cem
+ cell_norm = np.modf(cell_norm_c + cell_norm_m)
+ cell_n = np.uint64(cell_norm[1]) + np.uint64(cell_norm[0] >= u)
- end_e += int((cell_n_c+cell_n_m).sum())
+ number_of_photons[ibegin:iend] = cell_n
+
+ end_e += int(cell_n.sum())
if end_e > self.photons_per_chunk:
raise RuntimeError("Number of photons generated for this chunk "+
"exceeds photons_per_chunk (%d)! " % self.photons_per_chunk +
"Increase photons_per_chunk!")
- energies[start_e:end_e] = _generate_energies(cell_n_c, cell_n_m, cumspec_c, cumspec_m, energy)
-
+ if self.method == "invert_cdf":
+ cumspec_c = np.cumsum(cspec.d)
+ cumspec_m = np.cumsum(mspec.d)
+ cumspec_c = np.insert(cumspec_c, 0, 0.0)
+ cumspec_m = np.insert(cumspec_m, 0, 0.0)
+
+ ei = start_e
+ for cn, Z in zip(number_of_photons[ibegin:iend], metalZ[ibegin:iend]):
+ if cn == 0: continue
+ if self.method == "invert_cdf":
+ cumspec = cumspec_c + Z*cumspec_m
+ cumspec /= cumspec[-1]
+ randvec = np.random.uniform(size=cn)
+ randvec.sort()
+ cell_e = np.interp(randvec, cumspec, ebins)
+ elif self.method == "accept_reject":
+ tot_spec = cspec.d+Z*mspec.d
+ tot_spec /= tot_spec.sum()
+ eidxs = np.random.choice(nchan, size=cn, p=tot_spec)
+ cell_e = emid[eidxs]
+ energies[ei:ei+cn] = cell_e
+ cell_counter += 1
+ pbar.update(cell_counter)
+ ei += cn
+
start_e = end_e
- pbar.update(iend)
-
- pbar.finish()
-
active_cells = number_of_photons > 0
idxs = idxs[active_cells]
- mylog.info("Number of photons generated for this chunk: %d" % int(number_of_photons.sum()))
- mylog.info("Number of cells with photons: %d" % int(active_cells.sum()))
-
photons["NumberOfPhotons"].append(number_of_photons[active_cells])
photons["Energy"].append(ds.arr(energies[:end_e].copy(), "keV"))
photons["x"].append((chunk["x"][idxs]-src_ctr[0]).in_units("kpc"))
@@ -213,23 +229,17 @@
photons["vz"].append(chunk["velocity_z"][idxs].in_units("km/s"))
photons["dx"].append(chunk["dx"][idxs].in_units("kpc"))
+ pbar.finish()
+
for key in photons:
if len(photons[key]) > 0:
photons[key] = uconcatenate(photons[key])
+ elif key == "NumberOfPhotons":
+ photons[key] = np.array([])
+ else:
+ photons[key] = YTArray([], photon_units[key])
+
+ mylog.info("Number of photons generated: %d" % int(np.sum(photons["NumberOfPhotons"])))
+ mylog.info("Number of cells with photons: %d" % len(photons["x"]))
return photons
-
-def _generate_energies(cell_n_c, cell_n_m, counts_c, counts_m, energy):
- energies = np.array([])
- for cn_c, cn_m in zip(cell_n_c, cell_n_m):
- if cn_c > 0:
- randvec_c = np.random.uniform(size=cn_c)
- randvec_c.sort()
- cell_e_c = np.interp(randvec_c, counts_c, energy)
- energies = np.append(energies, cell_e_c)
- if cn_m > 0:
- randvec_m = np.random.uniform(size=cn_m)
- randvec_m.sort()
- cell_e_m = np.interp(randvec_m, counts_m, energy)
- energies = np.append(energies, cell_e_m)
- return energies
diff -r c66deb7d69cab3114e9563d6cf9ed57521675ff4 -r 98c8bdfbc1e00d4ee2388663a7f4072bbf98fc57 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
@@ -1,5 +1,10 @@
"""
Classes for generating lists of photons and detected events
+
+The SciPy Proceeding that describes this module in detail may be found at:
+
+http://conference.scipy.org/proceedings/scipy2014/zuhone.html
+
The algorithms used here are based off of the method used by the
PHOX code (http://www.mpa-garching.mpg.de/~kdolag/Phox/),
developed by Veronica Biffi and Klaus Dolag. References for
@@ -25,15 +30,21 @@
from yt.utilities.physical_constants import clight
from yt.utilities.cosmology import Cosmology
from yt.utilities.orientation import Orientation
+from yt.utilities.fits_image import assert_same_wcs
from yt.utilities.parallel_tools.parallel_analysis_interface import \
communication_system, parallel_root_only, get_mpi_type, \
parallel_capable
from yt.units.yt_array import YTQuantity, YTArray, uconcatenate
import h5py
from yt.utilities.on_demand_imports import _astropy
+import warnings
comm = communication_system.communicators[-1]
+axes_lookup = {"x":("y","z"),
+ "y":("z","x"),
+ "z":("x","y")}
+
def parse_value(value, default_units):
if isinstance(value, YTQuantity):
return value.in_units(default_units)
@@ -50,10 +61,10 @@
self.cosmo = cosmo
self.p_bins = p_bins
self.num_cells = len(photons["x"])
-
+
def keys(self):
return self.photons.keys()
-
+
def items(self):
ret = []
for k, v in self.photons.items():
@@ -62,7 +73,7 @@
else:
ret.append((k,v))
return ret
-
+
def values(self):
ret = []
for k, v in self.photons.items():
@@ -71,7 +82,7 @@
else:
ret.append(v)
return ret
-
+
def __getitem__(self, key):
if key == "Energy":
return [self.photons["Energy"][self.p_bins[i]:self.p_bins[i+1]]
@@ -127,9 +138,9 @@
p_bins = np.cumsum(photons["NumberOfPhotons"])
p_bins = np.insert(p_bins, 0, [np.uint64(0)])
-
+
photons["Energy"] = YTArray(f["/energy"][start_e:end_e], "keV")
-
+
f.close()
cosmo = Cosmology(hubble_constant=parameters["HubbleConstant"],
@@ -137,7 +148,7 @@
omega_lambda=parameters["OmegaLambda"])
return cls(photons, parameters, cosmo, p_bins)
-
+
@classmethod
def from_scratch(cls, data_source, redshift, area,
exp_time, photon_model, parameters=None,
@@ -277,18 +288,25 @@
D_A = parse_value(dist, "Mpc")
redshift = 0.0
- if center == "c":
+ if center in ("center", "c"):
parameters["center"] = ds.domain_center
- elif center == "max":
+ elif center in ("max", "m"):
parameters["center"] = ds.find_max("density")[-1]
elif iterable(center):
if isinstance(center, YTArray):
parameters["center"] = center.in_units("code_length")
+ elif isinstance(center, tuple):
+ if center[0] == "min":
+ parameters["center"] = ds.find_min(center[1])[-1]
+ elif center[0] == "max":
+ parameters["center"] = ds.find_max(center[1])[-1]
+ else:
+ raise RuntimeError
else:
parameters["center"] = ds.arr(center, "code_length")
elif center is None:
parameters["center"] = data_source.get_field_parameter("center")
-
+
parameters["FiducialExposureTime"] = parse_value(exp_time, "s")
parameters["FiducialArea"] = parse_value(area, "cm**2")
parameters["FiducialRedshift"] = redshift
@@ -308,32 +326,32 @@
dimension = max(dimension, int(width/delta_min))
parameters["Dimension"] = 2*dimension
parameters["Width"] = 2.*width.in_units("kpc")
-
+
photons = photon_model(data_source, parameters)
mylog.info("Finished generating photons.")
p_bins = np.cumsum(photons["NumberOfPhotons"])
p_bins = np.insert(p_bins, 0, [np.uint64(0)])
-
+
return cls(photons, parameters, cosmo, p_bins)
-
+
def write_h5_file(self, photonfile):
"""
Write the photons to the HDF5 file *photonfile*.
"""
if parallel_capable:
-
+
mpi_long = get_mpi_type("int64")
mpi_double = get_mpi_type("float64")
-
+
local_num_cells = len(self.photons["x"])
sizes_c = comm.comm.gather(local_num_cells, root=0)
-
- local_num_photons = self.photons["NumberOfPhotons"].sum()
+
+ local_num_photons = np.sum(self.photons["NumberOfPhotons"])
sizes_p = comm.comm.gather(local_num_photons, root=0)
-
+
if comm.rank == 0:
num_cells = sum(sizes_c)
num_photons = sum(sizes_p)
@@ -362,24 +380,24 @@
dx = np.empty([])
n_ph = np.empty([])
e = np.empty([])
-
- comm.comm.Gatherv([self.photons["x"].ndarray_view(), local_num_cells, mpi_double],
+
+ comm.comm.Gatherv([self.photons["x"].d, local_num_cells, mpi_double],
[x, (sizes_c, disps_c), mpi_double], root=0)
- comm.comm.Gatherv([self.photons["y"].ndarray_view(), local_num_cells, mpi_double],
+ comm.comm.Gatherv([self.photons["y"].d, local_num_cells, mpi_double],
[y, (sizes_c, disps_c), mpi_double], root=0)
- comm.comm.Gatherv([self.photons["z"].ndarray_view(), local_num_cells, mpi_double],
+ comm.comm.Gatherv([self.photons["z"].d, local_num_cells, mpi_double],
[z, (sizes_c, disps_c), mpi_double], root=0)
- comm.comm.Gatherv([self.photons["vx"].ndarray_view(), local_num_cells, mpi_double],
+ comm.comm.Gatherv([self.photons["vx"].d, local_num_cells, mpi_double],
[vx, (sizes_c, disps_c), mpi_double], root=0)
- comm.comm.Gatherv([self.photons["vy"].ndarray_view(), local_num_cells, mpi_double],
+ comm.comm.Gatherv([self.photons["vy"].d, local_num_cells, mpi_double],
[vy, (sizes_c, disps_c), mpi_double], root=0)
- comm.comm.Gatherv([self.photons["vz"].ndarray_view(), local_num_cells, mpi_double],
+ comm.comm.Gatherv([self.photons["vz"].d, local_num_cells, mpi_double],
[vz, (sizes_c, disps_c), mpi_double], root=0)
- comm.comm.Gatherv([self.photons["dx"].ndarray_view(), local_num_cells, mpi_double],
+ comm.comm.Gatherv([self.photons["dx"].d, local_num_cells, mpi_double],
[dx, (sizes_c, disps_c), mpi_double], root=0)
comm.comm.Gatherv([self.photons["NumberOfPhotons"], local_num_cells, mpi_long],
[n_ph, (sizes_c, disps_c), mpi_long], root=0)
- comm.comm.Gatherv([self.photons["Energy"].ndarray_view(), local_num_photons, mpi_double],
+ comm.comm.Gatherv([self.photons["Energy"].d, local_num_photons, mpi_double],
[e, (sizes_p, disps_p), mpi_double], root=0)
else:
@@ -393,13 +411,13 @@
dx = self.photons["dx"].d
n_ph = self.photons["NumberOfPhotons"]
e = self.photons["Energy"].d
-
+
if comm.rank == 0:
-
+
f = h5py.File(photonfile, "w")
# Scalars
-
+
f.create_dataset("fid_area", data=float(self.parameters["FiducialArea"]))
f.create_dataset("fid_exp_time", data=float(self.parameters["FiducialExposureTime"]))
f.create_dataset("fid_redshift", data=self.parameters["FiducialRedshift"])
@@ -409,7 +427,7 @@
f.create_dataset("fid_d_a", data=float(self.parameters["FiducialAngularDiameterDistance"]))
f.create_dataset("dimension", data=self.parameters["Dimension"])
f.create_dataset("width", data=float(self.parameters["Width"]))
-
+
# Arrays
f.create_dataset("x", data=x)
@@ -426,19 +444,22 @@
comm.barrier()
- def project_photons(self, L, area_new=None, exp_time_new=None,
+ def project_photons(self, normal, area_new=None, exp_time_new=None,
redshift_new=None, dist_new=None,
absorb_model=None, psf_sigma=None,
sky_center=None, responses=None,
- convolve_energies=False, no_shifting=False):
+ convolve_energies=False, no_shifting=False,
+ north_vector=None):
r"""
Projects photons onto an image plane given a line of sight.
Parameters
----------
- L : array_like
- Normal vector to the plane of projection.
+ normal : character or array_like
+ Normal vector to the plane of projection. If "x", "y", or "z", will
+ assume to be along that axis (and will probably be faster). Otherwise,
+ should be an off-axis normal vector, e.g [1.0,2.0,-3.0]
area_new : float, optional
New value for the effective area of the detector. If *responses*
are specified the value of this keyword is ignored.
@@ -463,7 +484,11 @@
If this is set, the photon energies will be convolved with the RMF.
no_shifting : boolean, optional
If set, the photon energies will not be Doppler shifted.
-
+ north_vector : a sequence of floats
+ A vector defining the 'up' direction. This option sets the orientation of
+ the plane of projection. If not set, an arbitrary grid-aligned north_vector
+ is chosen. Ignored in the case where a particular axis (e.g., "x", "y", or
+ "z") is explicitly specified.
Examples
--------
>>> L = np.array([0.1,-0.2,0.3])
@@ -475,7 +500,7 @@
if redshift_new is not None and dist_new is not None:
mylog.error("You may specify a new redshift or distance, "+
"but not both!")
-
+
if sky_center is None:
sky_center = YTArray([30.,45.], "degree")
else:
@@ -486,37 +511,36 @@
if psf_sigma is not None:
psf_sigma = parse_value(psf_sigma, "degree")
- L = np.array(L)
- L /= np.sqrt(np.dot(L, L))
- vecs = np.identity(3)
- t = np.cross(L, vecs).sum(axis=1)
- ax = t.argmax()
- north = np.cross(L, vecs[ax,:]).ravel()
- orient = Orientation(L, north_vector=north)
-
- x_hat = orient.unit_vectors[0]
- y_hat = orient.unit_vectors[1]
- z_hat = orient.unit_vectors[2]
+ if not isinstance(normal, string_types):
+ L = np.array(normal)
+ orient = Orientation(L, north_vector=north_vector)
+ x_hat = orient.unit_vectors[0]
+ y_hat = orient.unit_vectors[1]
+ z_hat = orient.unit_vectors[2]
n_ph = self.photons["NumberOfPhotons"]
n_ph_tot = n_ph.sum()
-
+
eff_area = None
parameters = {}
-
+
if responses is not None:
responses = ensure_list(responses)
parameters["ARF"] = responses[0]
if len(responses) == 2:
parameters["RMF"] = responses[1]
area_new = parameters["ARF"]
-
+
+ zobs0 = self.parameters["FiducialRedshift"]
+ D_A0 = self.parameters["FiducialAngularDiameterDistance"]
+ scale_factor = 1.0
+
if (exp_time_new is None and area_new is None and
redshift_new is None and dist_new is None):
my_n_obs = n_ph_tot
- zobs = self.parameters["FiducialRedshift"]
- D_A = self.parameters["FiducialAngularDiameterDistance"]
+ zobs = zobs0
+ D_A = D_A0
else:
if exp_time_new is None:
Tratio = 1.
@@ -545,8 +569,8 @@
Aratio = parse_value(area_new, "cm**2")/self.parameters["FiducialArea"]
if redshift_new is None and dist_new is None:
Dratio = 1.
- zobs = self.parameters["FiducialRedshift"]
- D_A = self.parameters["FiducialAngularDiameterDistance"]
+ zobs = zobs0
+ D_A = D_A0
else:
if redshift_new is None:
zobs = 0.0
@@ -554,28 +578,17 @@
else:
zobs = redshift_new
D_A = self.cosmo.angular_diameter_distance(0.0,zobs).in_units("Mpc")
- fid_D_A = self.parameters["FiducialAngularDiameterDistance"]
- Dratio = fid_D_A*fid_D_A*(1.+self.parameters["FiducialRedshift"]**3) / \
+ scale_factor = (1.+zobs0)/(1.+zobs)
+ Dratio = D_A0*D_A0*(1.+zobs0)**3 / \
(D_A*D_A*(1.+zobs)**3)
fak = Aratio*Tratio*Dratio
if fak > 1:
- raise ValueError("Spectrum scaling factor = %g, cannot be greater than unity." % (fak))
+ raise ValueError("Spectrum scaling factor = %g, cannot be greater than unity." % fak)
my_n_obs = np.uint64(n_ph_tot*fak)
n_obs_all = comm.mpi_allreduce(my_n_obs)
if comm.rank == 0:
mylog.info("Total number of photons to use: %d" % (n_obs_all))
-
- x = np.random.uniform(low=-0.5,high=0.5,size=my_n_obs)
- y = np.random.uniform(low=-0.5,high=0.5,size=my_n_obs)
- z = np.random.uniform(low=-0.5,high=0.5,size=my_n_obs)
-
- if not no_shifting:
- vz = self.photons["vx"]*z_hat[0] + \
- self.photons["vy"]*z_hat[1] + \
- self.photons["vz"]*z_hat[2]
- shift = -vz.in_cgs()/clight
- shift = np.sqrt((1.-shift)/(1.+shift))
if my_n_obs == n_ph_tot:
idxs = np.arange(my_n_obs,dtype='uint64')
@@ -584,32 +597,57 @@
obs_cells = np.searchsorted(self.p_bins, idxs, side='right')-1
delta = dx[obs_cells]
- x *= delta
- y *= delta
- z *= delta
- x += self.photons["x"][obs_cells]
- y += self.photons["y"][obs_cells]
- z += self.photons["z"][obs_cells]
+ if isinstance(normal, string_types):
+
+ xsky = np.random.uniform(low=-0.5,high=0.5,size=my_n_obs)
+ ysky = np.random.uniform(low=-0.5,high=0.5,size=my_n_obs)
+ xsky *= delta
+ ysky *= delta
+ xsky += self.photons[axes_lookup[normal][0]][obs_cells]
+ ysky += self.photons[axes_lookup[normal][1]][obs_cells]
+
+ if not no_shifting:
+ vz = self.photons["v%s" % normal]
+
+ else:
+ x = np.random.uniform(low=-0.5,high=0.5,size=my_n_obs)
+ y = np.random.uniform(low=-0.5,high=0.5,size=my_n_obs)
+ z = np.random.uniform(low=-0.5,high=0.5,size=my_n_obs)
+
+ if not no_shifting:
+ vz = self.photons["vx"]*z_hat[0] + \
+ self.photons["vy"]*z_hat[1] + \
+ self.photons["vz"]*z_hat[2]
+
+ x *= delta
+ y *= delta
+ z *= delta
+ x += self.photons["x"][obs_cells]
+ y += self.photons["y"][obs_cells]
+ z += self.photons["z"][obs_cells]
+
+ xsky = x*x_hat[0] + y*x_hat[1] + z*x_hat[2]
+ ysky = x*y_hat[0] + y*y_hat[1] + z*y_hat[2]
+
if no_shifting:
eobs = self.photons["Energy"][idxs]
else:
+ shift = -vz.in_cgs()/clight
+ shift = np.sqrt((1.-shift)/(1.+shift))
eobs = self.photons["Energy"][idxs]*shift[obs_cells]
+ eobs *= scale_factor
- xsky = x*x_hat[0] + y*x_hat[1] + z*x_hat[2]
- ysky = x*y_hat[0] + y*y_hat[1] + z*y_hat[2]
- eobs /= (1.+zobs)
-
if absorb_model is None:
not_abs = np.ones(eobs.shape, dtype='bool')
else:
mylog.info("Absorbing.")
- absorb_model.prepare()
+ absorb_model.prepare_spectrum()
emid = absorb_model.emid
aspec = absorb_model.get_spectrum()
absorb = np.interp(eobs, emid, aspec, left=0.0, right=0.0)
randvec = aspec.max()*np.random.random(eobs.shape)
not_abs = randvec < absorb
-
+
if eff_area is None:
detected = np.ones(eobs.shape, dtype='bool')
else:
@@ -618,14 +656,14 @@
earea = np.interp(eobs, earf, eff_area, left=0.0, right=0.0)
randvec = eff_area.max()*np.random.random(eobs.shape)
detected = randvec < earea
-
+
detected = np.logical_and(not_abs, detected)
-
+
events = {}
dx_min = self.parameters["Width"]/self.parameters["Dimension"]
dtheta = YTQuantity(np.rad2deg(dx_min/D_A), "degree")
-
+
events["xpix"] = xsky[detected]/dx_min.v + 0.5*(nx+1)
events["ypix"] = ysky[detected]/dx_min.v + 0.5*(nx+1)
events["eobs"] = eobs[detected]
@@ -635,15 +673,15 @@
if psf_sigma is not None:
events["xpix"] += np.random.normal(sigma=psf_sigma/dtheta)
events["ypix"] += np.random.normal(sigma=psf_sigma/dtheta)
-
+
num_events = len(events["xpix"])
-
+
if comm.rank == 0: mylog.info("Total number of observed photons: %d" % num_events)
if "RMF" in parameters and convolve_energies:
events, info = self._convolve_with_rmf(parameters["RMF"], events)
for k, v in info.items(): parameters[k] = v
-
+
if exp_time_new is None:
parameters["ExposureTime"] = self.parameters["FiducialExposureTime"]
else:
@@ -657,7 +695,7 @@
parameters["sky_center"] = sky_center
parameters["pix_center"] = np.array([0.5*(nx+1)]*2)
parameters["dtheta"] = dtheta
-
+
return EventList(events, parameters)
def _normalize_arf(self, respfile):
@@ -688,7 +726,7 @@
eidxs = np.argsort(events["eobs"])
- phEE = events["eobs"][eidxs].ndarray_view()
+ phEE = events["eobs"][eidxs].d
phXX = events["xpix"][eidxs]
phYY = events["ypix"][eidxs]
@@ -758,7 +796,7 @@
self.num_events = events["xpix"].shape[0]
self.wcs = _astropy.pywcs.WCS(naxis=2)
self.wcs.wcs.crpix = parameters["pix_center"]
- self.wcs.wcs.crval = parameters["sky_center"].ndarray_view()
+ self.wcs.wcs.crval = parameters["sky_center"].d
self.wcs.wcs.cdelt = [-parameters["dtheta"].value, parameters["dtheta"].value]
self.wcs.wcs.ctype = ["RA---TAN","DEC--TAN"]
self.wcs.wcs.cunit = ["deg"]*2
@@ -785,8 +823,9 @@
return self.events.__repr__()
def __add__(self, other):
- keys1 = self.parameters.keys()
- keys2 = other.parameters.keys()
+ assert_same_wcs(self.wcs, other.wcs)
+ keys1 = list(self.parameters.keys())
+ keys2 = list(other.parameters.keys())
keys1.sort()
keys2.sort()
if keys1 != keys2:
@@ -809,9 +848,9 @@
return EventList(events, self.parameters)
def filter_events(self, region):
- """
- Filter events using a ds9 region. Requires the pyregion package.
- Returns a new EventList.
+ """
+ Filter events using a ds9 region. Requires the pyregion package.
+ Returns a new EventList.
"""
import pyregion
import os
@@ -836,7 +875,7 @@
"""
events = {}
parameters = {}
-
+
f = h5py.File(h5file, "r")
parameters["ExposureTime"] = YTQuantity(f["/exp_time"].value, "s")
@@ -923,7 +962,7 @@
cols = []
col1 = pyfits.Column(name='ENERGY', format='E', unit='eV',
- array=self.events["eobs"].in_units("eV").ndarray_view())
+ array=self.events["eobs"].in_units("eV").d)
col2 = pyfits.Column(name='X', format='D', unit='pixel',
array=self.events["xpix"])
col3 = pyfits.Column(name='Y', format='D', unit='pixel',
@@ -942,7 +981,7 @@
cols.append(col4)
coldefs = pyfits.ColDefs(cols)
- tbhdu = pyfits.new_table(coldefs)
+ tbhdu = pyfits.BinTableHDU.from_columns(coldefs)
tbhdu.update_ext_name("EVENTS")
tbhdu.header["MTYPE1"] = "sky"
@@ -1006,29 +1045,29 @@
"""
pyfits = _astropy.pyfits
if isinstance(self.parameters["Area"], string_types):
- mylog.error("Writing SIMPUT files is only supported if you didn't convolve with an ARF.")
- raise TypeError("Writing SIMPUT files is only supported if you didn't convolve with an ARF.")
+ mylog.error("Writing SIMPUT files is only supported if you didn't convolve with responses.")
+ raise TypeError("Writing SIMPUT files is only supported if you didn't convolve with responses.")
if emin is None:
- emin = self.events["eobs"].min().value*0.95
+ emin = self.events["eobs"].min().value
if emax is None:
- emax = self.events["eobs"].max().value*1.05
+ emax = self.events["eobs"].max().value
- idxs = np.logical_and(self.events["eobs"].ndarray_view() >= emin,
- self.events["eobs"].ndarray_view() <= emax)
+ idxs = np.logical_and(self.events["eobs"].d >= emin,
+ self.events["eobs"].d <= emax)
flux = np.sum(self.events["eobs"][idxs].in_units("erg")) / \
self.parameters["ExposureTime"]/self.parameters["Area"]
col1 = pyfits.Column(name='ENERGY', format='E',
- array=self["eobs"].ndarray_view())
+ array=self["eobs"].d)
col2 = pyfits.Column(name='DEC', format='D',
- array=self["ysky"].ndarray_view())
+ array=self["ysky"].d)
col3 = pyfits.Column(name='RA', format='D',
- array=self["xsky"].ndarray_view())
+ array=self["xsky"].d)
coldefs = pyfits.ColDefs([col1, col2, col3])
- tbhdu = pyfits.new_table(coldefs)
+ tbhdu = pyfits.BinTableHDU.from_columns(coldefs)
tbhdu.update_ext_name("PHLIST")
tbhdu.header["HDUCLASS"] = "HEASARC/SIMPUT"
@@ -1056,7 +1095,7 @@
coldefs = pyfits.ColDefs([col1, col2, col3, col4, col5, col6, col7, col8])
- wrhdu = pyfits.new_table(coldefs)
+ wrhdu = pyfits.BinTableHDU.from_columns(coldefs)
wrhdu.update_ext_name("SRC_CAT")
wrhdu.header["HDUCLASS"] = "HEASARC"
@@ -1104,14 +1143,14 @@
f.create_dataset("/xpix", data=self.events["xpix"])
f.create_dataset("/ypix", data=self.events["ypix"])
- f.create_dataset("/xsky", data=self.events["xsky"].ndarray_view())
- f.create_dataset("/ysky", data=self.events["ysky"].ndarray_view())
- f.create_dataset("/eobs", data=self.events["eobs"].ndarray_view())
+ f.create_dataset("/xsky", data=self.events["xsky"].d)
+ f.create_dataset("/ysky", data=self.events["ysky"].d)
+ f.create_dataset("/eobs", data=self.events["eobs"].d)
if "PI" in self.events:
f.create_dataset("/pi", data=self.events["PI"])
if "PHA" in self.events:
- f.create_dataset("/pha", data=self.events["PHA"])
- f.create_dataset("/sky_center", data=self.parameters["sky_center"].ndarray_view())
+ f.create_dataset("/pha", data=self.events["PHA"])
+ f.create_dataset("/sky_center", data=self.parameters["sky_center"].d)
f.create_dataset("/pix_center", data=self.parameters["pix_center"])
f.create_dataset("/dtheta", data=float(self.parameters["dtheta"]))
@@ -1138,11 +1177,11 @@
if emin is None:
mask_emin = np.ones((self.num_events), dtype='bool')
else:
- mask_emin = self.events["eobs"].ndarray_view() > emin
+ mask_emin = self.events["eobs"].d > emin
if emax is None:
mask_emax = np.ones((self.num_events), dtype='bool')
else:
- mask_emax = self.events["eobs"].ndarray_view() < emax
+ mask_emax = self.events["eobs"].d < emax
mask = np.logical_and(mask_emin, mask_emax)
@@ -1175,40 +1214,37 @@
hdu.writeto(imagefile, clobber=clobber)
@parallel_root_only
- def write_spectrum(self, specfile, energy_bins=False, emin=0.1,
- emax=10.0, nchan=2000, clobber=False):
+ def write_spectrum(self, specfile, bin_type="channel", emin=0.1,
+ emax=10.0, nchan=2000, clobber=False, energy_bins=False):
r"""
Bin event energies into a spectrum and write it to a FITS binary table. Can bin
- on energy or channel, the latter only if the photons were convolved with an
- RMF. In that case, the spectral binning will be determined by the RMF binning.
+ on energy or channel. In that case, the spectral binning will be determined by
+ the RMF binning.
Parameters
----------
specfile : string
The name of the FITS file to be written.
+ bin_type : string, optional
+ Bin on "energy" or "channel". If an RMF is detected, channel information will be
+ imported from it.
+ emin : float, optional
+ The minimum energy of the spectral bins in keV. Only used if binning without an RMF.
+ emax : float, optional
+ The maximum energy of the spectral bins in keV. Only used if binning without an RMF.
+ nchan : integer, optional
+ The number of channels. Only used if binning without an RMF.
energy_bins : boolean, optional
- Bin on energy or channel.
- emin : float, optional
- The minimum energy of the spectral bins in keV. Only used if binning on energy.
- emax : float, optional
- The maximum energy of the spectral bins in keV. Only used if binning on energy.
- nchan : integer, optional
- The number of channels. Only used if binning on energy.
+ Bin on energy or channel. Deprecated in favor of *bin_type*.
"""
+ if energy_bins:
+ bin_type = "energy"
+ warnings.warn("The energy_bins keyword is deprecated. Please use "
+ "the bin_type keyword instead. Setting bin_type == 'energy'.")
pyfits = _astropy.pyfits
- if energy_bins:
- spectype = "energy"
- espec = self.events["eobs"].d
- range = (emin, emax)
- spec, ee = np.histogram(espec, bins=nchan, range=range)
- bins = 0.5*(ee[1:]+ee[:-1])
- else:
- if "ChannelType" not in self.parameters:
- mylog.error("These events were not convolved with an RMF. Set energy_bins=True.")
- raise KeyError
+ if bin_type == "channel" and "ChannelType" in self.parameters:
spectype = self.parameters["ChannelType"]
- espec = self.events[spectype]
f = pyfits.open(self.parameters["RMF"])
nchan = int(f[1].header["DETCHANS"])
try:
@@ -1217,18 +1253,22 @@
mylog.warning("Cannot determine minimum allowed value for channel. " +
"Setting to 0, which may be wrong.")
cmin = int(0)
- try:
- cmax = int(f[1].header["TLMAX4"])
- except KeyError:
- mylog.warning("Cannot determine maximum allowed value for channel. " +
- "Setting to DETCHANS, which may be wrong.")
- cmax = int(nchan)
f.close()
minlength = nchan
if cmin == 1: minlength += 1
spec = np.bincount(self.events[spectype],minlength=minlength)
if cmin == 1: spec = spec[1:]
- bins = np.arange(nchan).astype("int32")+cmin
+ bins = (np.arange(nchan)+cmin).astype("int32")
+ else:
+ espec = self.events["eobs"].d
+ erange = (emin, emax)
+ spec, ee = np.histogram(espec, bins=nchan, range=erange)
+ if bin_type == "energy":
+ bins = 0.5*(ee[1:]+ee[:-1])
+ spectype = "energy"
+ else:
+ bins = (np.arange(nchan)+1).astype("int32")
+ spectype = "pi"
col1 = pyfits.Column(name='CHANNEL', format='1J', array=bins)
col2 = pyfits.Column(name=spectype.upper(), format='1D', array=bins.astype("float64"))
@@ -1240,46 +1280,45 @@
tbhdu = pyfits.new_table(coldefs)
tbhdu.update_ext_name("SPECTRUM")
- if not energy_bins:
- tbhdu.header["DETCHANS"] = spec.shape[0]
- tbhdu.header["TOTCTS"] = spec.sum()
- tbhdu.header["EXPOSURE"] = float(self.parameters["ExposureTime"])
- tbhdu.header["LIVETIME"] = float(self.parameters["ExposureTime"])
- tbhdu.header["CONTENT"] = spectype
- tbhdu.header["HDUCLASS"] = "OGIP"
- tbhdu.header["HDUCLAS1"] = "SPECTRUM"
- tbhdu.header["HDUCLAS2"] = "TOTAL"
- tbhdu.header["HDUCLAS3"] = "TYPE:I"
- tbhdu.header["HDUCLAS4"] = "COUNT"
- tbhdu.header["HDUVERS"] = "1.1.0"
- tbhdu.header["HDUVERS1"] = "1.1.0"
- tbhdu.header["CHANTYPE"] = spectype
- tbhdu.header["BACKFILE"] = "none"
- tbhdu.header["CORRFILE"] = "none"
- tbhdu.header["POISSERR"] = True
- if "RMF" in self.parameters:
- tbhdu.header["RESPFILE"] = self.parameters["RMF"]
- else:
- tbhdu.header["RESPFILE"] = "none"
- if "ARF" in self.parameters:
- tbhdu.header["ANCRFILE"] = self.parameters["ARF"]
- else:
- tbhdu.header["ANCRFILE"] = "none"
- if "Mission" in self.parameters:
- tbhdu.header["MISSION"] = self.parameters["Mission"]
- else:
- tbhdu.header["MISSION"] = "none"
- if "Telescope" in self.parameters:
- tbhdu.header["TELESCOP"] = self.parameters["Telescope"]
- else:
- tbhdu.header["TELESCOP"] = "none"
- if "Instrument" in self.parameters:
- tbhdu.header["INSTRUME"] = self.parameters["Instrument"]
- else:
- tbhdu.header["INSTRUME"] = "none"
- tbhdu.header["AREASCAL"] = 1.0
- tbhdu.header["CORRSCAL"] = 0.0
- tbhdu.header["BACKSCAL"] = 1.0
+ tbhdu.header["DETCHANS"] = spec.shape[0]
+ tbhdu.header["TOTCTS"] = spec.sum()
+ tbhdu.header["EXPOSURE"] = float(self.parameters["ExposureTime"])
+ tbhdu.header["LIVETIME"] = float(self.parameters["ExposureTime"])
+ tbhdu.header["CONTENT"] = spectype
+ tbhdu.header["HDUCLASS"] = "OGIP"
+ tbhdu.header["HDUCLAS1"] = "SPECTRUM"
+ tbhdu.header["HDUCLAS2"] = "TOTAL"
+ tbhdu.header["HDUCLAS3"] = "TYPE:I"
+ tbhdu.header["HDUCLAS4"] = "COUNT"
+ tbhdu.header["HDUVERS"] = "1.1.0"
+ tbhdu.header["HDUVERS1"] = "1.1.0"
+ tbhdu.header["CHANTYPE"] = spectype
+ tbhdu.header["BACKFILE"] = "none"
+ tbhdu.header["CORRFILE"] = "none"
+ tbhdu.header["POISSERR"] = True
+ if "RMF" in self.parameters:
+ tbhdu.header["RESPFILE"] = self.parameters["RMF"]
+ else:
+ tbhdu.header["RESPFILE"] = "none"
+ if "ARF" in self.parameters:
+ tbhdu.header["ANCRFILE"] = self.parameters["ARF"]
+ else:
+ tbhdu.header["ANCRFILE"] = "none"
+ if "Mission" in self.parameters:
+ tbhdu.header["MISSION"] = self.parameters["Mission"]
+ else:
+ tbhdu.header["MISSION"] = "none"
+ if "Telescope" in self.parameters:
+ tbhdu.header["TELESCOP"] = self.parameters["Telescope"]
+ else:
+ tbhdu.header["TELESCOP"] = "none"
+ if "Instrument" in self.parameters:
+ tbhdu.header["INSTRUME"] = self.parameters["Instrument"]
+ else:
+ tbhdu.header["INSTRUME"] = "none"
+ tbhdu.header["AREASCAL"] = 1.0
+ tbhdu.header["CORRSCAL"] = 0.0
+ tbhdu.header["BACKSCAL"] = 1.0
hdulist = pyfits.HDUList([pyfits.PrimaryHDU(), tbhdu])
diff -r c66deb7d69cab3114e9563d6cf9ed57521675ff4 -r 98c8bdfbc1e00d4ee2388663a7f4072bbf98fc57 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
@@ -32,17 +32,17 @@
self.ebins = np.linspace(self.emin, self.emax, nchan+1)
self.de = np.diff(self.ebins)
self.emid = 0.5*(self.ebins[1:]+self.ebins[:-1])
-
- def prepare(self):
+
+ def prepare_spectrum(self):
pass
-
+
def get_spectrum(self):
pass
-
+
class XSpecThermalModel(SpectralModel):
r"""
Initialize a thermal gas emission model from PyXspec.
-
+
Parameters
----------
@@ -54,19 +54,25 @@
The maximum energy for the spectral model.
nchan : integer
The number of channels in the spectral model.
+ settings : dictionary, optional
+ A dictionary of key, value pairs (must both be strings)
+ that can be used to set various options in XSPEC.
Examples
--------
>>> mekal_model = XSpecThermalModel("mekal", 0.05, 50.0, 1000)
"""
- def __init__(self, model_name, emin, emax, nchan, thermal_broad=False):
+ def __init__(self, model_name, emin, emax, nchan,
+ thermal_broad=False, settings=None):
self.model_name = model_name
self.thermal_broad = thermal_broad
- SpectralModel.__init__(self, emin, emax, nchan)
-
- def prepare(self):
+ if settings is None: settings = {}
+ self.settings = settings
+ super(XSpecThermalModel, self).__init__(emin, emax, nchan)
+
+ def prepare_spectrum(self, zobs):
"""
- Prepare the thermal model for execution.
+ Prepare the thermal model for execution given a redshift *zobs* for the spectrum.
"""
import xspec
xspec.Xset.chatter = 0
@@ -79,9 +85,11 @@
else:
self.norm = 1.0e-14
self.thermal_comp.norm = 1.0
- self.thermal_comp.Redshift = 0.0
+ self.thermal_comp.Redshift = zobs
if self.thermal_broad:
xspec.Xset.addModelString("APECTHERMAL","yes")
+ for k,v in self.settings.items():
+ xspec.Xset.addModelString(k,v)
def get_spectrum(self, kT):
"""
@@ -89,18 +97,20 @@
"""
self.thermal_comp.kT = kT
self.thermal_comp.Abundanc = 0.0
- cosmic_spec = self.norm*np.array(self.model.values(0))
+ cosmic_spec = np.array(self.model.values(0))
if self.model_name == "bremss":
metal_spec = np.zeros(self.nchan)
else:
self.thermal_comp.Abundanc = 1.0
- metal_spec = self.norm*np.array(self.model.values(0)) - cosmic_spec
+ metal_spec = np.array(self.model.values(0)) - cosmic_spec
+ cosmic_spec *= self.norm
+ metal_spec *= self.norm
return YTArray(cosmic_spec, "cm**3/s"), YTArray(metal_spec, "cm**3/s")
-
+
class XSpecAbsorbModel(SpectralModel):
r"""
Initialize an absorption model from PyXspec.
-
+
Parameters
----------
@@ -114,19 +124,25 @@
The maximum energy for the spectral model.
nchan : integer, optional
The number of channels in the spectral model.
+ settings : dictionary, optional
+ A dictionary of key, value pairs (must both be strings)
+ that can be used to set various options in XSPEC.
Examples
--------
>>> abs_model = XSpecAbsorbModel("wabs", 0.1)
"""
- def __init__(self, model_name, nH, emin=0.01, emax=50.0, nchan=100000):
+ def __init__(self, model_name, nH, emin=0.01, emax=50.0,
+ nchan=100000, settings=None):
self.model_name = model_name
self.nH = nH
- SpectralModel.__init__(self, emin, emax, nchan)
-
- def prepare(self):
+ if settings is None: settings = {}
+ self.settings = settings
+ super(XSpecAbsorbModel, self).__init__(emin, emax, nchan)
+
+ def prepare_spectrum(self):
"""
- Prepare the absorption model for execution.
+ Prepare the absorption model for execution given a redshift *zobs* for the spectrum.
"""
import xspec
xspec.Xset.chatter = 0
@@ -135,6 +151,8 @@
self.model = xspec.Model(self.model_name+"*powerlaw")
self.model.powerlaw.norm = self.nchan/(self.emax.value-self.emin.value)
self.model.powerlaw.PhoIndex = 0.0
+ for k,v in self.settings.items():
+ xspec.Xset.addModelString(k,v)
def get_spectrum(self):
"""
@@ -150,7 +168,7 @@
available at http://www.atomdb.org. This code borrows heavily from Python
routines used to read the APEC tables developed by Adam Foster at the
CfA (afoster at cfa.harvard.edu).
-
+
Parameters
----------
@@ -183,7 +201,7 @@
self.apec_prefix+"_coco.fits")
self.linefile = os.path.join(self.apec_root,
self.apec_prefix+"_line.fits")
- SpectralModel.__init__(self, emin, emax, nchan)
+ super(TableApecModel, self).__init__(emin, emax, nchan)
self.wvbins = hc/self.ebins[::-1].d
# H, He, and trace elements
self.cosmic_elem = [1,2,3,4,5,9,11,15,17,19,21,22,23,24,25,27,29,30]
@@ -196,8 +214,8 @@
32.0650,35.4530,39.9480,39.0983,40.0780,
44.9559,47.8670,50.9415,51.9961,54.9380,
55.8450,58.9332,58.6934,63.5460,65.3800])
-
- def prepare(self):
+
+ def prepare_spectrum(self, zobs):
"""
Prepare the thermal model for execution.
"""
@@ -216,17 +234,18 @@
self.dTvals = np.diff(self.Tvals)
self.minlam = self.wvbins.min()
self.maxlam = self.wvbins.max()
-
+ self.scale_factor = 1.0/(1.+zobs)
+
def _make_spectrum(self, element, tindex):
-
+
tmpspec = np.zeros(self.nchan)
-
+
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]
vec = np.zeros(self.nchan)
- E0 = hc/self.line_handle[tindex].data.field('lambda')[i]
+ E0 = hc/self.line_handle[tindex].data.field('lambda')[i]*self.scale_factor
amp = self.line_handle[tindex].data.field('epsilon')[i]
ebins = self.ebins.d
if self.thermal_broad:
@@ -246,19 +265,19 @@
return tmpspec
else:
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]
+ e_cont = self.coco_handle[tindex].data.field('E_Cont')[ind][:n_cont]*self.scale_factor
continuum = self.coco_handle[tindex].data.field('Continuum')[ind][:n_cont]
tmpspec += np.interp(self.emid.d, e_cont, continuum)*self.de.d
-
+
n_pseudo = self.coco_handle[tindex].data.field('N_Pseudo')[ind]
- e_pseudo = self.coco_handle[tindex].data.field('E_Pseudo')[ind][:n_pseudo]
+ e_pseudo = self.coco_handle[tindex].data.field('E_Pseudo')[ind][:n_pseudo]*self.scale_factor
pseudo = self.coco_handle[tindex].data.field('Pseudo')[ind][:n_pseudo]
-
+
tmpspec += np.interp(self.emid.d, e_pseudo, pseudo)*self.de.d
-
+
return tmpspec
def get_spectrum(self, kT):
@@ -288,7 +307,7 @@
class TableAbsorbModel(SpectralModel):
r"""
Initialize an absorption model from a table stored in an HDF5 file.
-
+
Parameters
----------
@@ -299,7 +318,7 @@
Examples
--------
- >>> abs_model = XSpecAbsorbModel("wabs", 0.1)
+ >>> abs_model = XSpecAbsorbModel("abs_table.h5", 0.1)
"""
def __init__(self, filename, nH):
@@ -312,15 +331,15 @@
self.sigma = YTArray(f["cross_section"][:], "cm**2")
nchan = self.sigma.shape[0]
f.close()
- SpectralModel.__init__(self, emin, emax, nchan)
+ super(TableAbsorbModel, self).__init__(emin, emax, nchan)
self.nH = YTQuantity(nH*1.0e22, "cm**-2")
-
- def prepare(self):
+
+ def prepare_spectrum(self):
"""
Prepare the absorption model for execution.
"""
pass
-
+
def get_spectrum(self):
"""
Get the absorption spectrum.
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