diff options
Diffstat (limited to 'plot_quad_classic.py')
| -rw-r--r-- | plot_quad_classic.py | 164 |
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 |