Improve coldens fluctu histo (pdf)

disk_postprocess.py

@@ -872,7 +872,9 @@ def disk_pdf(
     interactive=False,
     nb_bin_hist=50,
     tag="",
-    rad_min=0.1,
+    rad_min=0.08,
+    rad_max=0.35,
+    interpol=True,
 ):
 
     # Load property file
@@ -915,12 +917,34 @@ def disk_pdf(
     bins = np.zeros(rr.shape, dtype=int)
     for r in rad_bins[1:]:
         bins = bins + (rr >= r).astype(int)
+    bins_flat = bins.flatten()
+    rr_flat = rr.flatten()
+
+    # Retrieve the mean of the bi
+    if interpol:  # use linear interpolation to improve accuracy
+        rad_bins_ext = np.append(
+            rad_bins,
+            np.sqrt(
+                (im_extent[1] - pos_star[0]) ** 2 + (im_extent[3] - pos_star[1]) ** 2
+            ),
+        )
+        coldens_rad_ext = np.append(coldens_rad, coldens_rad[-1])
+        coldens_mean = (rad_bins_ext[bins_flat + 1] - rr_flat) * coldens_rad[bins_flat]
+        coldens_mean = (
+            coldens_mean
+            + (rr_flat - rad_bins[bins_flat]) * coldens_rad_ext[bins_flat + 1]
+        )
+        coldens_mean = coldens_mean / (
+            rad_bins_ext[bins_flat + 1] - rad_bins[bins_flat]
+        )
+    else:
+        coldens_mean = coldens_rad[bins_flat]
 
-    coldens_mean = coldens_rad[bins.flatten()]
     coldens_mean_map = np.reshape(coldens_mean, coldens_map.shape)
 
-    #  Kill fluctuations for r < rad_min
+    #  Kill fluctuations for r < rad_min and r > rad_max
     coldens_map[rr < rad_min] = coldens_mean_map[rr < rad_min]
+    coldens_map[rr > rad_max] = coldens_mean_map[rr > rad_max]
 
     # Histogramm : density fluctuation distribution function
     dcoldens = np.log10(coldens_map.flatten() / coldens_mean)
@@ -929,7 +953,9 @@ def disk_pdf(
     P.xlabel(r"$\log(N / \bar{N})$")
     P.ylabel("Probability density")
     P.title(title)
-    P.hist(dcoldens, nb_bin_hist, range=(-1, 2))
+    P.hist(
+        dcoldens[(rr_flat > rad_min) & (rr_flat < rad_max)], nb_bin_hist, range=(-1, 3)
+    )
 
     if interactive:
         pass