make and save cuts for multiple outputs

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: