Commits

Britton Smith  committed 05e7150

Adding cosmological units to cooling rate script, allowing it to work
properly for redshifts other than 0.

  • Participants
  • Parent commits 24c84c7

Comments (0)

Files changed (1)

File src/python/examples/cooling_rate.py

 import copy
+import numpy as np
 import sys
 
 def check_convergence(fc1, fc2, fields=None, tol=0.10):
             return False
     return True
 
+def set_cosmology_units(my_chemistry, hubble_constant=0.704,
+                        omega_matter=0.268, omega_lambda=0.732,
+                        current_redshift=0.0, initial_redshift=0.0,
+                        comoving_box_size=1.0):
+    my_chemistry.a_units = 1. / (1. + initial_redshift)
+    my_chemistry.density_units = 1.8788e-29 * omega_matter * \
+      np.power(hubble_constant,2) * np.power(1 + current_redshift,3)
+    my_chemistry.length_units = 3.085678e24 * comoving_box_size / \
+      hubble_constant / (1. + current_redshift)
+    my_chemistry.time_units = 2.519445e17 / np.sqrt(omega_matter) / \
+     hubble_constant / np.power(1 + initial_redshift, 1.5)
+    temperature_units = 1.81723e6 * \
+      np.power(comoving_box_size, 2) * omega_matter * (1. + initial_redshift)
+    velocity_units = 1.22475e7 * comoving_box_size * \
+     np.sqrt(omega_matter)* np.sqrt(1. + initial_redshift)
+    return (temperature_units, velocity_units)
+
 kboltz      = 1.3806504e-16
 mass_h      = 1.67262171e-24   
 mass_e      = 9.10938215e-28
 my_chemistry.with_radiative_cooling = 0
 my_chemistry.primordial_chemistry = 3
 my_chemistry.metal_cooling = 1
-my_chemistry.UVbackground = 1;
+my_chemistry.UVbackground = 1
 my_chemistry.grackle_data_file = "CloudyData_UVB=HM2012.h5"
 my_chemistry.include_metal_heating = 1
 
-my_chemistry.comoving_coordinates = 0
-my_chemistry.density_units = mass_h
-my_chemistry.length_units = 1.0e16
-my_chemistry.time_units = yr_to_s * 1e6
-my_chemistry.a_units = 1.0
+my_chemistry.comoving_coordinates = 1
+current_redshift = 0.0
+a_value = 1.0 / (1.0 + current_redshift)
+temperature_units, velocity_units = \
+  set_cosmology_units(my_chemistry, 
+                      current_redshift=current_redshift)
 
-energy_units = (my_chemistry.length_units /
-                my_chemistry.time_units)**2.0
-
-a_value = 1.0
+print my_chemistry.density_units
+print my_chemistry.length_units
+print my_chemistry.a_units
+print my_chemistry.time_units
+print temperature_units
+print velocity_units
 
 my_value = my_chemistry.initialize(a_value)
 if not my_value:
 
 n_points = 200
 fc = FluidContainer(my_chemistry, n_points)
-fc["density"][:] = 1.0
+fc["density"][:] = 1. * mass_h / my_chemistry.density_units
 fc["HI"][:] = 0.76 * fc["density"]
 fc["HII"][:] = tiny_number * fc["density"]
 fc["HM"][:] = tiny_number * fc["density"]
 fc["de"][:] = tiny_number * fc["density"]
 fc["metal"][:] = 0.02041 * fc["density"]
 
-temperature_units = mass_h * ((my_chemistry.length_units/
-                              my_chemistry.time_units)**2) / kboltz
-
 initial_energy = np.logspace(1, 9, n_points) / temperature_units
-
 fc["energy"] = np.copy(initial_energy)
 fc["x-velocity"][:] = 0.0
 fc["y-velocity"][:] = 0.0
 
 fc_last = copy.deepcopy(fc)
 
-dt = 0.1
+dt = 1e4 * yr_to_s / my_chemistry.time_units
 my_time = 0.0
 while True:
     print "t = %.2f Myr." % (my_time * my_chemistry.time_units /
 dbase1   = my_chemistry.density_units*(a_value*my_chemistry.a_units)**3
 cool_unit = (my_chemistry.a_units**5 * xbase1**2 * mass_h**2) / \
   (my_chemistry.time_units**3 * dbase1)
-
-cooling_rate = cool_unit * fc["density"] * fc["energy"] / \
-  fc["cooling_time"]
+  
+cooling_rate = cool_unit * fc["energy"] / \
+  fc["cooling_time"] / fc["density"]
   
 from matplotlib import pyplot
 pyplot.loglog(fc['temperature'], cooling_rate)
 pyplot.xlabel('T [K]')
-pyplot.ylabel('$\\Lambda$ [erg s$^{-1}$ cm$^{-3}$]')
-pyplot.ylim(1e-30, 1e-21)
+pyplot.ylabel('$\\Lambda$ [erg s$^{-1}$ cm$^{3}$]')
 output_file = 'cooling_rate.png'
 print "Writing %s." % output_file
 pyplot.savefig(output_file)