[yt-svn] commit/yt: 10 new changesets

commits-noreply at bitbucket.org commits-noreply at bitbucket.org
Mon Jun 1 08:38:17 PDT 2015


10 new commits in yt:

https://bitbucket.org/yt_analysis/yt/commits/5e319f9fa9d6/
Changeset:   5e319f9fa9d6
Branch:      yt
User:        brittonsmith
Date:        2015-05-20 21:51:30+00:00
Summary:     Removing hard-coded unit conversions from absorption spectrum generator.
Affected #:  2 files

diff -r 4cba8216c42868a394e99e3b10cedd31bec9e704 -r 5e319f9fa9d67997fd9c1ec0fe2ecbc293d071b4 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
@@ -171,26 +171,28 @@
 
     ## shift lam0 by deltav
     if deltav is not None:
-        lam1 = lam0 * (1 + deltav / c)
+        lam1 = lam0 * (1 + deltav / ccgs)
     elif delta_lambda is not None:
         lam1 = lam0 + delta_lambda
     else:
         lam1 = lam0
 
     ## 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
+    vdop = vkms.in_units("cm/s")      # in cm/s
+    lam0cgs = lam0.in_units("cm")   # rest wavelength in cm
+    lam1cgs = lam1.in_units("cm")   # line wavelength in cm
+    nu1 = (ccgs / lam1cgs).in_units("1/s")              # line freq in Hz
+    nudop = (vdop / ccgs * nu1).in_units("1/s")         # doppler width in Hz
+    lamdop = (vdop / ccgs * lam1).in_units("angstrom")       # doppler width in Ang
 
+    #import pdb ; pdb.set_trace()
+        
     ## 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)
+    nua = (ccgs / lambda_bins).in_units("1/s") # frequency vector (Hz)
 
     ## tau_0
     tau_X = np.sqrt(np.pi) * e**2 / (me * ccgs) * \

diff -r 4cba8216c42868a394e99e3b10cedd31bec9e704 -r 5e319f9fa9d67997fd9c1ec0fe2ecbc293d071b4 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
+    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,7 +77,8 @@
            mass of atom in amu.
         """
         self.line_list.append({'label': label, 'field_name': field_name,
-                               'wavelength': wavelength, 'f_value': f_value,
+                               'wavelength': YTQuantity(wavelength, "angstrom"),
+                               'f_value': f_value,
                                'gamma': gamma, 'atomic_mass': atomic_mass,
                                'label_threshold': label_threshold})
 
@@ -209,16 +211,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']) /
+                                  (amu_cgs * 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 +241,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 +263,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 +274,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 +283,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 +298,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 +310,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 +322,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)


https://bitbucket.org/yt_analysis/yt/commits/287d5d45255b/
Changeset:   287d5d45255b
Branch:      yt
User:        brittonsmith
Date:        2015-05-21 17:48:17+00:00
Summary:     Cleaning up voigt profile generator and making it use units more properly.
Affected #:  1 file

diff -r 5e319f9fa9d67997fd9c1ec0fe2ecbc293d071b4 -r 287d5d45255bc8fd3a2ea707e5cdb8dd89daeb2f 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,69 +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 / ccgs)
+    ## 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.in_units("cm/s")      # in cm/s
-    lam0cgs = lam0.in_units("cm")   # rest wavelength in cm
-    lam1cgs = lam1.in_units("cm")   # line wavelength in cm
-    nu1 = (ccgs / lam1cgs).in_units("1/s")              # line freq in Hz
-    nudop = (vdop / ccgs * nu1).in_units("1/s")         # doppler width in Hz
-    lamdop = (vdop / ccgs * lam1).in_units("angstrom")       # 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
 
-    #import pdb ; pdb.set_trace()
-        
     ## 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).in_units("1/s") # 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)


https://bitbucket.org/yt_analysis/yt/commits/400801f9a26b/
Changeset:   400801f9a26b
Branch:      yt
User:        brittonsmith
Date:        2015-05-21 18:13:04+00:00
Summary:     Fixing a tab.
Affected #:  1 file

diff -r 287d5d45255bc8fd3a2ea707e5cdb8dd89daeb2f -r 400801f9a26b821b668e99c355503527eaabeb8f 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
@@ -193,7 +193,7 @@
                                      column_density[lixel] / continuum['normalization']
                 self.tau_field[left_index[lixel]:right_index[lixel]] += line_tau
                 pbar.update(i)
-            pbar.finish()
+        pbar.finish()
 
     def _add_lines_to_spectrum(self, field_data, use_peculiar_velocity):
         """


https://bitbucket.org/yt_analysis/yt/commits/515a01ad24f1/
Changeset:   515a01ad24f1
Branch:      yt
User:        brittonsmith
Date:        2015-05-21 19:20:37+00:00
Summary:     It was right all along.
Affected #:  1 file

diff -r 400801f9a26b821b668e99c355503527eaabeb8f -r 515a01ad24f13f22619c3226c8f7d72e2d1a945e 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
@@ -193,7 +193,7 @@
                                      column_density[lixel] / continuum['normalization']
                 self.tau_field[left_index[lixel]:right_index[lixel]] += line_tau
                 pbar.update(i)
