Improved the plotting routine to add a measurement of the \sigma of the error and a histogram.

examples/plot.py

@@ -3,6 +3,7 @@
 import matplotlib
 matplotlib.use("Agg")
 from pylab import *
+import numpy
 # tex stuff                                                                                                                                
 #rc('font',**{'family':'serif','serif':['Palatino']})                                                                                      
 params = {'axes.labelsize': 16,
@@ -33,42 +34,150 @@ rc('font', family='serif')
 import sys
 import os
 
+# Read Quickshed accelerations
 data=loadtxt("particle_dump.dat")
 id = data[:,0]
-accx_e=data[:,1]
-accy_e=data[:,2]
-accz_e=data[:,3]
+accx_e=data[:,4]
+accy_e=data[:,5]
+accz_e=data[:,6]
 
-accx_bh=data[:,4]
-accy_bh=data[:,5]
-accz_bh=data[:,6]
+accx_bh=data[:,7]
+accy_bh=data[:,8]
+accz_bh=data[:,9]
 
-accx_new=data[:,7]
-accy_new=data[:,8]
-accz_new=data[:,9]
+accx_new=data[:,10]
+accy_new=data[:,11]
+accz_new=data[:,12]
 
+# Sort accelerations
+rank = argsort(id)
+id = id[rank]
+accx_e = accx_e[rank]
+accy_e = accy_e[rank]
+accz_e = accz_e[rank]
 
+accx_bh = accx_bh[rank]
+accy_bh = accy_bh[rank]
+accz_bh = accz_bh[rank]
 
+accx_new = accx_new[rank]
+accy_new = accy_new[rank]
+accz_new = accz_new[rank]
+
+
+# Read Gadget accelerations
+data=loadtxt("particle_dump_gadget.dat")
+id = data[:,0]
+accx_g=data[:,4]
+accy_g=data[:,5]
+accz_g=data[:,6]
+
+# Sort accelerations
+rank = argsort(id)
+id = id[rank]
+accx_g = accx_g[rank]
+accy_g = accy_g[rank]
+accz_g = accz_g[rank]
+
+# Build error ------------------------------------------------
+
+errx_bh = (accx_bh - accx_e )/abs(accx_e) 
+erry_bh = (accy_bh - accy_e )/abs(accy_e) 
+errz_bh = (accz_bh - accz_e )/abs(accz_e) 
+
+errx_new = (accx_new - accx_e )/abs(accx_e) 
+erry_new = (accy_new - accy_e )/abs(accy_e) 
+errz_new = (accz_new - accz_e )/abs(accz_e) 
+
+errx_g = (accx_g - accx_e )/abs(accx_e) 
+erry_g = (accy_g - accy_e )/abs(accy_e) 
+errz_g = (accz_g - accz_e )/abs(accz_e) 
+
+
+# Statistics
+meanx_bh = mean(errx_bh)
+stdx_bh = sqrt(mean(errx_bh**2) - meanx_bh**2)
+meany_bh = mean(erry_bh)
+stdy_bh = std(erry_bh)
+meanz_bh = mean(errz_bh)
+stdz_bh = std(errz_bh)
+
+meanx_new = mean(errx_new)
+stdx_new = std(errx_new)
+meany_new = mean(erry_new)
+stdy_new = std(erry_new)
+meanz_new = mean(errz_new)
+stdz_new = std(errz_new)
+
+meanx_g = mean(errx_g)
+stdx_g = std(errx_g)
+meany_g = mean(erry_g)
+stdy_g = std(erry_g)
+meanz_g = mean(errz_g)
+stdz_g = std(errz_g)
+
+# Plot -------------------------------------------------------
 figure(frameon=True)
 
-subplot(311)
-plot(id, (accx_bh - accx_e )/abs(accx_e) + 1, 'rx')
-plot(id, (accx_new - accx_e )/abs(accx_e)  + 1 , 'b.')
+subplot(311, title="Acceleration along X")
+plot(id, errx_g , 'gs')
+plot(id, errx_bh , 'rx')
+plot(id, errx_new , 'b.')
+
+text(id[-1], 0.18, "Gadget: $%5.3f\\pm%5.3f$ \n B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanx_g, stdx_g, meanx_bh, stdx_bh, meanx_new, stdx_new), backgroundcolor="w", va="top", ha="right" )
 
-ylim(0.8, 1.2)
+
+ylim(-0.2, 0.2)
+xlim(0,id[-1])
 grid()
 
-subplot(312)
-plot(id, (accy_bh - accy_e )/abs(accy_e ) + 1, 'rx')
-plot(id, (accy_new - accy_e )/abs(accy_e ) + 1, 'b.')
-ylim(0.8, 1.2)
+subplot(312, title="Acceleration along Y")
+plot(id, erry_g , 'gs')
+plot(id, erry_bh , 'rx')
+plot(id, erry_new , 'b.')
+text(id[-1], 0.18, "Gadget: $%5.3f\\pm%5.3f$ \n B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meany_g, stdy_g, meany_bh, stdy_bh, meany_new, stdy_new), backgroundcolor="w", va="top", ha="right" )
+
+ylim(-0.2, 0.2)
+xlim(0,id[-1])
+
 grid()
 
-subplot(313)
-plot(id, (accz_bh - accz_e )/abs(accz_e ) + 1, 'rx')
-plot(id, (accz_new - accz_e )/abs(accz_e)  + 1, 'b.')
-ylim(0.8, 1.2)
+subplot(313, title="Acceleration along Z")
+plot(id, errz_g , 'gs')
+plot(id, errz_bh , 'rx')
+plot(id, errz_new , 'b.')
+text(id[-1], 0.18, "Gadget: $%5.3f\\pm%5.3f$ \n B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meany_g, stdy_g, meany_bh, stdy_bh, meany_new, stdy_new), backgroundcolor="w", va="top", ha="right" )
+
+ylim(-0.2, 0.2)
+xlim(0,id[-1])
 grid()
 
 savefig("accelerations.png")
 
+
+
+
+# Plot -------------------------------------------------------
+bins = linspace(-3, 3, 10000)
+
+
+figure(frameon=True)
+subplot(311, title="Acceleration along X")
+hist(errx_g, bins=bins, normed=1, histtype='step', rwidth=0.01, color='g')
+hist(errx_bh, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r')
+hist(errx_new, bins=bins, normed=1, histtype='step', rwidth=0.01, color='b')
+xlim(-0.03, 0.03)
+
+subplot(312, title="Acceleration along Y")
+hist(erry_g, bins=bins, normed=1, histtype='step', rwidth=0.01, color='g')
+hist(erry_bh, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r')
+hist(erry_new, bins=bins, normed=1, histtype='step', rwidth=0.01, color='b')
+xlim(-0.03, 0.03)
+
+subplot(313, title="Acceleration along Z")
+hist(errz_g, bins=bins, normed=1, histtype='step', rwidth=0.01, color='g')
+hist(errz_bh, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r')
+hist(errz_new, bins=bins, normed=1, histtype='step', rwidth=0.01, color='b')
+xlim(-0.03, 0.03)
+
+savefig("histogram.png")