[yt-svn] commit/yt: ngoldbaum: Merged in chummels/yt (pull request #2323)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Thu Aug 18 08:58:51 PDT 2016
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/df6fc87fa819/
Changeset: df6fc87fa819
Branch: yt
User: ngoldbaum
Date: 2016-08-18 15:58:12+00:00
Summary: Merged in chummels/yt (pull request #2323)
AbsorptionSpectrum continuum improvements and fixes
Affected #: 5 files
diff -r 8de8a3c315df8cbf36140c6cc6cd0d4e7a5e748e -r df6fc87fa819404a937c497f8187afa123e92bac doc/source/analyzing/analysis_modules/absorption_spectrum.rst
--- a/doc/source/analyzing/analysis_modules/absorption_spectrum.rst
+++ b/doc/source/analyzing/analysis_modules/absorption_spectrum.rst
@@ -116,7 +116,12 @@
Continuum features with optical depths that follow a power law can also be
added. Like adding lines, you must specify details like the wavelength
and the field in the dataset and LightRay that is tied to this feature.
-Below, we will add H Lyman continuum.
+The wavelength refers to the location at which the continuum begins to be
+applied to the dataset, and as it moves to lower wavelength values, the
+optical depth value decreases according to the defined power law. The
+normalization value is the column density of the linked field which results
+in an optical depth of 1 at the defined wavelength. Below, we add the hydrogen
+Lyman continuum.
.. code-block:: python
@@ -131,7 +136,7 @@
Making the Spectrum
^^^^^^^^^^^^^^^^^^^
-Once all the lines and continuum are added, it is time to make a spectrum out
+Once all the lines and continuua are added, it is time to make a spectrum out
of some light ray data.
.. code-block:: python
diff -r 8de8a3c315df8cbf36140c6cc6cd0d4e7a5e748e -r df6fc87fa819404a937c497f8187afa123e92bac tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -67,11 +67,12 @@
local_ytdata_000:
- yt/frontends/ytdata
- local_absorption_spectrum_001:
+ local_absorption_spectrum_003:
- yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py:test_absorption_spectrum_non_cosmo
- yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py:test_absorption_spectrum_cosmo
- yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py:test_absorption_spectrum_non_cosmo_sph
- yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py:test_absorption_spectrum_cosmo_sph
+ - yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py:test_absorption_spectrum_with_continuum
local_axialpix_001:
- yt/geometry/coordinates/tests/test_axial_pixelization.py:test_axial_pixelization
diff -r 8de8a3c315df8cbf36140c6cc6cd0d4e7a5e748e -r df6fc87fa819404a937c497f8187afa123e92bac yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
@@ -273,42 +273,85 @@
def _add_continua_to_spectrum(self, field_data, use_peculiar_velocity,
observing_redshift=0.):
"""
- Add continuum features to the spectrum.
+ Add continuum features to the spectrum. Continuua are recorded as
+ a name, associated field, wavelength, normalization value, and index.
+ Continuua are applied at and below the denoted wavelength, where the
+ optical depth decreases as a power law of desired index. For positive
+ index values, this means optical depth is highest at the denoted
+ wavelength, and it drops with shorter and shorter wavelengths.
+ Consequently, transmitted flux undergoes a discontinuous cutoff at the
+ denoted wavelength, and then slowly increases with decreasing wavelength
+ according to the power law.
"""
# Change the redshifts of continuum sources to account for the
# redshift at which the observer sits
redshift, redshift_eff = self._apply_observing_redshift(field_data,
use_peculiar_velocity, observing_redshift)
- # Only add continuum features down to tau of 1.e-4.
- min_tau = 1.e-3
+ # min_tau is the minimum optical depth value that warrants
+ # accounting for an absorber. for a single absorber, noticeable
+ # continuum effects begin for tau = 1e-3 (leading to transmitted
+ # flux of e^-tau ~ 0.999). but we apply a cutoff to remove
+ # absorbers with insufficient column_density to contribute
+ # significantly to a continuum (see below). because lots of
+ # low column density absorbers can add up to a significant
+ # continuum effect, we normalize min_tau by the n_absorbers.
+ n_absorbers = field_data['dl'].size
+ min_tau = 1.e-3/n_absorbers
for continuum in self.continuum_list:
- column_density = field_data[continuum['field_name']] * field_data['dl']
+
+ # Normalization is in cm**-2, so column density must be as well
+ column_density = (field_data[continuum['field_name']] *
+ field_data['dl']).in_units('cm**-2')
# redshift_eff field combines cosmological and velocity redshifts
if use_peculiar_velocity:
delta_lambda = continuum['wavelength'] * redshift_eff
else:
delta_lambda = continuum['wavelength'] * redshift
+
+ # right index of continuum affected area is wavelength itself
this_wavelength = delta_lambda + continuum['wavelength']
- right_index = np.digitize(this_wavelength, self.lambda_field).clip(0, self.n_lambda)
+ right_index = np.digitize(this_wavelength,
+ self.lambda_field).clip(0, self.n_lambda)
+ # left index of continuum affected area wavelength at which
+ # optical depth reaches tau_min
left_index = np.digitize((this_wavelength *
- np.power((min_tau * continuum['normalization'] /
- column_density), (1. / continuum['index']))),
- self.lambda_field).clip(0, self.n_lambda)
+ np.power((min_tau * continuum['normalization'] /
+ column_density),
+ (1. / continuum['index']))),
+ self.lambda_field).clip(0, self.n_lambda)
+ # Only calculate the effects of continuua where normalized
+ # column_density is greater than min_tau
+ # because lower column will not have significant contribution
valid_continuua = np.where(((column_density /
continuum['normalization']) > min_tau) &
(right_index - left_index > 1))[0]
+ if valid_continuua.size == 0:
+ mylog.info("Not adding continuum %s: insufficient column density" %
+ continuum['label'])
+ return
+
pbar = get_pbar("Adding continuum - %s [%f A]: " % \
(continuum['label'], continuum['wavelength']),
valid_continuua.size)
+
+ # Tau value is (wavelength / continuum_wavelength)**index /
+ # (column_dens / norm)
+ # i.e. a power law decreasing as wavelength decreases
+
+ # Step through the absorber list and add continuum tau for each to
+ # the total optical depth for all wavelengths
for i, lixel in enumerate(valid_continuua):
- line_tau = np.power((self.lambda_field[left_index[lixel]:right_index[lixel]] /
- this_wavelength[lixel]), continuum['index']) * \
- column_density[lixel] / continuum['normalization']
- self.tau_field[left_index[lixel]:right_index[lixel]] += line_tau
+ cont_tau = \
+ np.power((self.lambda_field[left_index[lixel] :
+ right_index[lixel]] /
+ this_wavelength[lixel]), \
+ continuum['index']) * \
+ (column_density[lixel] / continuum['normalization'])
+ self.tau_field[left_index[lixel]:right_index[lixel]] += cont_tau
pbar.update(i)
pbar.finish()
diff -r 8de8a3c315df8cbf36140c6cc6cd0d4e7a5e748e -r df6fc87fa819404a937c497f8187afa123e92bac yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
@@ -33,7 +33,7 @@
COSMO_PLUS_SINGLE = "enzo_cosmology_plus/RD0009/RD0009"
GIZMO_PLUS = "gizmo_cosmology_plus/N128L16.param"
GIZMO_PLUS_SINGLE = "gizmo_cosmology_plus/snap_N128L16_151.hdf5"
-
+ISO_GALAXY = "IsolatedGalaxy/galaxy0030/galaxy0030"
@requires_file(COSMO_PLUS)
@requires_answer_testing()
@@ -360,3 +360,64 @@
# clean up
os.chdir(curdir)
shutil.rmtree(tmpdir)
+
+ at requires_file(ISO_GALAXY)
+ at requires_answer_testing()
+def test_absorption_spectrum_with_continuum():
+ """
+ This test generates an absorption spectrum from a simple light ray on a
+ grid dataset and adds Lyman alpha and Lyman continuum to it
+ """
+
+ # Set up in a temp dir
+ tmpdir = tempfile.mkdtemp()
+ curdir = os.getcwd()
+ os.chdir(tmpdir)
+
+ ds = load(ISO_GALAXY)
+ lr = LightRay(ds)
+
+ ray_start = ds.domain_left_edge
+ ray_end = ds.domain_right_edge
+ lr.make_light_ray(start_position=ray_start, end_position=ray_end,
+ fields=['temperature', 'density', 'H_number_density'],
+ data_filename='lightray.h5')
+
+ sp = AbsorptionSpectrum(800.0, 1300.0, 5001)
+
+ my_label = 'HI Lya'
+ field = 'H_number_density'
+ wavelength = 1215.6700 # Angstromss
+ f_value = 4.164E-01
+ gamma = 6.265e+08
+ mass = 1.00794
+
+ sp.add_line(my_label, field, wavelength, f_value,
+ gamma, mass, label_threshold=1.e10)
+
+ my_label = 'Ly C'
+ field = 'H_number_density'
+ wavelength = 912.323660 # Angstroms
+ normalization = 1.6e17
+ index = 3.0
+
+ sp.add_continuum(my_label, field, wavelength, normalization, index)
+
+ wavelength, flux = sp.make_spectrum('lightray.h5',
+ output_file='spectrum.h5',
+ line_list_file='lines.txt',
+ use_peculiar_velocity=True)
+
+ # load just-generated hdf5 file of spectral data (for consistency)
+ data = h5.File('spectrum.h5', 'r')
+
+ for key in data.keys():
+ func = lambda x=key: data[x][:]
+ func.__name__ = "{}_continuum".format(key)
+ test = GenericArrayTest(None, func)
+ test_absorption_spectrum_with_continuum.__name__ = test.description
+ yield test
+
+ # clean up
+ os.chdir(curdir)
+ shutil.rmtree(tmpdir)
diff -r 8de8a3c315df8cbf36140c6cc6cd0d4e7a5e748e -r df6fc87fa819404a937c497f8187afa123e92bac yt/utilities/answer_testing/framework.py
--- a/yt/utilities/answer_testing/framework.py
+++ b/yt/utilities/answer_testing/framework.py
@@ -819,6 +819,10 @@
kwargs = self.kwargs
return self.array_func(*args, **kwargs)
def compare(self, new_result, old_result):
+ if not isinstance(new_result, dict):
+ new_result = {'answer': new_result}
+ old_result = {'answer': old_result}
+
assert_equal(len(new_result), len(old_result),
err_msg="Number of outputs not equal.",
verbose=True)
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