-        pbar.finish()
+            pbar.finish()
 
     def _add_lines_to_spectrum(self, field_data, use_peculiar_velocity):
         """


https://bitbucket.org/yt_analysis/yt/commits/2276a782446a/
Changeset:   2276a782446a
Branch:      yt
User:        brittonsmith
Date:        2015-05-24 15:56:09+00:00
Summary:     Atomic masses now are YTQuantities.
Affected #:  1 file

diff -r 515a01ad24f13f22619c3226c8f7d72e2d1a945e -r 2276a782446a539fc9f925f9ef97cd8a87f58156 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
@@ -23,7 +23,7 @@
 from yt.funcs import get_pbar, mylog
 from yt.units.yt_array import YTArray, YTQuantity
 from yt.utilities.physical_constants import \
-    amu_cgs, boltzmann_constant_cgs, \
+    boltzmann_constant_cgs, \
     speed_of_light_cgs
 from yt.utilities.on_demand_imports import _astropy
 
@@ -79,7 +79,8 @@
         self.line_list.append({'label': label, 'field_name': field_name,
                                'wavelength': YTQuantity(wavelength, "angstrom"),
                                'f_value': f_value,
-                               'gamma': gamma, 'atomic_mass': atomic_mass,
+                               'gamma': gamma,
+                               'atomic_mass': YTQuantity(atomic_mass, "amu"),
                                'label_threshold': label_threshold})
 
     def add_continuum(self, label, field_name, wavelength,
@@ -213,7 +214,7 @@
                     field_data['velocity_los'] / speed_of_light_cgs
             thermal_b =  np.sqrt((2 * boltzmann_constant_cgs *
                                   field_data['temperature']) /
-                                  (amu_cgs * line['atomic_mass']))
+                                  line['atomic_mass'])
             center_bins = np.digitize((delta_lambda + line['wavelength']),
                                       self.lambda_bins)
 


https://bitbucket.org/yt_analysis/yt/commits/d26b1b67dee3/
Changeset:   d26b1b67dee3
Branch:      yt
User:        brittonsmith
Date:        2015-05-24 19:55:26+00:00
Summary:     Updating cookbook recipe.
Affected #:  1 file

diff -r 2276a782446a539fc9f925f9ef97cd8a87f58156 -r d26b1b67dee311dfe1b6b4f4384d1eb036c03ea3 doc/source/cookbook/fit_spectrum.py
--- a/doc/source/cookbook/fit_spectrum.py
+++ b/doc/source/cookbook/fit_spectrum.py
@@ -100,6 +100,6 @@
 order_fits = ['OVI', 'HI']
 
 # Fit spectrum and save fit
-fitted_lines, fitted_flux = generate_total_fit(wavelength,
+fitted_lines, fitted_flux = generate_total_fit(wavelength.d,
                                                flux, order_fits, species_dicts,
                                                output_file='spectrum_fit.h5')


https://bitbucket.org/yt_analysis/yt/commits/5e426c5aaf38/
Changeset:   5e426c5aaf38
Branch:      yt
User:        brittonsmith
Date:        2015-05-29 19:46:17+00:00
Summary:     Undoing change to cookbook.
Affected #:  1 file

diff -r d26b1b67dee311dfe1b6b4f4384d1eb036c03ea3 -r 5e426c5aaf389407192b07f07f0a7cb0d9cbabba doc/source/cookbook/fit_spectrum.py
--- a/doc/source/cookbook/fit_spectrum.py
+++ b/doc/source/cookbook/fit_spectrum.py
@@ -100,6 +100,6 @@
 order_fits = ['OVI', 'HI']
 
 # Fit spectrum and save fit
-fitted_lines, fitted_flux = generate_total_fit(wavelength.d,
+fitted_lines, fitted_flux = generate_total_fit(wavelength,
                                                flux, order_fits, species_dicts,
                                                output_file='spectrum_fit.h5')


https://bitbucket.org/yt_analysis/yt/commits/ee1cb9fa80bd/
Changeset:   ee1cb9fa80bd
Branch:      yt
User:        brittonsmith
Date:        2015-05-29 19:50:48+00:00
Summary:     Stripping units on wavelength array if they're there.
Affected #:  1 file

diff -r 5e426c5aaf389407192b07f07f0a7cb0d9cbabba -r ee1cb9fa80bd207896d77470a1365469521a7b56 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
@@ -79,6 +79,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 = {}
 


https://bitbucket.org/yt_analysis/yt/commits/691da9acc759/
Changeset:   691da9acc759
Branch:      yt
User:        brittonsmith
Date:        2015-05-29 20:49:05+00:00
Summary:     Importing YTArray and cleaning up imports.
Affected #:  1 file

diff -r ee1cb9fa80bd207896d77470a1365469521a7b56 -r 691da9acc75937e1ba62aa3620c7858130d82d7c 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
 
@@ -1011,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()
-


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