[yt-svn] commit/yt-doc: jzuhone: New cookbook scripts:
Bitbucket
commits-noreply at bitbucket.org
Tue Jul 24 11:13:24 PDT 2012
1 new commit in yt-doc:
https://bitbucket.org/yt_analysis/yt-doc/changeset/b0b01f8d5424/
changeset: b0b01f8d5424
user: jzuhone
date: 2012-07-24 20:03:56
summary: New cookbook scripts:
free_free_field.py: Defines a field for emission of free-free radiation that uses several field parameters.
hse_field.py: A complex field definition that uses the computation of the gravitational potential and pressure gradients to compute the "hydrostatic disequilibrium parameter".
multi_plot_slice_and_proj.py: Takes a single FLASH dataset and combines slices and projections to make a 2x3 multiplot.
rad_velocity.py: Computes the bulk velocity of a sphere, subtracts this velocity from the velocity field of the sphere, and plots a radial profile of the resulting radial velocity.
radial_profile_styles.py: Creates a radial profile of density and then exhibits various linestyles in Matplotlib.
save_profiles.py: Creates radial profiles and saves them to disk, with an ASCII example and a HDF5 example. Shows how to read them back and plot them for good measure.
affected #: 6 files
diff -r 3eaca37968de2b33efd75ed4867172260a29c843 -r b0b01f8d54244542cc60b1e5928d6e93298b2c0d source/cookbook/free_free_field.py
--- /dev/null
+++ b/source/cookbook/free_free_field.py
@@ -0,0 +1,110 @@
+"""
+This is a rather complicated example of how to use field parameters to alter a field.
+The field in question is the monochromatic free-free emission from a hot, optically thin plasma (such as in galaxy clusters) in erg/s/cm^3/keV, or optionally in photons/s/cm^3/keV. (see Rybicki and Lightman 1979). The following field parameters may be specified:
+Z = the mean charge number of the ions in the plasma
+mue = the mean molecular weight for the electrons
+mui = the mean molecular weight for the ions
+Ephoton = the photon energy in keV
+photom_emissivity = logical flag determining whether we use energy or photon emission
+
+In order to get the total luminosity for a data object, we need to add a new derived quantity that multiplies the emission by the cell volumes. We do that here as well.
+"""
+from yt.mods import *
+from yt.utilities.physical_constants import mp # Need to grab the proton mass from the
+ # constants database
+
+# Define the emission field
+
+keVtoerg = 1.602e-9 # Convert energy in keV to energy in erg
+KtokeV = 8.617e-08 # Convert degrees Kelvin to degrees keV
+sqrt3 = na.sqrt(3.)
+expgamma = 1.78107241799 # Exponential of Euler's constant
+
+def _FreeFree_Emission(field, data) :
+
+ if data.has_field_parameter("Z") :
+ Z = data.get_field_parameter("Z")
+ else :
+ Z = 1.077 # Primordial H/He plasma
+
+ if data.has_field_parameter("mue") :
+ mue = data.get_field_parameter("mue")
+ else :
+ mue = 1./0.875 # Primordial H/He plasma
+
+ if data.has_field_parameter("mui") :
+ mui = data.get_field_parameter("mui")
+ else :
+ mui = 1./0.8125 # Primordial H/He plasma
+
+ if data.has_field_parameter("Ephoton") :
+ Ephoton = data.get_field_parameter("Ephoton")
+ else :
+ Ephoton = 1.0 # in keV
+
+ if data.has_field_parameter("photon_emission") :
+ photon_emission = data.get_field_parameter("photon_emission")
+ else :
+ photon_emission = False # Flag for energy or photon emission
+
+ n_e = data["Density"]/(mue*mp)
+ n_i = data["Density"]/(mui*mp)
+ kT = data["Temperature"]*KtokeV
+
+ # Compute the Gaunt factor
+
+ g_ff = na.zeros(kT.shape)
+ g_ff[Ephoton/kT > 1.] = na.sqrt((3./na.pi)*kT[Ephoton/kT > 1.]/Ephoton)
+ g_ff[Ephoton/kT < 1.] = (sqrt3/na.pi)*na.log((4./expgamma) *
+ kT[Ephoton/kT < 1.]/Ephoton)
+
+ eps_E = 1.64e-20*Z*Z*n_e*n_i/na.sqrt(data["Temperature"]) * \
+ na.exp(-Ephoton/kT)*g_ff
+
+ if photon_emission: eps_E /= (Ephoton*keVtoerg)
+
+ return eps_E
+
+add_field("FreeFree_Emission", function=_FreeFree_Emission)
+
+# Define the luminosity derived quantity
+
+def _FreeFreeLuminosity(data) :
+ return (data["FreeFree_Emission"]*data["CellVolume"]).sum()
+
+def _combFreeFreeLuminosity(data, luminosity) :
+ return luminosity.sum()
+
+add_quantity("FreeFree_Luminosity", function=_FreeFreeLuminosity,
+ combine_function=_combFreeFreeLuminosity, n_ret = 1)
+
+pf = load("GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150")
+
+sphere = pf.h.sphere(pf.domain_center, (100., "kpc"))
+
+# Print out the total luminosity at 1 keV for the sphere
+
+print "L_E (1 keV, primordial) = ", sphere.quantities["FreeFree_Luminosity"]()
+
+# The defaults for the field assume a H/He primordial plasma. Let's set the appropriate
+# parameters for a pure hydrogen plasma.
+
+sphere.set_field_parameter("mue", 1.0)
+sphere.set_field_parameter("mui", 1.0)
+sphere.set_field_parameter("Z", 1.0)
+
+print "L_E (1 keV, pure hydrogen) = ", sphere.quantities["FreeFree_Luminosity"]()
+
+# Now let's print the luminosity at an energy of E = 10 keV
+
+sphere.set_field_parameter("Ephoton", 10.0)
+
+print "L_E (10 keV, pure hydrogen) = ", sphere.quantities["FreeFree_Luminosity"]()
+
+# Finally, let's set the flag for photon emission, to get the total number of photons
+# emitted at this energy:
+
+sphere.set_field_parameter("photon_emission", True)
+
+print "L_ph (10 keV, pure hydrogen) = ", sphere.quantities["FreeFree_Luminosity"]()
+
diff -r 3eaca37968de2b33efd75ed4867172260a29c843 -r b0b01f8d54244542cc60b1e5928d6e93298b2c0d source/cookbook/hse_field.py
--- /dev/null
+++ b/source/cookbook/hse_field.py
@@ -0,0 +1,179 @@
+from yt.mods import *
+
+# Define the components of the gravitational acceleration vector field by taking the
+# gradient of the gravitational potential
+
+def _Grav_Accel_x(field, data) :
+
+ # We need to set up stencils
+
+ sl_left = slice(None,-2,None)
+ sl_right = slice(2,None,None)
+ div_fac = 2.0
+
+ dx = div_fac * data['dx'].flat[0]
+
+ gx = data["Grav_Potential"][sl_right,1:-1,1:-1]/dx
+ gx -= data["Grav_Potential"][sl_left, 1:-1,1:-1]/dx
+
+ new_field = na.zeros(data["Grav_Potential"].shape, dtype='float64')
+ new_field[1:-1,1:-1,1:-1] = -gx
+
+ return new_field
+
+def _Grav_Accel_y(field, data) :
+
+ # We need to set up stencils
+
+ sl_left = slice(None,-2,None)
+ sl_right = slice(2,None,None)
+ div_fac = 2.0
+
+ dy = div_fac * data['dy'].flat[0]
+
+ gy = data["Grav_Potential"][1:-1,sl_right,1:-1]/dy
+ gy -= data["Grav_Potential"][1:-1,sl_left ,1:-1]/dy
+
+ new_field = na.zeros(data["Grav_Potential"].shape, dtype='float64')
+ new_field[1:-1,1:-1,1:-1] = -gy
+
+ return new_field
+
+def _Grav_Accel_z(field, data) :
+
+ # We need to set up stencils
+
+ sl_left = slice(None,-2,None)
+ sl_right = slice(2,None,None)
+ div_fac = 2.0
+
+ dz = div_fac * data['dz'].flat[0]
+
+ gz = data["Grav_Potential"][1:-1,1:-1,sl_right]/dz
+ gz -= data["Grav_Potential"][1:-1,1:-1,sl_left ]/dz
+
+ new_field = na.zeros(data["Grav_Potential"].shape, dtype='float64')
+ new_field[1:-1,1:-1,1:-1] = -gz
+
+ return new_field
+
+# Define the components of the pressure gradient field
+
+def _Grad_Pressure_x(field, data) :
+
+ # We need to set up stencils
+
+ sl_left = slice(None,-2,None)
+ sl_right = slice(2,None,None)
+ div_fac = 2.0
+
+ dx = div_fac * data['dx'].flat[0]
+
+ px = data["Pressure"][sl_right,1:-1,1:-1]/dx
+ px -= data["Pressure"][sl_left, 1:-1,1:-1]/dx
+
+ new_field = na.zeros(data["Pressure"].shape, dtype='float64')
+ new_field[1:-1,1:-1,1:-1] = px
+
+ return new_field
+
+def _Grad_Pressure_y(field, data) :
+
+ # We need to set up stencils
+
+ sl_left = slice(None,-2,None)
+ sl_right = slice(2,None,None)
+ div_fac = 2.0
+
+ dy = div_fac * data['dy'].flat[0]
+
+ py = data["Pressure"][1:-1,sl_right,1:-1]/dy
+ py -= data["Pressure"][1:-1,sl_left ,1:-1]/dy
+
+ new_field = na.zeros(data["Pressure"].shape, dtype='float64')
+ new_field[1:-1,1:-1,1:-1] = py
+
+ return new_field
+
+def _Grad_Pressure_z(field, data) :
+
+ # We need to set up stencils
+
+ sl_left = slice(None,-2,None)
+ sl_right = slice(2,None,None)
+ div_fac = 2.0
+
+ dz = div_fac * data['dz'].flat[0]
+
+ pz = data["Pressure"][1:-1,1:-1,sl_right]/dz
+ pz -= data["Pressure"][1:-1,1:-1,sl_left ]/dz
+
+ new_field = na.zeros(data["Pressure"].shape, dtype='float64')
+ new_field[1:-1,1:-1,1:-1] = pz
+
+ return new_field
+
+# Define the "degree of hydrostatic equilibrium" field
+
+def _HSE(field, data) :
+
+ gx = data["Density"]*data["Grav_Accel_x"]
+ gy = data["Density"]*data["Grav_Accel_y"]
+ gz = data["Density"]*data["Grav_Accel_z"]
+
+ hx = data["Grad_Pressure_x"] - gx
+ hy = data["Grad_Pressure_y"] - gy
+ hz = data["Grad_Pressure_z"] - gz
+
+ h = na.sqrt((hx*hx+hy*hy+hz*hz)/(gx*gx+gy*gy+gz*gz))
+
+ return h
+
+# Now add the fields to the database
+
+add_field("Grav_Accel_x", function=_Grav_Accel_x, take_log=False,
+ validators=[ValidateSpatial(1,["Grav_Potential"])])
+
+add_field("Grav_Accel_y", function=_Grav_Accel_y, take_log=False,
+ validators=[ValidateSpatial(1,["Grav_Potential"])])
+
+add_field("Grav_Accel_z", function=_Grav_Accel_z, take_log=False,
+ validators=[ValidateSpatial(1,["Grav_Potential"])])
+
+add_field("Grad_Pressure_x", function=_Grad_Pressure_x, take_log=False,
+ validators=[ValidateSpatial(1,["Pressure"])])
+
+add_field("Grad_Pressure_y", function=_Grad_Pressure_y, take_log=False,
+ validators=[ValidateSpatial(1,["Pressure"])])
+
+add_field("Grad_Pressure_z", function=_Grad_Pressure_z, take_log=False,
+ validators=[ValidateSpatial(1,["Pressure"])])
+
+add_field("HSE", function=_HSE, take_log=False)
+
+# Open two files, one at the beginning and the other at a later time when there's a
+# lot of sloshing going on.
+
+pfi = load("sloshing_low_res_hdf5_plt_cnt_0000")
+pff = load("sloshing_low_res_hdf5_plt_cnt_0350")
+
+# Sphere objects centered at the cluster potential minimum with a radius of 200 kpc
+
+sphere_i = pfi.h.sphere(pfi.domain_center, (200, "kpc"))
+sphere_f = pff.h.sphere(pff.domain_center, (200, "kpc"))
+
+# Average "degree of hydrostatic equilibrium" in these spheres
+
+hse_i = sphere_i.quantities["WeightedAverageQuantity"]("HSE", "CellMass")
+hse_f = sphere_f.quantities["WeightedAverageQuantity"]("HSE", "CellMass")
+
+print "Degree of hydrostatic equilibrium initially: ", hse_i
+print "Degree of hydrostatic equilibrium later: ", hse_f
+
+# Just for good measure, take slices through the center of the domain of the two files
+
+slc_i = SlicePlot(pfi, 2, ["Density","HSE"], center=pfi.domain_center, width=(1.0, "mpc"))
+slc_f = SlicePlot(pff, 2, ["Density","HSE"], center=pff.domain_center, width=(1.0, "mpc"))
+
+slc_i.save("initial")
+slc_f.save("final")
diff -r 3eaca37968de2b33efd75ed4867172260a29c843 -r b0b01f8d54244542cc60b1e5928d6e93298b2c0d source/cookbook/multi_plot_slice_and_proj.py
--- /dev/null
+++ b/source/cookbook/multi_plot_slice_and_proj.py
@@ -0,0 +1,74 @@
+"""
+Title: Basic 2x3 Multi Plot
+Description: This is a simple recipe to show how to open a dataset and then
+ plot a slice through it, centered at its most dense point.
+Outputs:
+"""
+from yt.mods import * # set up our namespace
+import matplotlib.colorbar as cb
+from matplotlib.colors import LogNorm
+
+fn = "GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150" # parameter file to load
+orient = 'horizontal'
+
+pf = load(fn) # load data
+
+# There's a lot in here:
+# From this we get a containing figure, a list-of-lists of axes into which we
+# can place plots, and some axes that we'll put colorbars.
+# We feed it:
+# Number of plots on the x-axis, number of plots on the y-axis, and how we
+# want our colorbars oriented. (This governs where they will go, too.
+# bw is the base-width in inches, but 4 is about right for most cases.
+fig, axes, colorbars = get_multi_plot(3, 2, colorbar=orient, bw = 4)
+
+slc = pf.h.slice(2, 0.0, fields=["Density","Temperature","z-velocity"],
+ center=pf.domain_center)
+proj = pf.h.proj(2, "Density", weight_field="Density", center=pf.domain_center)
+
+slc_frb = slc.to_frb((1.0, "mpc"), 512)
+proj_frb = proj.to_frb((1.0, "mpc"), 512)
+
+dens_axes = [axes[0][0], axes[1][0]]
+temp_axes = [axes[0][1], axes[1][1]]
+vels_axes = [axes[0][2], axes[1][2]]
+
+for dax, tax, vax in zip(dens_axes, temp_axes, vels_axes) :
+
+ dax.xaxis.set_visible(False)
+ dax.yaxis.set_visible(False)
+ tax.xaxis.set_visible(False)
+ tax.yaxis.set_visible(False)
+ vax.xaxis.set_visible(False)
+ vax.yaxis.set_visible(False)
+
+plots = [dens_axes[0].imshow(slc_frb["Density"], origin='lower', norm=LogNorm()),
+ dens_axes[1].imshow(proj_frb["Density"], origin='lower', norm=LogNorm()),
+ temp_axes[0].imshow(slc_frb["Temperature"], origin='lower'),
+ temp_axes[1].imshow(proj_frb["Temperature"], origin='lower'),
+ vels_axes[0].imshow(slc_frb["z-velocity"], origin='lower'),
+ vels_axes[1].imshow(proj_frb["z-velocity"], origin='lower')]
+
+plots[0].set_clim((1.0e-27,1.0e-25))
+plots[0].set_cmap("bds_highcontrast")
+plots[1].set_clim((1.0e-27,1.0e-25))
+plots[1].set_cmap("bds_highcontrast")
+plots[2].set_clim((1.0e7,1.0e8))
+plots[2].set_cmap("hot")
+plots[3].set_clim((1.0e7,1.0e8))
+plots[3].set_cmap("hot")
+plots[4].set_clim((-1.5e7,1.5e7))
+plots[4].set_cmap("gist_rainbow")
+plots[5].set_clim((-1.5e7,1.5e7))
+plots[5].set_cmap("gist_rainbow")
+
+titles=[r'$\mathrm{Density}\ (\mathrm{g\ cm^{-3}})$',
+ r'$\mathrm{Temperature}\ (\mathrm{K})$',
+ r'$\mathrm{z-velocity}\ (\mathrm{cm\ s^{-1}})$']
+
+for p, cax, t in zip(plots[0:6:2], colorbars, titles):
+ cbar = fig.colorbar(p, cax=cax, orientation=orient)
+ cbar.set_label(t)
+
+# And now we're done!
+fig.savefig("%s_3x2" % pf)
diff -r 3eaca37968de2b33efd75ed4867172260a29c843 -r b0b01f8d54244542cc60b1e5928d6e93298b2c0d source/cookbook/rad_velocity.py
--- /dev/null
+++ b/source/cookbook/rad_velocity.py
@@ -0,0 +1,44 @@
+from yt.mods import *
+import matplotlib.pyplot as plt
+
+pf = load("GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150")
+
+# Get the first sphere
+
+sphere0 = pf.h.sphere(pf.domain_center, (500., "kpc"))
+
+# Compute the bulk velocity from the cells in this sphere
+
+bulk_vel = sphere0.quantities["BulkVelocity"]()
+
+# Get the second sphere
+
+sphere1 = pf.h.sphere(pf.domain_center, (500., "kpc"))
+
+# Set the bulk velocity field parameter
+sphere1.set_field_parameter("bulk_velocity", bulk_vel)
+
+# Radial profile without correction
+
+rad_profile0 = BinnedProfile1D(sphere0, 100, "Radiuskpc", 0.0, 500., log_space=False)
+rad_profile0.add_fields("RadialVelocity")
+
+# Radial profile with correction for bulk velocity
+
+rad_profile1 = BinnedProfile1D(sphere1, 100, "Radiuskpc", 0.0, 500., log_space=False)
+rad_profile1.add_fields("RadialVelocity")
+
+# Make a plot using matplotlib
+
+fig = plt.figure()
+ax = fig.add_subplot(111)
+
+# Here we scale the velocities by 1.0e5 to get into km/s
+ax.plot(rad_profile0["Radiuskpc"], rad_profile0["RadialVelocity"]/1.0e5,
+ rad_profile1["Radiuskpc"], rad_profile1["RadialVelocity"]/1.0e5)
+
+ax.set_xlabel(r"$\mathrm{r\ (kpc)}$")
+ax.set_ylabel(r"$\mathrm{v_r\ (km/s)}$")
+ax.legend(["Without Correction", "With Correction"])
+
+fig.savefig("%s_profiles.png" % pf)
\ No newline at end of file
diff -r 3eaca37968de2b33efd75ed4867172260a29c843 -r b0b01f8d54244542cc60b1e5928d6e93298b2c0d source/cookbook/radial_profile_styles.py
--- /dev/null
+++ b/source/cookbook/radial_profile_styles.py
@@ -0,0 +1,60 @@
+from yt.mods import *
+import matplotlib.pyplot as plt
+
+pf = load("GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150")
+
+# Get a sphere object
+
+sphere = pf.h.sphere(pf.domain_center, (500., "kpc"))
+
+# Bin up the data from the sphere into a radial profile
+
+rad_profile = BinnedProfile1D(sphere, 100, "Radiuskpc", 0.0, 500., log_space=False)
+rad_profile.add_fields("Density","Temperature")
+
+# Make plots using matplotlib
+
+fig = plt.figure()
+ax = fig.add_subplot(111)
+
+# Plot the density as a log-log plot using the default settings
+dens_plot = ax.loglog(rad_profile["Radiuskpc"], rad_profile["Density"])
+
+# Here we set the labels of the plot axes
+
+ax.set_xlabel(r"$\mathrm{r\ (kpc)}$")
+ax.set_ylabel(r"$\mathrm{\rho\ (g\ cm^{-3})}$")
+
+# Save the default plot
+
+fig.savefig("density_profile_default.png" % pf)
+
+# The "dens_plot" object is a list of plot objects. In our case we only have one,
+# so we index the list by '0' to get it.
+
+# Plot using dashed red lines
+
+dens_plot[0].set_linestyle("--")
+dens_plot[0].set_color("red")
+
+fig.savefig("density_profile_dashed_red.png")
+
+# Increase the line width and add points in the shape of x's
+
+dens_plot[0].set_linewidth(5)
+dens_plot[0].set_marker("x")
+dens_plot[0].set_markersize(10)
+
+fig.savefig("density_profile_thick_with_xs.png")
+
+# Now get rid of the line on the axes plot
+
+ax.lines = []
+
+# Since the rad_profile object also includes the standard deviation in each bin,
+# we'll use these as errorbars. We have to make a new plot for this:
+
+dens_err_plot = ax.errorbar(rad_profile["Radiuskpc"], rad_profile["Density"],
+ yerr=rad_profile["Density_std"])
+
+fig.savefig("density_profile_with_errorbars.png")
\ No newline at end of file
diff -r 3eaca37968de2b33efd75ed4867172260a29c843 -r b0b01f8d54244542cc60b1e5928d6e93298b2c0d source/cookbook/save_profiles.py
--- /dev/null
+++ b/source/cookbook/save_profiles.py
@@ -0,0 +1,69 @@
+from yt.mods import *
+import matplotlib.pyplot as plt
+import h5py
+
+pf = load("GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150")
+
+# Get a sphere
+
+sp = pf.h.sphere(pf.domain_center, (500., "kpc"))
+
+# Radial profile from the sphere
+
+rad_profile = BinnedProfile1D(sp, 100, "Radiuskpc", 0.0, 500., log_space=False)
+
+# Adding density and temperature fields to the profile
+
+rad_profile.add_fields(["Density","Temperature"])
+
+# Write profiles to ASCII file
+
+rad_profile.write_out("%s_profile.dat" % pf, bin_style="center")
+
+# Write profiles to HDF5 file
+
+rad_profile.write_out_h5("%s_profile.h5" % pf, bin_style="center")
+
+# Now we will show how using NumPy, h5py, and Matplotlib the data in these
+# files may be plotted.
+
+# Plot density from ASCII file
+
+# Open the text file using NumPy's "loadtxt" method. In order to get the
+# separate columns into separate NumPy arrays, it is essential to set unpack=True.
+
+r, dens, std_dens, temp, std_temp = \
+ na.loadtxt("sloshing_nomag2_hdf5_plt_cnt_0150_profile.dat", unpack=True)
+
+fig1 = plt.figure()
+
+ax = fig1.add_subplot(111)
+ax.plot(r, dens)
+ax.set_xlabel(r"$\mathrm{r\ (kpc)}$")
+ax.set_ylabel(r"$\mathrm{\rho\ (g\ cm^{-3})}$")
+ax.set_title("Density vs. Radius")
+fig1.savefig("%s_dens.png" % pf)
+
+# Plot temperature from HDF5 file
+
+# Get the file handle
+
+f = h5py.File("%s_profile.h5" % pf, "r")
+
+# Get the radius and temperature arrays from the file handle
+
+r = f["/Radiuskpc-1d"].attrs["x-axis-Radiuskpc"][:]
+temp = f["/Radiuskpc-1d/Temperature"][:]
+
+# Close the file handle
+
+f.close()
+
+fig2 = plt.figure()
+
+ax = fig2.add_subplot(111)
+ax.plot(r, temp)
+ax.set_xlabel(r"$\mathrm{r\ (kpc)}$")
+ax.set_ylabel(r"$\mathrm{T\ (K)}$")
+ax.set_title("Temperature vs. Radius")
+fig2.savefig("%s_temp.png" % pf)
Repository URL: https://bitbucket.org/yt_analysis/yt-doc/
--
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