1. Britton Smith
  2. light_ray_tools

Commits

Britton Smith  committed 2e7a4de

Adding option to get b values for lines.

  • Participants
  • Parent commits 6bfebf3
  • Branches default

Comments (0)

Files changed (3)

File get_dndz.py

View file
-from lightraytools import *
    files = glob("%s/*.h5" % data_dir)
    # Resolution here is the spectral resolution.  This will determine 
    # the number of ray elements blended together to form one absorber.
    # Range is the range in log column densities to be included. 
    # Format for range is, for example, range=[14, 20].
    # If range is set to None, as it is now, all absorbers are included.
    dn, dn_bins, total_rays, total_path, field_stats = dndl_rays(files, field, bins=bins, 
                                                                 path='dz', resolution=5000,
                                                                 mean_fields=mean_fields,
                                                                 variance_fields=variance_fields,
                                                                 filename=filename)
    output_file = "%s.h5" % field
    # This is the field in the ray data that has the ion fraction.
    fraction_field = "%sIon_Fraction" % (field[0:field.find('NumberDensity')])
    # Call the above wrapper for getting dn/dz.
    # stats_fields is a list of fields for which we want the mean properties for each column density bin.
    get_dndz(data_dir, field, bins=bins, 
             mean_fields=['Metallicity', 'Temperature', 'Density', fraction_field],
             variance_fields=['los_velocity'],
             filename=output_file)
+from light_ray_tools import *
    files = glob("%s/*.h5" % data_dir)
    # Resolution here is the spectral resolution.  This will determine 
    # the number of ray elements blended together to form one absorber.
    # Range is the range in log column densities to be included. 
    # Format for range is, for example, range=[14, 20].
    # If range is set to None, as it is now, all absorbers are included.
    dn, dn_bins, total_rays, total_path, field_stats = dndl_rays(files, field, bins=bins, 
                                                                 path='dz', resolution=5000,
                                                                 mean_fields=mean_fields,
                                                                 variance_fields=variance_fields,
                                                                 filename=filename)
    output_file = "%s.h5" % field
    # This is the field in the ray data that has the ion fraction.
    fraction_field = "%sIon_Fraction" % (field[0:field.find('NumberDensity')])
    # Call the above wrapper for getting dn/dz.
    # stats_fields is a list of fields for which we want the mean properties for each column density bin.
    get_dndz(data_dir, field, bins=bins, 
             mean_fields=['Metallicity', 'Temperature', 'Density', fraction_field],
             variance_fields=['los_velocity'], get_b_values=b_constant,
             filename=output_file)

File ray_dndl.py

View file
 
 def dndl_rays(lightray_files, field, path='dz', bins=None, 
               mean_fields=None, variance_fields=None,
-              filename=None, **kwargs):
+              filename=None, get_b_values=False, **kwargs):
     "Calculate dn/dl for multiple light rays."
 
     dn_total = na.zeros((bins.size-1), dtype=float)
     for file in lightray_files:
         dn, dn_bins, dpath, f_stats = dndl_ray(file, field, path=path, bins=bins, 
                                                mean_fields=mean_fields, variance_fields=variance_fields,
-                                               **kwargs)
+                                               get_b_values=get_b_values, **kwargs)
         bins = dn_bins
 
-        if field_stats is None:
-            field_stats = {}
-            for s_field in mean_fields + variance_fields + [field]:
-                field_stats[s_field] = [[] for bin in range(dn_bins.size - 1)]
-
+        if get_b_values:
+            mean_fields.append('b_value')
+        field_stats = {}
         for s_field in mean_fields + variance_fields + [field]:
+            field_stats[s_field] = [[] for bin in range(dn_bins.size - 1)]
             for bin in range(dn_bins.size - 1):
                 field_stats[s_field][bin].extend(f_stats[s_field][bin])
 
     for s_field in mean_fields + variance_fields + [field]:
         final_field_stats[s_field] = na.zeros(((dn_bins.size - 1),2))
         for bin in range(dn_bins.size - 1):
-            if s_field == field:
+            if s_field == field or s_field == 'b_value':
                 final_field_stats[s_field][bin][0] = na.array(field_stats[s_field][bin]).mean()
                 final_field_stats[s_field][bin][1] = na.array(field_stats[s_field][bin]).std()
             else:
     return (dn, dn_bins_total, len(lightray_files), total_path, final_field_stats)
 
 def dndl_ray(lightray_file, field, path='dz', resolution=None, cosmology=None, 
-             mean_fields=None, variance_fields=None, **kwargs):
+             mean_fields=None, variance_fields=None, get_b_values=False, **kwargs):
     "Calculte dn/dl for a single ray."
 
     field_stats = None
         ray_data[field] *= ray_data['dl']
     else:        
         new_data = ray_resample(resolution, ray_data, fields=[field], weight_field=field, 
-                                mean_fields=mean_fields, variance_fields=variance_fields)
+                                mean_fields=mean_fields, variance_fields=variance_fields,
+                                get_b_values=get_b_values)
         ray_data.update(new_data)
         del new_data
     del ray_data['dl']
     dn, dn_bins = na.histogram(ray_data[field], **kwargs)
 
     field_stats = {field: []}
+    if get_b_values:
+        field_stats['b_value'] = []
     for s_field in mean_fields + variance_fields:
         field_stats[s_field] = []
     digits = na.digitize(ray_data[field], dn_bins)
         field_stats[field].append(ray_data[field][indices])
         for s_field in mean_fields + variance_fields:
             field_stats[s_field].append(ray_data[s_field][indices])
+        if get_b_values:
+            field_stats['b_value'].append(ray_data['b_value'][indices])
 
     print "%s - %s: path (%s): %.2f, min = %.2f, max = %.2f." % \
         (lightray_file, field, path, dpath, ray_data[field].min(), ray_data[field].max())

File ray_resample.py

View file
 
 def ray_resample(resolution, ray_data, z_field='z', dz_field='dz', dl_field='dl', 
                  fields=None, mean_fields=None, variance_fields=None,
-                 weight_field=None, distance_field=None):
+                 weight_field=None, distance_field=None, get_b_values=False):
     """
     Resample light ray at a constant resolution.
     dz = (z + 1) / R
         new_data[field] = na.zeros(new_data[z_field].size, dtype=new_data[z_field].dtype)
     if distance_field is not None:
         new_data[distance_field] = na.zeros(new_data[z_field].size, dtype=new_data[z_field].dtype)
+    if get_b_values:
+        new_data['b_value'] = na.zeros(new_data[z_field].size, dtype=new_data[z_field].dtype)
 
     for new_index in xrange(new_data[z_field].size):
         for field in fields:
             new_data[distance_field][new_index] = distance_field.take(map(lambda a: a[0], 
                                                                           index_list[new_index])).min()
 
+        if get_b_values:
+            print "Thermal %e." % (na.sqrt(get_b_values * new_data['Temperature'][new_index]))
+            print "Vel %e." % new_data['los_velocity'][new_index]
+            new_data['b_value'][new_index] = 1.e-5 * na.sqrt(get_b_values * new_data['Temperature'][new_index] +
+                                                             new_data['los_velocity'][new_index]**2)
+            print new_data['b_value'][new_index]
+
     return new_data