summaryrefslogtreecommitdiffstats
path: root/plot_quad_classic.py
diff options
context:
space:
mode:
Diffstat (limited to 'plot_quad_classic.py')
-rw-r--r--plot_quad_classic.py164
1 files changed, 164 insertions, 0 deletions
diff --git a/plot_quad_classic.py b/plot_quad_classic.py
new file mode 100644
index 0000000..a541242
--- /dev/null
+++ b/plot_quad_classic.py
@@ -0,0 +1,164 @@
+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)
+
+# Determine Density Limits
+first = True
+for iout in iouts:
+ print "Scanning MinMax %i/%i." % (iout, iouts[-1])
+ disk = DiskXYZ.Disk(iout)
+ disk2 = DiskRTZ.Disk(iout)
+ disk.load_npz_minmax()
+ disk.load_npz_stats()
+ disk2.load_npz_minmax()
+ disk2.load_npz_stats()
+
+ if first:
+ rho_r_lo = disk.rho_r_no0_min; rho_r_hi = disk.rho_r_no0_max
+ rho_xy_lo = disk.rho_xy_min; rho_xy_hi = disk.rho_xy_max
+ ylo_pos = disk2.mdotr_r_g0_min; yhi_pos = disk2.mdotr_r_g0_max
+ ylo_neg = disk2.mdotr_r_l0_min; yhi_neg = disk2.mdotr_r_l0_max
+ # note, could be zero. but this is the ONLY chance it can be zero
+ infall_rate_r_max = disk.infall_rate * \
+ disk2.mass_r_max / disk2.total_mass
+ infall_rate_r_min = disk.infall_rate * \
+ disk2.mass_r_min / disk2.total_mass
+ first = False
+
+ else:
+ if disk.rho_r_no0_min < rho_r_lo: rho_r_lo = disk.rho_r_no0_min
+ if disk.rho_r_no0_max > rho_r_hi: rho_r_hi = disk.rho_r_no0_max
+ if disk.rho_xy_min < rho_xy_lo: rho_xy_lo = disk.rho_xy_min
+ if disk.rho_xy_max > rho_xy_hi: rho_xy_hi = disk.rho_xy_max
+
+ infall_rate_r_max_loc = disk.infall_rate * \
+ disk2.mass_r_max / disk2.total_mass
+ infall_rate_r_min_loc = disk.infall_rate * \
+ disk2.mass_r_min / disk2.total_mass
+
+ # we don't want to record =0 infall rates
+ if infall_rate_r_max_loc > 0:
+ if infall_rate_r_max_loc > infall_rate_r_max:
+ infall_rate_r_max = infall_rate_r_max_loc
+ if infall_rate_r_min_loc > 0:
+ if infall_rate_r_min_loc < infall_rate_r_min:
+ infall_rate_r_min = infall_rate_r_min_loc
+
+ if np.isnan(ylo_pos) and not np.isnan(disk2.mdotr_r_g0_min):
+ ylo_pos = disk2.mdotr_r_g0_min
+ if np.isnan(yhi_pos) and not np.isnan(disk2.mdotr_r_g0_max):
+ yhi_pos = disk2.mdotr_r_g0_max
+ if np.isnan(ylo_neg) and not np.isnan(disk2.mdotr_r_l0_min):
+ ylo_neg = disk2.mdotr_r_l0_min
+ if np.isnan(yhi_neg) and not np.isnan(disk2.mdotr_r_l0_max):
+ yhi_neg = disk2.mdotr_r_l0_max
+
+ if disk2.mdotr_r_g0_min < ylo_pos: ylo_pos = disk2.mdotr_r_g0_min
+ if disk2.mdotr_r_g0_max > yhi_pos: yhi_pos = disk2.mdotr_r_g0_max
+ if disk2.mdotr_r_l0_min < ylo_neg: ylo_neg = disk2.mdotr_r_l0_min
+ if disk2.mdotr_r_l0_max > yhi_neg: yhi_neg = disk2.mdotr_r_l0_max
+
+# rho limits
+rho_r_lim = [0.5*rho_r_lo, 2.0*rho_r_hi]
+rho_xy_lim = [np.log10(rho_xy_lo), np.log10(rho_xy_hi)]
+
+# mdotr limits (copied from plot_mdotr_r.py)
+lim_pos = np.array([ylo_pos, yhi_pos])
+lim_neg = np.array([ylo_neg, yhi_neg])
+lim_neg = -lim_neg[::-1]
+lim_tot = np.array([np.minimum(lim_neg[0], lim_pos[0]), \
+ np.maximum(lim_neg[1], lim_pos[1])])
+
+# adjust limits in we display case of >0 infall_rates
+# can be zero ONLY if the initial value was zero
+# in this case, disregard it
+if infall_rate_r_min > 0:
+ lim_tot[0] = np.minimum(lim_tot[0], infall_rate_r_min)
+if infall_rate_r_max > 0:
+ lim_tot[1] = np.maximum(lim_tot[1], infall_rate_r_max)
+
+# Make Plots
+for iout in iouts:
+ print "Rendering Figure %i/%i." % (iout, iouts[-1])
+ disk = DiskXYZ.Disk(iout)
+ disk2 = DiskRTZ.Disk(iout)
+ disk.load_npz()
+ disk2.load_npz()
+
+ # Compute Infall Rate per Radial Bin
+ infall_rate_r = disk.infall_rate * disk2.mass_r / disk2.total_mass
+
+ # Soften
+ if disk2.mdotr_r.shape[0] > 7:
+ kernel = np.array([1,1,1,1,1], dtype='float64') / 5.
+ disk2.mdotr_r = np.convolve(disk2.mdotr_r, kernel, mode='same')
+
+ # split out mdotr_r>0 and <0
+ # copied from plot_dotr_r.py
+ # @todo should probably push this into some sort of wrapper function...
+ posbool = disk2.mdotr_r>0; negbool = disk2.mdotr_r<0
+ disk_pos = deepcopy(disk2); disk_neg = deepcopy(disk2)
+ disk_abs = deepcopy(disk2)
+ disk_pos.r["r"] = disk_pos.r["r"][posbool]
+ disk_neg.r["r"] = disk_neg.r["r"][negbool]
+ disk_pos.mdotr_r = disk_pos.mdotr_r[posbool]
+ disk_neg.mdotr_r = -disk_neg.mdotr_r[negbool]
+ disk_abs.mdotr_r = np.abs(disk_abs.mdotr_r)
+
+ # resume figuring
+ 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_rho_r(ax1, ylim=rho_r_lim)
+ disk.plot_rho_xy(ax2, clim=rho_xy_lim)
+
+ # mdotr fluff
+ # copied from plot_dotr_r.py
+ ax3.hold(True)
+ h0 = disk_abs.plot_mdotr_r(ax_in=ax3, ylog=True)
+ h1 = disk_pos.plot_mdotr_r(ax_in=ax3, ylog=True)
+ h2 = disk_neg.plot_mdotr_r(ax_in=ax3, ylog=True)
+ h3 = ax3.plot(disk2.r["r"], infall_rate_r, 'mo-')
+ ax3.hold(False)
+ ax3.set_title('Radial Mass Flux')
+ ax3.grid(True)
+ h1.set_label('Outward'); h2.set_label('Inward')
+ h0.set_linestyle('-'); h0.set_color('k'), h0.set_marker('')
+ h1.set_linestyle(''); h2.set_linestyle('')
+ h1.set_marker('s'); h2.set_marker('s')
+ h1.set_color('b'); h2.set_color('r')
+ ax3.set_xlim([np.min(disk2.r["r"]), np.max(disk2.r["r"])])
+ ax3.legend(loc='best')
+ ax3.set_ylim(lim_tot)
+
+ # resume figuring
+ disk.plot_Q_r(ax4)
+ plt.suptitle("t=%0.2f [code] / t=%0.2f [yr] / iout=%05d / nstep=%06d / M_infall=%1.2e [Mstar/yr]" % \
+ (disk.info["time"], disk.info["time"]/2./pi, iout, disk.info["nstep_coarse"], disk.infall_rate))
+ plt.savefig('Quad_Classic_%05d.png' % iout)
+ plt.clf()
+ plt.close()
+ del disk
+ del disk2