diff options
| -rw-r--r-- | DiskXYZ.py | 89 | ||||
| -rw-r--r-- | plot_quad_mag.py | 45 |
2 files changed, 134 insertions, 0 deletions
@@ -209,6 +209,15 @@ class DiskIo(DiskBase): self.angular_momentum_lost_to_halo = npz["angular_momentum_lost_to_halo"] self.infall_rate = npz["infall_rate"] self.info = npz["info"][()] + self.h_xy = npz["h_xy"] + self.T_mid_xy = npz["T_mid_xy"] + self.Omega_xy = npz["Omega_xy"] + self.rho_mid_xy = npz["rho_mid_xy"] + self.n_mid_xy = npz["n_mid_xy"] + self.xe_mid_xy = npz["xe_mid_xy"] + self.etaB_mid_xy = npz["etaB_mid_xy"] + self.ua2_mid_xy = npz["ua2_mid_xy"] + self.ReM_mid_xy = npz["ReM_mid_xy"] def load_stats(self): """Load Stats.""" @@ -510,6 +519,86 @@ class DiskReduce(DiskReduceBase): class DiskPlots(DiskReduce): """Plotting Routines.""" + def plot_ReM_mid_xy(self, ax_in=None, clim=None): + ext = [np.min(self.xy0["x"]), np.max(self.xy0["x"]), \ + np.min(self.xy0["y"]), np.max(self.xy0["y"])] + if not ax_in: + fig = plt.figure() + ax = fig.add_subplot(1,1,1) + else: + ax = ax_in + im = ax.imshow(np.rot90(np.log10(self.ReM_mid_xy.reshape(self.nx, self.ny))), \ + extent=ext, interpolation='none') + ax.set_xlabel('X [AU]') + ax.set_ylabel('Y [AU]') + ax.set_title('Log10 Magnetic Reynolds Number [-]') + ax.grid() + if clim: + im.set_clim(clim) + plt.colorbar(im, ax=ax) + if not ax_in: + plt.show() + + def plot_xe_mid_xy(self, ax_in=None, clim=None): + ext = [np.min(self.xy0["x"]), np.max(self.xy0["x"]), \ + np.min(self.xy0["y"]), np.max(self.xy0["y"])] + if not ax_in: + fig = plt.figure() + ax = fig.add_subplot(1,1,1) + else: + ax = ax_in + im = ax.imshow(np.rot90(np.log10(self.xe_mid_xy.reshape(self.nx, self.ny))), \ + extent=ext, interpolation='none') + ax.set_xlabel('X [AU]') + ax.set_ylabel('Y [AU]') + ax.set_title('Log10 Midplane Electron Fraction [-]') + ax.grid() + if clim: + im.set_clim(clim) + plt.colorbar(im, ax=ax) + if not ax_in: + plt.show() + + def plot_n_mid_xy(self, ax_in=None, clim=None): + ext = [np.min(self.xy0["x"]), np.max(self.xy0["x"]), \ + np.min(self.xy0["y"]), np.max(self.xy0["y"])] + if not ax_in: + fig = plt.figure() + ax = fig.add_subplot(1,1,1) + else: + ax = ax_in + im = ax.imshow(np.rot90(np.log10(self.n_mid_xy.reshape(self.nx, self.ny))), \ + extent=ext, interpolation='none') + ax.set_xlabel('X [AU]') + ax.set_ylabel('Y [AU]') + ax.set_title('Log10 Midplane H2 Number Density [1/cm^3]') + ax.grid() + if clim: + im.set_clim(clim) + plt.colorbar(im, ax=ax) + if not ax_in: + plt.show() + + def plot_T_mid_xy(self, ax_in=None, clim=None): + ext = [np.min(self.xy0["x"]), np.max(self.xy0["x"]), \ + np.min(self.xy0["y"]), np.max(self.xy0["y"])] + if not ax_in: + fig = plt.figure() + ax = fig.add_subplot(1,1,1) + else: + ax = ax_in + im = ax.imshow(np.rot90(self.T_mid_xy.reshape(self.nx, self.ny)), \ + extent=ext, interpolation='none') + ax.set_xlabel('X [AU]') + ax.set_ylabel('Y [AU]') + ax.set_title('Midplane Temperature [K]') + ax.grid() + if clim: + im.set_clim(clim) + plt.colorbar(im, ax=ax) + if not ax_in: + plt.show() + def plot_rho_xy(self, ax_in=None, clim=None): ext = [np.min(self.xy0["x"]), np.max(self.xy0["x"]), \ np.min(self.xy0["y"]), np.max(self.xy0["y"])] diff --git a/plot_quad_mag.py b/plot_quad_mag.py new file mode 100644 index 0000000..1067592 --- /dev/null +++ b/plot_quad_mag.py @@ -0,0 +1,45 @@ +import matplotlib as mpl; mpl.use('agg') +import matplotlib.pyplot as plt +import DiskXYZ +import DiskRTZ +import numpy as np +import argparse +import sys +from math import pi +from copy import deepcopy + +# Parse Arguments +parser = argparse.ArgumentParser() +parser.add_argument("imin", type=int, help='First Output') +parser.add_argument("imax", type=int, help='Last Output') +args = parser.parse_args() + +# Sanity Checks +if args.imin > args.imax: + print("Cannot Work With Imin > Imax. Use -h for help.") + sys.exit(-1) + +# Build Output Range +iouts = range(args.imin, args.imax+1) + +# Make Plots +for iout in iouts: + print "Rendering Figure %i/%i." % (iout, iouts[-1]) + disk = DiskXYZ.Disk(iout) + disk.load_npz() + + fig = plt.figure(figsize=(16.0, 12.0)) + ax1 = fig.add_subplot(2,2,1) + ax2 = fig.add_subplot(2,2,2) + ax3 = fig.add_subplot(2,2,3) + ax4 = fig.add_subplot(2,2,4) + disk.plot_T_mid_xy(ax1) + disk.plot_n_mid_xy(ax2) + disk.plot_xe_mid_xy(ax3) + + disk.plot_ReM_mid_xy(ax4) + plt.suptitle("t=%0.2f [code] / t=%0.2f [yr] / %05d [iout] / M_infall=%1.2e [Mstar/yr]" % \ + (disk.info["time"], disk.info["time"]/2./pi, iout, disk.infall_rate)) + fig.savefig('Quad_Mag_%05d.png' % iout) + + del disk |