plot dM/dt for masses in plot_stats

1 files changed, 32 insertions, 4 deletions

diff --git a/plot_stats.py b/plot_stats.py
index 1802d8d..45eddcd 100644
--- a/plot_stats.py
+++ b/plot_stats.py
@@ -62,33 +62,61 @@ for iout in iouts:
     mass_total[ii] = disk.mass_halo
     ii += 1
 
+# Calculate Time Derivatives
+# @todo - What about negative derivatives?
+# @todo - Get them in the plot somehow; color code \pm?
+delta_mass_total_lost_to_star = np.diff(mass_total_lost_to_star) / \
+                                np.diff(time)
+delta_mass_total_lost_to_halo = np.diff(mass_total_lost_to_halo) / \
+                                np.diff(time)
+delta_mass_disk = np.diff(mass_disk) / np.diff(time)
+delta_mass_halo = np.diff(mass_halo) / np.diff(time)
+delta_mass_total = np.diff(mass_total) / np.diff(time)
+
 # Make Plots
-fig = plt.figure()
-ax = fig.add_subplot(1,1,1)
+fig = plt.figure(figsize=(16.0, 6.0))
+ax = fig.add_subplot(1,2,1)
 ax2 = ax.twiny()
+ax3 = fig.add_subplot(1,2,2)
+ax4 = ax3.twiny()
 if args.disk:
     ax.plot(time/2./pi, mass_disk, 'bs-', label='Disk Mass')
+    ax3.plot(time[1::]/2./pi, delta_mass_disk, 'bs-', label='d(Disk Mass)/dt')
 if args.halo:
     ax.plot(time/2./pi, mass_halo, 'rs-', label='Halo Mass')
+    ax3.plot(time[1::]/2./pi, delta_mass_halo, 'rs-', label='d(Halo Mass)/dt')
 if args.box:
     ax.plot(time/2./pi, mass_total, 'gs-', label='Box Mass')
+    ax3.plot(time[1::]/2./pi, delta_mass_total, 'gs-', label='d(Box Mass)/dt')
 if args.lost_halo:
     ax.plot(time/2./pi, mass_total_lost_to_halo, 'ms-', \
         label='Total Mass Lost to Halo')
+    ax3.plot(time[1::]/2./pi, delta_mass_total_lost_to_halo, 'ms-', \
+             label='d(Total Mass Lost to Halo)/dt')
 if args.lost_star:
     ax.plot(time/2./pi, mass_total_lost_to_star, 'ks-', \
         label='Total Mass Lost to Star')
+    ax3.plot(time[1::]/2./pi, 2.*pi*delta_mass_total_lost_to_star, 'ks-', \
+             label='d(Total Mass Lost to Star)/dt')
 fig.suptitle('Mass Budget')
 ax2.plot(time, np.ones_like(time))
 ax2.cla()
 ax2.set_xlabel('Time [code]')
-ax.grid()
-ax.set_yscale('log')
+ax4.plot(time, np.ones_like(time))
+ax4.cla()
+ax4.set_xlabel('Time [code]')
+ax.grid(); ax.set_yscale('log')
+ax3.grid(); ax3.set_yscale('log')
 ax2.set_ylim(ax.get_ylim())
 ax2.set_yscale(ax.get_yscale())
+ax4.set_ylim(ax3.get_ylim())
+ax4.set_yscale(ax3.get_yscale())
 ax.set_xlabel('Time [yr]')
 ax.set_ylabel('Mass [Mstar]')
 ax.legend(loc='best')
+ax3.set_xlabel('Time [yr]')
+ax3.set_ylabel('d(Mass)/dt [Mstar/yr]')
+ax3.legend(loc='best')
 if args.save:
     plt.savefig('Stats.pdf')
 if args.show: