aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorVolker Hoffmann <volker@cheleb.net>2013-04-02 13:19:27 +0200
committerVolker Hoffmann <volker@cheleb.net>2013-04-02 13:19:27 +0200
commitc6f57b21e95a1eb14116fa01da73a173ad50bf59 (patch)
tree68f8c70f9b4f9a692bb0ada92fc94c005fe7e85c
parentcf4c74a2a73e41bd5d044257a58ffa9c6e9465cb (diff)
compute and plot MdotR properly (Mdot=int_A rho*v*dA, i.e. flux through cross section)
-rw-r--r--comp_MdotR.py23
-rw-r--r--plot_MdotR.py7
-rw-r--r--plots2low.py15
3 files changed, 17 insertions, 28 deletions
diff --git a/comp_MdotR.py b/comp_MdotR.py
index cbeabfb..63ce460 100644
--- a/comp_MdotR.py
+++ b/comp_MdotR.py
@@ -23,10 +23,9 @@ def main():
# Compute Requested Outputs
iouts = outputs(iout)
- # Loop Over Outputs, Skip First
- for iiout in iouts[1:]:
- print "Computing Mass Flow Profile for Outputs %i -> %i." % \
- ( (iiout - 1), iiout )
+ # Loop Over Outputs
+ for iiout in iouts:
+ print "Computing Mass Flow Profile for Output %i." % iiout
do_comp_MdotR(iiout)
@@ -35,21 +34,15 @@ Load Data. Compute MdotR.
"""
def do_comp_MdotR(iout):
- npz1 = np.load("SigmaR_%05d.npz" % (iout - 1))
- npz2 = np.load("SigmaR_%05d.npz" % iout)
+ npzSigma = np.load("SigmaR_%05d.npz" % iout)
+ npzVeloc = np.load("VelocitiesR_%05d.npz" % iout)
- dSigmaR = npz2["SigmaR"] - npz1["SigmaR"]
- dt = npz2["tout"] - npz1["tout"]
- dr = npz2["rbins"][1] - npz2["rbins"][0]
- # MdotR = pi * (rhi**2. - rlo**2.) * dSigmaR / dt
- # rhi**2. - rlo**2. = 2 * r * dr
- MdotR = pi * 2. * npz1["rbins"] * dr * dSigmaR / dt
+ MdotR = npzSigma["SigmaR"] * npzVeloc["vr_r"] * 2. * pi * npzVeloc["rbins"]
np.savez('MdotR_%05d.npz' % iout, \
MdotR=MdotR,\
- rbins=npz2["rbins"],\
- tout1=npz1["tout"],\
- tout2=npz2["tout"])
+ rbins=npzVeloc["rbins"],\
+ tout=npzVeloc["tout"])
"""
diff --git a/plot_MdotR.py b/plot_MdotR.py
index cd32716..df830ef 100644
--- a/plot_MdotR.py
+++ b/plot_MdotR.py
@@ -29,9 +29,6 @@ pgroup.add_argument("--save", action="store_true", \
args = parser.parse_args()
# Sanity Checks
-if args.imin < 2:
- print("Require Imin>=2.")
- sys.exit(-1)
if args.imin > args.imax:
print("Cannot Work With Imin > Imax. Use -h for help.")
sys.exit(-1)
@@ -92,8 +89,8 @@ else:
vmin = np.array([])
vmax = np.array([])
for iout in iouts:
- print "Writing Mdot(R) for Output %i->%i." % (iout-1, iout)
- vmin, vmax = plots2.MdotR(iout0=iout, iout1=2,\
+ print "Writing Mdot(R) for Output %i." % iout
+ vmin, vmax = plots2.MdotR(iout0=iout, iout1=1,\
ymin=ymin, ymax=ymax,\
vmin=vmin, vmax=vmax,\
save=True, show=False)
diff --git a/plots2low.py b/plots2low.py
index 4d1788e..77c7bef 100644
--- a/plots2low.py
+++ b/plots2low.py
@@ -93,26 +93,25 @@ def MdotR(iout0, iout1, ymin=np.nan, ymax=np.nan,\
ax = fig.add_subplot(nrow, ncol, isub)
# Plot to Axis
- ax.plot(npz0["rbins"], np.abs(npz0["MdotR"])*2.*pi, 'bs-',\
+ ax.plot(npz0["rbins"], np.abs(npz0["MdotR"]), 'bs-',\
label='Current')
- ax.plot(npz1["rbins"], np.abs(npz1["MdotR"]*2.*pi), 'ks-', linewidth=0.5,\
+ ax.plot(npz1["rbins"], np.abs(npz1["MdotR"]), 'ks-', linewidth=0.5,\
label='Reference')
if showminmax:
- ax.plot(npz0["rbins"], vmin*2.*pi, 'gs-', linewidth=0.5, label='Minimum')
- ax.plot(npz0["rbins"], vmax*2.*pi, 'rs-', linewidth=0.5, label='Maximum')
+ ax.plot(npz0["rbins"], vmin, 'gs-', linewidth=0.5, label='Minimum')
+ ax.plot(npz0["rbins"], vmax, 'rs-', linewidth=0.5, label='Maximum')
# Axis Properties
ax.set_yscale('log')
ax.set_xlim([0,4])
ax.set_xlabel('Radius [AU]')
ax.set_ylabel('Mass Flow [Mstar/yr]')
- ax.set_title('Mass Flow / dt=%.2f-%.2f [code] / dt=%.2f-%.2f [yr]' % \
- (npz0["tout2"], npz0["tout1"],\
- npz0["tout2"]/2./pi, npz0["tout1"]/2./pi))
+ ax.set_title('Mass Flow / t=%.2f [code] / t=%.2f [yr]' % \
+ (npz0["tout"], npz0["tout"]/2./pi))
ax.legend(loc='lower center')
ax.grid(True)
if not (np.isnan(ymin) and np.isnan(ymax)):
- ax.set_ylim([ymin*2.*pi,ymax*2.*pi])
+ ax.set_ylim([ymin,ymax])
# Close Files
npz0.close()