convert cuts to cgs; add temperature

1 files changed, 55 insertions, 8 deletions

diff --git a/mkcuts_dvp.py b/mkcuts_dvp.py
index 508c5b8..5d1cd42 100644
--- a/mkcuts_dvp.py
+++ b/mkcuts_dvp.py
@@ -3,17 +3,19 @@ Make X and Z Cuts:
 - Density
 - Velocity
 - Pressure
+- Temperature (Sound Speed)
 """
 
 import matplotlib as mpl
 from pymses import RamsesOutput
 from pymses.filters import CellsToPoints
-from pymses.analysis.visualization import Camera, ScalarOperator, SliceMap
+from pymses.utils import constants as C
+from pymses.analysis.visualization import Camera, ScalarOperator, SliceMap, \
+                                          FractionOperator
 import matplotlib.pyplot as plt
 import numpy as np
 import sys
 
-
 """
 Main Routine. Define Cameras. Take Slices. Plot.
 """
@@ -50,6 +52,7 @@ def main():
     op_rho = ScalarOperator(func_rho)
     op_vel = ScalarOperator(func_vel)
     op_pre = ScalarOperator(func_pre)
+    op_cs2 = FractionOperator(func_pre, func_rho)
 
     # Compute Slice Maps
     rhoZ = SliceMap(source, camZ, op_rho, z=0.0)
@@ -58,6 +61,33 @@ def main():
     velX = SliceMap(source, camX, op_vel, z=0.0)
     preZ = SliceMap(source, camZ, op_pre, z=0.0)
     preX = SliceMap(source, camX, op_pre, z=0.0)
+    cs2Z = SliceMap(source, camZ, op_cs2, z=0.0)
+    cs2X = SliceMap(source, camX, op_cs2, z=0.0)
+
+    # Convert Human Mind Parsable Units Pt 1 - (CGS)
+    rhoZ *= output.info["unit_density"].express(C.g_cc)
+    rhoX *= output.info["unit_density"].express(C.g_cc)
+    velZ *= output.info["unit_velocity"].express(C.cm / C.s)
+    velX *= output.info["unit_velocity"].express(C.cm / C.s)
+    preZ *= output.info["unit_pressure"].express(C.barye)
+    preX *= output.info["unit_pressure"].express(C.barye)
+
+    # Convert Human Mind Parsable Units Pt 2 - (km/s)
+    # Sound Speed
+    factor = output.info["unit_pressure"].express(C.barye) / \
+             output.info["unit_density"].express(C.g_cc)
+    cs2Z *= factor / 100.**2. / 1000.**2.
+    cs2X *= factor / 100.**2. / 1000.**2.
+    # Total Flow Speed
+    velZ *= 1.0 / 100. / 1000.
+    velX *= 1.0 / 100. / 1000.
+
+    # Compute Temperature
+    gamma = 1.4
+    TX = (cs2X * 100.**2. * 1000.**2.) * 2. * C.mH.express(C.kg) / \
+         gamma / C.kB.express(C.cm**2 * C.kg / C.s**2 / C.K)
+    TZ = (cs2Z * 100.**2. * 1000.**2.) * 2. * C.mH.express(C.kg) / \
+         gamma / C.kB.express(C.cm**2 * C.kg / C.s**2 / C.K)
 
     # Log10?
     rhoZ = np.log10(rhoZ)
@@ -95,7 +125,7 @@ def main():
     plt.imshow(rhoZ, cmap=rho_cm, extent=ext, interpolation='none')
     plt.colorbar()
     plt.grid()
-    plt.title('Log10 Density Cut, t=%.2f' % output.info["time"])
+    plt.title('Log10 Density Cut [g/cc], t=%.2f' % output.info["time"])
     plt.xlabel('X')
     plt.ylabel('Y')
 
@@ -103,7 +133,7 @@ def main():
     plt.imshow(rhoX, cmap=rho_cm, extent=ext, interpolation='none')
     plt.colorbar()
     plt.grid()
-    plt.title('Log10 Density Cut, t=%.2f' % output.info["time"])
+    plt.title('Log10 Density Cut [g/cc], t=%.2f' % output.info["time"])
     plt.xlabel('Y')
     plt.ylabel('Z')
 
@@ -112,7 +142,7 @@ def main():
     plt.imshow(velZ, cmap=vel_cm, extent=ext)
     plt.colorbar()
     plt.grid()
-    plt.title('Angular Velocity Cut, t=%.2f' % output.info["time"])
+    plt.title('Total Flow Speed Cut [km/s], t=%.2f' % output.info["time"])
     plt.xlabel('X')
     plt.ylabel('Y')
 
@@ -120,7 +150,7 @@ def main():
     plt.imshow(velX, cmap=vel_cm, extent=ext)
     plt.colorbar()
     plt.grid()
-    plt.title('Angular XY-Velocity Cut, t=%.2f' % output.info["time"])
+    plt.title('Total Flow Speed Cut [km/s], t=%.2f' % output.info["time"])
     plt.xlabel('Y')
     plt.ylabel('Z')
 
@@ -129,7 +159,7 @@ def main():
     plt.imshow(preZ, cmap=pre_cm, extent=ext, interpolation='none')
     plt.colorbar()
     plt.grid()
-    plt.title('Log10 Pressure Cut, t=%.2f' % output.info["time"])
+    plt.title('Log10 Pressure Cut [debye], t=%.2f' % output.info["time"])
     plt.xlabel('X')
     plt.ylabel('Y')
 
@@ -137,7 +167,24 @@ def main():
     plt.imshow(preX, cmap=pre_cm, extent=ext, interpolation='none')
     plt.colorbar()
     plt.grid()
-    plt.title('Log10 Pressure Cut, t=%.2f' % output.info["time"])
+    plt.title('Log10 Pressure Cut [debye], t=%.2f' % output.info["time"])
+    plt.xlabel('Y')
+    plt.ylabel('Z')
+
+    # Plot Temperature
+    f7 = plt.figure(7)
+    plt.imshow(TZ, cmap=pre_cm, extent=ext, interpolation='none')
+    plt.colorbar()
+    plt.grid()
+    plt.title('Temperature Cut [K], t=%.2f' % output.info["time"])
+    plt.xlabel('X')
+    plt.ylabel('Y')
+
+    f8 = plt.figure(8)
+    plt.imshow(TX, cmap=pre_cm, extent=ext, interpolation='none')
+    plt.colorbar()
+    plt.grid()
+    plt.title('Temperature Cut [K], t=%.2f' % output.info["time"])
     plt.xlabel('Y')
     plt.ylabel('Z')