diff --git a/plots/densityprofile.png b/plots/densityprofile.png
index 3f3076d912434d96e1ec9421f0cb092266851530..ce7f7065bba88fbdb463434a86f3899c6d2f7701 100644
Binary files a/plots/densityprofile.png and b/plots/densityprofile.png differ
diff --git a/plots/densityprofile.py b/plots/densityprofile.py
index b134c7d9cd885ff93fcfe184b6e621d65997fd89..3c6560f80cc5230c3038a9f7cdeafe896dd58984 100644
--- a/plots/densityprofile.py
+++ b/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)
diff --git a/plots/densityprofile2.py b/plots/densityprofile2.py
deleted file mode 100644
index ce7851d75678079b7a0df00945bfc5740d0d4975..0000000000000000000000000000000000000000
--- a/plots/densityprofile2.py
+++ /dev/null
@@ -1,48 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-
-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)
-
-hist, bins = np.histogram(r, bins=nbins, range=(0,rmax), weights=m)
-
-for i in range(len(hist)):
-    volume = 4/3*np.pi*(bins[i+1]**3 - bins[i]**3)
-    hist[i] /= volume
-
-a = np.diff(bins)
-rhoHernquist = M / (2*np.pi) * a / (bins[1:]-dr/2) / (bins[1:]-dr/2 + a)**3
-
-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(bins[1:], 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)
-# plot density profile of 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.legend()
-plt.savefig('densityprofile.png', dpi=200)
-plt.show()
\ No newline at end of file
diff --git a/plots/olddensityprofile.py b/plots/olddensityprofile.py
new file mode 100644
index 0000000000000000000000000000000000000000..b134c7d9cd885ff93fcfe184b6e621d65997fd89
--- /dev/null
+++ b/plots/olddensityprofile.py
@@ -0,0 +1,57 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+data = np.loadtxt('../data/data.txt')
+
+x = data[:,2]
+y = data[:,3]
+z = data[:,4]
+
+m = data[:,1]
+M = np.sum(m)
+
+print("Number of particles:",len(m))
+
+maxr = 16
+dr = maxr/100
+
+a = dr
+
+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(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
+
+
+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(r, 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')
+# plot density profile of data 
+plt.loglog(r, rho, 'b.-', label='Density Profile from data')
+
+plt.xlabel('Radius r')
+plt.ylabel('Density rho(r)')
+
+plt.legend()
+plt.savefig('densityprofile.png', dpi=200)
+plt.show()
\ No newline at end of file