[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