now with correct a for Hernquist profile

plots/densityprofile.py

@@ -13,7 +13,17 @@ M = np.sum(m)
 
 print("Number of particles:",len(m))
 
-rmax = 32
+halfmassradius = 0
+tmpmass = 0
+rsorted = np.sort(r)
+for i in range(len(rsorted)):
+    tmpmass += m[i]
+    if tmpmass > M/2:
+        halfmassradius = rsorted[i]
+        break
+print("Half mass radius:",halfmassradius)
+
+rmax = 128
 dr = 0.075
 nbins = int(rmax/dr)
 
@@ -23,7 +33,7 @@ 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)
+a = halfmassradius/(1+np.sqrt(2))
 rhoHernquist = M / (2*np.pi) * a / (bins[1:]-dr/2) / (bins[1:]-dr/2 + a)**3
 
 print("Number of bins:",nbins)
@@ -33,16 +43,6 @@ print("Bin size:",dr)
 #     if hist[i] == 0:
 #         hist[i] = np.nan
 
-halfmassradius = 0
-tmpmass = 0
-rsorted = np.sort(r)
-for i in range(len(rsorted)):
-    tmpmass += m[i]
-    if tmpmass > M/2:
-        halfmassradius = rsorted[i]
-        break
-print("Half mass radius:",halfmassradius)
-
 # 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