summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorVolker Hoffmann <volker@cheleb.net>2013-08-16 16:08:55 +0200
committerVolker Hoffmann <volker@cheleb.net>2013-08-16 16:08:55 +0200
commitbf6c9cba14b1ca7b3b39c6898ae43e3c7c5e3d51 (patch)
treedf367b51a3c032d6c8c2eb8faca984a94f3be5ed
parent6852254705b5d035ef53f731bad9e3b0e1f967e0 (diff)
load and plot magneic Reynolds number (and others) in post
-rw-r--r--DiskXYZ.py89
-rw-r--r--plot_quad_mag.py45
2 files changed, 134 insertions, 0 deletions
diff --git a/DiskXYZ.py b/DiskXYZ.py
index f9f6134..5894978 100644
--- a/DiskXYZ.py
+++ b/DiskXYZ.py
@@ -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