From a587704ed50ef08728534925c9a39bf7d35575ca Mon Sep 17 00:00:00 2001
From: Noe Brucy <noe.brucy@cea.fr>
Date: Tue, 14 May 2019 16:57:16 +0200
Subject: [PATCH] Use center values for interpolation

---
 disk_postprocess.py | 41 +++++++++++++++++++----------------------
 1 file changed, 19 insertions(+), 22 deletions(-)

diff --git a/disk_postprocess.py b/disk_postprocess.py
index cacce52..191362f 100644
--- a/disk_postprocess.py
+++ b/disk_postprocess.py
@@ -675,7 +675,7 @@ def disk_prop(
         "time": time,
         "mass_disk": total_mass_disk,
         "mass_box": total_mass,
-        "rad": rad[0 : nb_bin - 1],
+        "rad": rad[:-1],
         "center": pos_star,
         "alpha_rey": alpha_rey_rad,
         "alpha_rey_mean": alpha_rey_mean,
@@ -873,8 +873,7 @@ def disk_pdf(
     nb_bin_hist=50,
     tag="",
     rad_min=0.08,
-    rad_max=0.35,
-    interpol=True,
+    rad_max=0.4,
 ):
 
     # Load property file
@@ -905,8 +904,18 @@ def disk_pdf(
 
     coldens_map = maps_disk["coldens_z"]
     im_extent = maps_disk["im_extent"]
+    # radius of the corner of the box plus a margin
+    rad_box = (
+        np.sqrt((im_extent[1] - pos_star[0]) ** 2 + (im_extent[3] - pos_star[1]) ** 2)
+        + 0.1
+    )
+    # radial bins
     rad_bins = prop_disk["rad"]
+    rad_bins = 0.5 * (rad_bins[0:-1] + rad_bins[1:])
+    rad_bins = np.concatenate(([0.0], rad_bins, [rad_box]))
     coldens_rad = prop_disk["coldens"]
+    # Add value for borders
+    coldens_rad = np.concatenate(([coldens_rad[0]], coldens_rad))
 
     x = np.linspace(im_extent[0], im_extent[1], coldens_map.shape[0])
     y = np.linspace(im_extent[2], im_extent[3], coldens_map.shape[1])
@@ -920,25 +929,13 @@ def disk_pdf(
     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]
+    # Retrieve the mean of the bins
+    # use linear interpolation to improve accuracy
+    coldens_mean = (rad_bins[bins_flat + 1] - rr_flat) * coldens_rad[bins_flat]
+    coldens_mean = (
+        coldens_mean + (rr_flat - rad_bins[bins_flat]) * coldens_rad[bins_flat + 1]
+    )
+    coldens_mean = coldens_mean / (rad_bins[bins_flat + 1] - rad_bins[bins_flat])
 
     coldens_mean_map = np.reshape(coldens_mean, coldens_map.shape)