aboutsummaryrefslogtreecommitdiffstats
path: root/mkcuts_dvpT.py
diff options
context:
space:
mode:
Diffstat (limited to 'mkcuts_dvpT.py')
-rw-r--r--mkcuts_dvpT.py326
1 files changed, 168 insertions, 158 deletions
diff --git a/mkcuts_dvpT.py b/mkcuts_dvpT.py
index 6bfeaad..b1f1d02 100644
--- a/mkcuts_dvpT.py
+++ b/mkcuts_dvpT.py
@@ -24,7 +24,8 @@ def main():
# Parse Arguments
parser = argparse.ArgumentParser()
- parser.add_argument("iout", type=int, help='Output', default=1)
+ parser.add_argument("ilo", type=int, help='First Output', default=1)
+ parser.add_argument("ihi", type=int, help='Last Output', default=1)
parser.add_argument("-d", action="store_true", help='Show Density Cuts')
parser.add_argument("-v", action="store_true", help='Show Velocity Cuts')
parser.add_argument("-p", action="store_true", help='Show Pressure Cuts')
@@ -40,20 +41,18 @@ def main():
if not args.save and not args.show:
print "At least one of --save or --show is required."
sys.exit(0)
+ if args.ilo != args.ihi and args.show:
+ print "Can only show single output."
+ sys.exit(0)
- # Give Feedback
- print "Creating Cuts for Output %i." % args.iout
- print ""
-
- # Link AMR Data
- output = RamsesOutput(".", args.iout)
- source = output.amr_source(["rho", "vel", "P"])
+ # Build Output List
+ iouts = range(args.ilo, args.ihi + 1)
# Define Cameras
camZ = Camera(center=[0.5, 0.5, 0.5], line_of_sight_axis='z',
- region_size=[1., 1.], up_vector='y', map_max_size=args.px, log_sensitive=True)
+ region_size=[1., 1.], up_vector='y', map_max_size=args.px, log_sensitive=True)
camX = Camera(center=[0.5, 0.5, 0.5], line_of_sight_axis='x',
- region_size=[1., 1.], up_vector='y', map_max_size=args.px, log_sensitive=True)
+ region_size=[1., 1.], up_vector='y', map_max_size=args.px, log_sensitive=True)
# Functions to Get Data
func_rho = lambda dset: dset["rho"]
@@ -67,158 +66,169 @@ def main():
op_pre = ScalarOperator(func_pre)
op_cs2 = FractionOperator(func_pre, func_rho)
- # Compute Slice Maps
- rhoZ = SliceMap(source, camZ, op_rho, z=0.0)
- rhoX = SliceMap(source, camX, op_rho, z=0.0)
- velZ = SliceMap(source, camZ, op_vel, z=0.0)
- 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)
- rhoX = np.log10(rhoX)
- preZ = np.log10(preZ)
- preX = np.log10(preX)
-
- # Mark Arrays vs. NaN
- rhoZ = np.ma.masked_array(rhoZ, mask=np.isnan(rhoZ))
- rhoX = np.ma.masked_array(rhoX, mask=np.isnan(rhoX))
- velZ = np.ma.masked_array(velZ, mask=np.isnan(velZ))
- velX = np.ma.masked_array(velX, mask=np.isnan(velX))
- preZ = np.ma.masked_array(preZ, mask=np.isnan(preZ))
- preX = np.ma.masked_array(preX, mask=np.isnan(preX))
-
- # Set Masked Colormap
- # cmap = mpl.cm.hot
- cmap = mpl.cm.jet
- cmap.set_bad([0.5, 0.5, 0.5])
-
- # Colormaps:
- # http://matplotlib.org/examples/pylab_examples/show_colormaps.html
- # rho_cm = 'hot'
- # vel_cm = 'hot'
- # pre_cm = 'hot'
- rho_cm = cmap
- vel_cm = cmap
- pre_cm = cmap
-
- # Compute Image Extent
- ext = output.info["boxlen"] * np.array([-0.5, 0.5, -0.5, 0.5])
-
- # Plot Density
- if args.d:
- f1 = plt.figure(1)
- plt.imshow(rhoZ, cmap=rho_cm, extent=ext, interpolation='none')
- plt.colorbar()
- plt.grid()
- plt.title('Log10 Density Cut [g/cc], t=%.2f' % output.info["time"])
- plt.xlabel('X [AU]')
- plt.ylabel('Y [AU]')
-
- f2 = plt.figure(2)
- plt.imshow(rhoX, cmap=rho_cm, extent=ext, interpolation='none')
- plt.colorbar()
- plt.grid()
- plt.title('Log10 Density Cut [g/cc], t=%.2f' % output.info["time"])
- plt.xlabel('Y [AU]')
- plt.ylabel('Z [AU]')
-
- # Plot Velocity
- if args.v:
- f3 = plt.figure(3)
- plt.imshow(velZ, cmap=vel_cm, extent=ext)
- plt.colorbar()
- plt.grid()
- plt.title('Total Flow Speed Cut [km/s], t=%.2f' % output.info["time"])
- plt.xlabel('X [AU]')
- plt.ylabel('Y [AU]')
-
- f4 = plt.figure(4)
- plt.imshow(velX, cmap=vel_cm, extent=ext)
- plt.colorbar()
- plt.grid()
- plt.title('Total Flow Speed Cut [km/s], t=%.2f' % output.info["time"])
- plt.xlabel('Y [AU]')
- plt.ylabel('Z [AU]')
-
- # Plot Pressure
- if args.p:
- f5 = plt.figure(5)
- plt.imshow(preZ, cmap=pre_cm, extent=ext, interpolation='none')
- plt.colorbar()
- plt.grid()
- plt.title('Log10 Pressure Cut [debye], t=%.2f' % output.info["time"])
- plt.xlabel('X [AU]')
- plt.ylabel('Y [AU]')
-
- f6 = plt.figure(6)
- plt.imshow(preX, cmap=pre_cm, extent=ext, interpolation='none')
- plt.colorbar()
- plt.grid()
- plt.title('Log10 Pressure Cut [debye], t=%.2f' % output.info["time"])
- plt.xlabel('Y [AU]')
- plt.ylabel('Z [AU]')
-
- # Plot Temperature
- if args.T:
- 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 [AU]')
- plt.ylabel('Y [AU]')
-
- 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 [AU]')
- plt.ylabel('Z [AU]')
-
- # Save Figures
- if args.save:
+ # The Master Loop
+ for iout in iouts:
+
+ # Give Feedback
+ print "Creating Cuts for Output %i." % iout
+ print ""
+
+ # Link AMR Data
+ output = RamsesOutput(".", iout)
+ source = output.amr_source(["rho", "vel", "P"])
+
+ # Compute Slice Maps
+ rhoZ = SliceMap(source, camZ, op_rho, z=0.0)
+ rhoX = SliceMap(source, camX, op_rho, z=0.0)
+ velZ = SliceMap(source, camZ, op_vel, z=0.0)
+ 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)
+ rhoX = np.log10(rhoX)
+ preZ = np.log10(preZ)
+ preX = np.log10(preX)
+
+ # Mark Arrays vs. NaN
+ rhoZ = np.ma.masked_array(rhoZ, mask=np.isnan(rhoZ))
+ rhoX = np.ma.masked_array(rhoX, mask=np.isnan(rhoX))
+ velZ = np.ma.masked_array(velZ, mask=np.isnan(velZ))
+ velX = np.ma.masked_array(velX, mask=np.isnan(velX))
+ preZ = np.ma.masked_array(preZ, mask=np.isnan(preZ))
+ preX = np.ma.masked_array(preX, mask=np.isnan(preX))
+
+ # Set Masked Colormap
+ # cmap = mpl.cm.hot
+ cmap = mpl.cm.jet
+ cmap.set_bad([0.5, 0.5, 0.5])
+
+ # Colormaps:
+ # http://matplotlib.org/examples/pylab_examples/show_colormaps.html
+ # rho_cm = 'hot'
+ # vel_cm = 'hot'
+ # pre_cm = 'hot'
+ rho_cm = cmap
+ vel_cm = cmap
+ pre_cm = cmap
+
+ # Compute Image Extent
+ ext = output.info["boxlen"] * np.array([-0.5, 0.5, -0.5, 0.5])
+
+ # Plot Density
if args.d:
- f1.savefig('cut_dxy_%05d.png' % args.iout)
- f2.savefig('cut_dzy_%05d.png' % args.iout)
+ f1 = plt.figure(1)
+ plt.imshow(rhoZ, cmap=rho_cm, extent=ext, interpolation='none')
+ plt.colorbar()
+ plt.grid()
+ plt.title('Log10 Density Cut [g/cc], t=%.2f' % output.info["time"])
+ plt.xlabel('X [AU]')
+ plt.ylabel('Y [AU]')
+
+ f2 = plt.figure(2)
+ plt.imshow(rhoX, cmap=rho_cm, extent=ext, interpolation='none')
+ plt.colorbar()
+ plt.grid()
+ plt.title('Log10 Density Cut [g/cc], t=%.2f' % output.info["time"])
+ plt.xlabel('Y [AU]')
+ plt.ylabel('Z [AU]')
+
+ # Plot Velocity
if args.v:
- f3.savefig('cut_vxy_%05d.png' % args.iout)
- f4.savefig('cut_vzy_%05d.png' % args.iout)
+ f3 = plt.figure(3)
+ plt.imshow(velZ, cmap=vel_cm, extent=ext)
+ plt.colorbar()
+ plt.grid()
+ plt.title('Total Flow Speed Cut [km/s], t=%.2f' % output.info["time"])
+ plt.xlabel('X [AU]')
+ plt.ylabel('Y [AU]')
+
+ f4 = plt.figure(4)
+ plt.imshow(velX, cmap=vel_cm, extent=ext)
+ plt.colorbar()
+ plt.grid()
+ plt.title('Total Flow Speed Cut [km/s], t=%.2f' % output.info["time"])
+ plt.xlabel('Y [AU]')
+ plt.ylabel('Z [AU]')
+
+ # Plot Pressure
if args.p:
- f5.savefig('cut_pxy_%05d.png' % args.iout)
- f6.savefig('cut_pzy_%05d.png' % args.iout)
+ f5 = plt.figure(5)
+ plt.imshow(preZ, cmap=pre_cm, extent=ext, interpolation='none')
+ plt.colorbar()
+ plt.grid()
+ plt.title('Log10 Pressure Cut [debye], t=%.2f' % output.info["time"])
+ plt.xlabel('X [AU]')
+ plt.ylabel('Y [AU]')
+
+ f6 = plt.figure(6)
+ plt.imshow(preX, cmap=pre_cm, extent=ext, interpolation='none')
+ plt.colorbar()
+ plt.grid()
+ plt.title('Log10 Pressure Cut [debye], t=%.2f' % output.info["time"])
+ plt.xlabel('Y [AU]')
+ plt.ylabel('Z [AU]')
+
+ # Plot Temperature
if args.T:
- f7.savefig('cut_Txy_%05d.png' % args.iout)
- f8.savefig('cut_Tzy_%05d.png' % args.iout)
+ 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 [AU]')
+ plt.ylabel('Y [AU]')
+
+ 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 [AU]')
+ plt.ylabel('Z [AU]')
+
+ # Save Figures
+ if args.save:
+ if args.d:
+ f1.savefig('cut_dxy_%05d.png' % iout)
+ f2.savefig('cut_dzy_%05d.png' % iout)
+ if args.v:
+ f3.savefig('cut_vxy_%05d.png' % iout)
+ f4.savefig('cut_vzy_%05d.png' % iout)
+ if args.p:
+ f5.savefig('cut_pxy_%05d.png' % iout)
+ f6.savefig('cut_pzy_%05d.png' % iout)
+ if args.T:
+ f7.savefig('cut_Txy_%05d.png' % iout)
+ f8.savefig('cut_Tzy_%05d.png' % iout)
# Show Figures?
if args.show: