compute mass flow form surface density (radial profile)

1 files changed, 57 insertions, 0 deletions

diff --git a/comp_MdotR.py b/comp_MdotR.py
new file mode 100644
index 0000000..86d8f4b
--- /dev/null
+++ b/comp_MdotR.py
@@ -0,0 +1,57 @@
+"""
+Compute Mass Flow Profile -- Mdot(R).
+"""
+
+import numpy as np
+import sys
+from helpers import outputs
+from math import pi
+
+
+"""
+Main Routine. Parse Arguments. Loop Desired Outputs.
+"""
+def main():
+
+    # Set Defaults
+    iout = 1
+
+    # Parse Arguments
+    if len(sys.argv) == 2:
+        iout = sys.argv[1]
+
+    # 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 )
+        do_comp_MdotR(iiout)
+
+
+"""
+Load Data. Compute MdotR.
+"""
+def do_comp_MdotR(iout):
+
+    npz1 = np.load("SigmaR_%05d.npz" % (iout - 1))
+    npz2 = np.load("SigmaR_%05d.npz" % iout)
+
+    dSigmaR = npz2["SigmaR"] - npz1["SigmaR"]
+    dt = npz2["tout"] - npz1["tout"]
+    dr = npz2["rbins"][1] - npz2["rbins"][0]
+    MdotR = pi * dr * dSigmaR / dt
+
+    np.savez('MdotR_%05d.npz' % iout, \
+        MdotR=MdotR,\
+        rbins=npz2["rbins"],\
+        tout1=npz1["tout"],\
+        tout2=npz2["tout"])
+
+
+"""
+Jump into Main().
+"""
+if __name__ == "__main__":
+    main()