diff --git a/disk_postprocess.py b/disk_postprocess.py index 62323b0..cacce52 100644 --- a/disk_postprocess.py +++ b/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