[yt-svn] commit/yt: xarthisius: Merged in brittonsmith/yt (pull request #1596)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon Jun 1 08:38:18 PDT 2015
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/1dee4eb06343/
Changeset: 1dee4eb06343
Branch: yt
User: xarthisius
Date: 2015-06-01 15:38:09+00:00
Summary: Merged in brittonsmith/yt (pull request #1596)
Fixing unit support in AbsorptionSpectrum
Affected #: 4 files
diff -r d6983e51e584eb0e996dfd2a4475ff2a8e0a6014 -r 1dee4eb0634354c90ace1ba14a6d2b4ed4237882 yt/analysis_modules/absorption_spectrum/absorption_line.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_line.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_line.py
@@ -16,12 +16,9 @@
import numpy as np
from yt.utilities.physical_constants import \
charge_proton_cgs, \
- cm_per_km, \
- km_per_cm, \
mass_electron_cgs, \
speed_of_light_cgs
-
def voigt(a,u):
"""
NAME:
@@ -140,67 +137,75 @@
k1 = k1.astype(np.float64).clip(0)
return k1
-def tau_profile(lam0, fval, gamma, vkms, column_density,
- deltav=None, delta_lambda=None,
+def tau_profile(lamba_0, f_value, gamma, v_doppler, column_density,
+ delta_v=None, delta_lambda=None,
lambda_bins=None, n_lambda=12000, dlambda=0.01):
- """
+ r"""
Create an optical depth vs. wavelength profile for an
absorption line using a voigt profile.
- :param lam0 (float): central wavelength (angstroms).
- :param fval (float): f-value.
- :param gamma (float): gamma value.
- :param vkms (float): doppler b-parameter.
- :param column_density (float): column density (cm^-2).
- :param deltav (float): velocity offset from lam0 (km/s).
- Default: None (no shift).
- :param delta_lambda (float): wavelength offset in angstroms.
- Default: None (no shift).
- :param lambda_bins (array): array of wavelengths in angstroms.
- Default: None
- :param n_lambda (float): size of lambda bins to create
- array if lambda_bins is None. Default: 12000
- :param dlambda (float): lambda bin width if lambda_bins is
- None. Default: 0.01
+
+ Parameters
+ ----------
+
+ lamba_0 : float YTQuantity in length units
+ central wavelength.
+ f_value : float
+ absorption line f-value.
+ gamma : float
+ absorption line gamma value.
+ v_doppler : float YTQuantity in velocity units
+ doppler b-parameter.
+ column_density : float YTQuantity in (length units)^-2
+ column density.
+ delta_v : float YTQuantity in velocity units
+ velocity offset from lamba_0.
+ Default: None (no shift).
+ delta_lambda : float YTQuantity in length units
+ wavelength offset.
+ Default: None (no shift).
+ lambda_bins : YTArray in length units
+ wavelength array for line deposition. If None, one will be
+ created using n_lambda and dlambda.
+ Default: None.
+ n_lambda : int
+ size of lambda bins to create if lambda_bins is None.
+ Default: 12000.
+ dlambda : float
+ lambda bin width in angstroms if lambda_bins is None.
+ Default: 0.01.
+
"""
- ## constants
- me = mass_electron_cgs # grams mass electron
- e = charge_proton_cgs # esu
- c = speed_of_light_cgs * km_per_cm # km/s
- ccgs = speed_of_light_cgs # cm/s
-
- ## shift lam0 by deltav
- if deltav is not None:
- lam1 = lam0 * (1 + deltav / c)
+ ## shift lamba_0 by delta_v
+ if delta_v is not None:
+ lam1 = lamba_0 * (1 + delta_v / speed_of_light_cgs)
elif delta_lambda is not None:
- lam1 = lam0 + delta_lambda
+ lam1 = lamba_0 + delta_lambda
else:
- lam1 = lam0
+ lam1 = lamba_0
## conversions
- vdop = vkms * cm_per_km # in cm/s
- lam0cgs = lam0 / 1.e8 # rest wavelength in cm
- lam1cgs = lam1 / 1.e8 # line wavelength in cm
- nu1 = ccgs / lam1cgs # line freq in Hz
- nudop = vdop / ccgs * nu1 # doppler width in Hz
- lamdop = vdop / ccgs * lam1 # doppler width in Ang
+ nu1 = speed_of_light_cgs / lam1 # line freq in Hz
+ nudop = v_doppler / speed_of_light_cgs * nu1 # doppler width in Hz
+ lamdop = v_doppler / speed_of_light_cgs * lam1 # doppler width in Ang
## create wavelength
if lambda_bins is None:
lambda_bins = lam1 + \
np.arange(n_lambda, dtype=np.float) * dlambda - \
- n_lambda * dlambda / 2 # wavelength vector (angstroms)
- nua = ccgs / (lambda_bins / 1.e8) # frequency vector (Hz)
+ n_lambda * dlambda / 2 # wavelength vector (angstroms)
+ nua = (speed_of_light_cgs / lambda_bins) # frequency vector (Hz)
## tau_0
- tau_X = np.sqrt(np.pi) * e**2 / (me * ccgs) * \
- column_density * fval / vdop
- tau0 = tau_X * lam0cgs
+ tau_X = np.sqrt(np.pi) * charge_proton_cgs**2 / \
+ (mass_electron_cgs * speed_of_light_cgs) * \
+ column_density * f_value / v_doppler
+ tau0 = tau_X * lamba_0
# dimensionless frequency offset in units of doppler freq
- x = (nua - nu1) / nudop
- a = gamma / (4 * np.pi * nudop) # damping parameter
- phi = voigt(a, x) # profile
- tauphi = tau0 * phi # profile scaled with tau0
+ x = ((nua - nu1) / nudop).in_units("")
+ a = (gamma / (4 * np.pi * nudop)).in_units("s") # damping parameter
+ phi = voigt(a, x) # line profile
+ tauphi = (tau0 * phi).in_units("") # profile scaled with tau0
return (lambda_bins, tauphi)
diff -r d6983e51e584eb0e996dfd2a4475ff2a8e0a6014 -r 1dee4eb0634354c90ace1ba14a6d2b4ed4237882 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
@@ -4,7 +4,7 @@
"""
-from __future__ import print_function
+
from __future__ import absolute_import
#-----------------------------------------------------------------------------
@@ -21,13 +21,12 @@
from .absorption_line import tau_profile
from yt.funcs import get_pbar, mylog
-from yt.units.yt_array import YTArray
+from yt.units.yt_array import YTArray, YTQuantity
from yt.utilities.physical_constants import \
- amu_cgs, boltzmann_constant_cgs, \
- speed_of_light_cgs, km_per_cm
+ boltzmann_constant_cgs, \
+ speed_of_light_cgs
from yt.utilities.on_demand_imports import _astropy
-speed_of_light_kms = speed_of_light_cgs * km_per_cm
pyfits = _astropy.pyfits
class AbsorptionSpectrum(object):
@@ -49,8 +48,10 @@
self.tau_field = None
self.flux_field = None
self.spectrum_line_list = None
- self.lambda_bins = np.linspace(lambda_min, lambda_max, n_lambda)
- self.bin_width = (lambda_max - lambda_min) / float(n_lambda - 1)
+ self.lambda_bins = YTArray(np.linspace(lambda_min, lambda_max, n_lambda),
+ "angstrom")
+ self.bin_width = YTQuantity((lambda_max - lambda_min) /
+ float(n_lambda - 1), "angstrom")
self.line_list = []
self.continuum_list = []
@@ -76,8 +77,10 @@
mass of atom in amu.
"""
self.line_list.append({'label': label, 'field_name': field_name,
- 'wavelength': wavelength, 'f_value': f_value,
- 'gamma': gamma, 'atomic_mass': atomic_mass,
+ 'wavelength': YTQuantity(wavelength, "angstrom"),
+ 'f_value': f_value,
+ 'gamma': gamma,
+ 'atomic_mass': YTQuantity(atomic_mass, "amu"),
'label_threshold': label_threshold})
def add_continuum(self, label, field_name, wavelength,
@@ -209,16 +212,15 @@
# include factor of (1 + z) because our velocity is in proper frame.
delta_lambda += line['wavelength'] * (1 + field_data['redshift']) * \
field_data['velocity_los'] / speed_of_light_cgs
- thermal_b = km_per_cm * np.sqrt((2 * boltzmann_constant_cgs *
- field_data['temperature']) /
- (amu_cgs * line['atomic_mass']))
- thermal_b.convert_to_cgs()
+ thermal_b = np.sqrt((2 * boltzmann_constant_cgs *
+ field_data['temperature']) /
+ line['atomic_mass'])
center_bins = np.digitize((delta_lambda + line['wavelength']),
self.lambda_bins)
# ratio of line width to bin width
width_ratio = ((line['wavelength'] + delta_lambda) * \
- thermal_b / speed_of_light_kms / self.bin_width).value
+ thermal_b / speed_of_light_cgs / self.bin_width).in_units("").d
if (width_ratio < 1.0).any():
mylog.warn(("%d out of %d line components are unresolved, " +
@@ -240,11 +242,13 @@
my_bin_ratio = spectrum_bin_ratio
while True:
lambda_bins, line_tau = \
- tau_profile(line['wavelength'], line['f_value'],
- line['gamma'], thermal_b[lixel],
- column_density[lixel],
- delta_lambda=delta_lambda[lixel],
- lambda_bins=self.lambda_bins[left_index[lixel]:right_index[lixel]])
+ tau_profile(
+ line['wavelength'], line['f_value'],
+ line['gamma'], thermal_b[lixel].in_units("km/s"),
+ column_density[lixel],
+ delta_lambda=delta_lambda[lixel],
+ lambda_bins=self.lambda_bins[left_index[lixel]:right_index[lixel]])
+
# Widen wavelength window until optical depth reaches a max value at the ends.
if (line_tau[0] < max_tau and line_tau[-1] < max_tau) or \
(left_index[lixel] <= 0 and right_index[lixel] >= self.n_lambda):
@@ -260,7 +264,7 @@
if line['label_threshold'] is not None and \
column_density[lixel] >= line['label_threshold']:
if use_peculiar_velocity:
- peculiar_velocity = km_per_cm * field_data['velocity_los'][lixel]
+ peculiar_velocity = field_data['velocity_los'][lixel].in_units("km/s")
else:
peculiar_velocity = 0.0
self.spectrum_line_list.append({'label': line['label'],
@@ -271,7 +275,7 @@
'redshift': field_data['redshift'][lixel],
'v_pec': peculiar_velocity})
pbar.update(i)
- pbar.finish()
+ pbar.finish()
del column_density, delta_lambda, thermal_b, \
center_bins, width_ratio, left_index, right_index
@@ -280,7 +284,7 @@
"""
Write out list of spectral lines.
"""
- print("Writing spectral line list: %s." % filename)
+ mylog.info("Writing spectral line list: %s." % filename)
self.spectrum_line_list.sort(key=lambda obj: obj['wavelength'])
f = open(filename, 'w')
f.write('#%-14s %-14s %-12s %-12s %-12s %-12s\n' %
@@ -295,7 +299,7 @@
"""
Write spectrum to an ascii file.
"""
- print("Writing spectrum to ascii file: %s." % filename)
+ mylog.info("Writing spectrum to ascii file: %s." % filename)
f = open(filename, 'w')
f.write("# wavelength[A] tau flux\n")
for i in range(self.lambda_bins.size):
@@ -307,7 +311,7 @@
"""
Write spectrum to a fits file.
"""
- print("Writing spectrum to fits file: %s." % filename)
+ mylog.info("Writing spectrum to fits file: %s." % filename)
col1 = pyfits.Column(name='wavelength', format='E', array=self.lambda_bins)
col2 = pyfits.Column(name='flux', format='E', array=self.flux_field)
cols = pyfits.ColDefs([col1, col2])
@@ -319,7 +323,7 @@
Write spectrum to an hdf5 file.
"""
- print("Writing spectrum to hdf5 file: %s." % filename)
+ mylog.info("Writing spectrum to hdf5 file: %s." % filename)
output = h5py.File(filename, 'w')
output.create_dataset('wavelength', data=self.lambda_bins)
output.create_dataset('tau', data=self.tau_field)
diff -r d6983e51e584eb0e996dfd2a4475ff2a8e0a6014 -r 1dee4eb0634354c90ace1ba14a6d2b4ed4237882 yt/analysis_modules/absorption_spectrum/absorption_spectrum_fit.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_spectrum_fit.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_spectrum_fit.py
@@ -1,9 +1,14 @@
-from __future__ import print_function
+import h5py
import numpy as np
-import h5py
-from yt.analysis_modules.absorption_spectrum.absorption_line \
- import voigt
-from yt.utilities.on_demand_imports import _scipy
+
+from yt.analysis_modules.absorption_spectrum.absorption_line import \
+ voigt
+from yt.funcs import \
+ mylog
+from yt.units.yt_array import \
+ YTArray
+from yt.utilities.on_demand_imports import \
+ _scipy
optimize = _scipy.optimize
@@ -79,6 +84,10 @@
absorption profiles. Same size as x.
"""
+ # convert to NumPy array if we have a YTArray
+ if isinstance(x, YTArray):
+ x = x.d
+
#Empty dictionary for fitted lines
allSpeciesLines = {}
@@ -1007,6 +1016,5 @@
f.create_dataset("{0}/b".format(ion),data=params['b'])
f.create_dataset("{0}/z".format(ion),data=params['z'])
f.create_dataset("{0}/complex".format(ion),data=params['group#'])
- print('Writing spectrum fit to {0}'.format(file_name))
+ mylog.info('Writing spectrum fit to {0}'.format(file_name))
f.close()
-
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