implemented calculation of half mass radius

plots/densityprofile.py

@@ -29,9 +29,19 @@ 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
+# for i in range(len(hist)):
+#     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)