prettier plot, rename script files

plots/densityprofile.py

@@ -6,51 +6,42 @@ data = np.loadtxt('../data/data.txt')
 x = data[:,2]
 y = data[:,3]
 z = data[:,4]
+r = np.sqrt(x**2 + y**2 + z**2)
 
 m = data[:,1]
 M = np.sum(m)
 
 print("Number of particles:",len(m))
 
-maxr = 16
-dr = maxr/100
+rmax = 32
+dr = 0.075
+nbins = int(rmax/dr)
 
-a = dr
+hist, bins = np.histogram(r, bins=nbins, range=(0,rmax), weights=m)
 
-r = np.arange(dr, maxr, dr)
-rhoHernquist = M / (2*np.pi) * a / r / (r + a)**4
-rho = np.zeros(len(r))
+for i in range(len(hist)):
+    volume = 4/3*np.pi*(bins[i+1]**3 - bins[i]**3)
+    hist[i] /= volume
 
-for i in range(len(x)):
-    ri2 = x[i]**2 + y[i]**2 + z[i]**2
-    for j in range(len(r)):
-        if (ri2 > r[j]**2 and ri2 < (r[j] + dr)**2):
-            rho[j] += m[i]
-            break
+a = np.diff(bins)
+rhoHernquist = M / (2*np.pi) * a / (bins[1:]-dr/2) / (bins[1:]-dr/2 + a)**3
 
-
-for i in range(len(r)-1):
-    volume = 4*np.pi*(r[i+1]**3 - r[i]**3)/3
-    rho[i] /= volume
-
-#rho /= 4*np.pi*r**2*dr
-
-# Normalize
-# rhoHernquist /= np.sum(rhoHernquist)
-# rho /= np.sum(rho)
-
-print("Number of bins:",len(rho))
+print("Number of bins:",nbins)
 print("Bin size:",dr)
 
+for i in range(len(hist)):
+    if hist[i] == 0:
+        hist[i] = np.nan
+
 # plot poissonian error bars
-plt.errorbar(r, rhoHernquist, yerr=np.sqrt(rhoHernquist), color='grey', fmt='.', label='Hernquist with Poisson error', markersize=0)
+plt.errorbar(bins[1:], rhoHernquist, yerr=np.sqrt(rhoHernquist), color='grey', fmt='.', label='Hernquist with Poisson error', markersize=0)
 # plot Hernquist profile
-plt.loglog(r, rhoHernquist, 'r.-', label='Hernquist Density Profile')
+plt.loglog(bins[1:]-dr/2, rhoHernquist, 'r-', label='Hernquist Density Profile', linewidth=1)
 # plot density profile of data 
-plt.loglog(r, rho, 'b.-', label='Density Profile from data')
+plt.loglog(bins[1:]-dr/2, hist, 'b-', label='Density Profile from data', linewidth=1)
 
-plt.xlabel('Radius r')
-plt.ylabel('Density rho(r)')
+plt.xlabel(r'Radius $r$')
+plt.ylabel(r'Density $\rho(r)$')
 
 plt.legend()
 plt.savefig('densityprofile.png', dpi=200)

plots/densityprofile2.py → plots/olddensityprofile.py

@@ -6,39 +6,48 @@ data = np.loadtxt('../data/data.txt')
 x = data[:,2]
 y = data[:,3]
 z = data[:,4]
-r = np.sqrt(x**2 + y**2 + z**2)
 
 m = data[:,1]
 M = np.sum(m)
 
 print("Number of particles:",len(m))
 
-rmax = 32
-dr = 0.075
-nbins = int(rmax/dr)
+maxr = 16
+dr = maxr/100
 
-hist, bins = np.histogram(r, bins=nbins, range=(0,rmax), weights=m)
+a = dr
 
-for i in range(len(hist)):
-    volume = 4/3*np.pi*(bins[i+1]**3 - bins[i]**3)
-    hist[i] /= volume
+r = np.arange(dr, maxr, dr)
+rhoHernquist = M / (2*np.pi) * a / r / (r + a)**4
+rho = np.zeros(len(r))
 
-a = np.diff(bins)
-rhoHernquist = M / (2*np.pi) * a / (bins[1:]-dr/2) / (bins[1:]-dr/2 + a)**3
+for i in range(len(x)):
+    ri2 = x[i]**2 + y[i]**2 + z[i]**2
+    for j in range(len(r)):
+        if (ri2 > r[j]**2 and ri2 < (r[j] + dr)**2):
+            rho[j] += m[i]
+            break
 
-print("Number of bins:",nbins)
-print("Bin size:",dr)
 
-for i in range(len(hist)):
-    if hist[i] == 0:
-        hist[i] = np.nan
+for i in range(len(r)-1):
+    volume = 4*np.pi*(r[i+1]**3 - r[i]**3)/3
+    rho[i] /= volume
+
+#rho /= 4*np.pi*r**2*dr
+
+# Normalize
+# rhoHernquist /= np.sum(rhoHernquist)
+# rho /= np.sum(rho)
+
+print("Number of bins:",len(rho))
+print("Bin size:",dr)
 
 # plot poissonian error bars
-plt.errorbar(bins[1:], rhoHernquist, yerr=np.sqrt(rhoHernquist), color='grey', fmt='.', label='Hernquist with Poisson error', markersize=0)
+plt.errorbar(r, rhoHernquist, yerr=np.sqrt(rhoHernquist), color='grey', fmt='.', label='Hernquist with Poisson error', markersize=0)
 # plot Hernquist profile
-plt.loglog(bins[1:]-dr/2, rhoHernquist, 'r-', label='Hernquist Density Profile', linewidth=1)
+plt.loglog(r, rhoHernquist, 'r.-', label='Hernquist Density Profile')
 # plot density profile of data 
-plt.loglog(bins[1:]-dr/2, hist, 'b-', label='Density Profile from data', linewidth=1)
+plt.loglog(r, rho, 'b.-', label='Density Profile from data')
 
 plt.xlabel('Radius r')
 plt.ylabel('Density rho(r)')