Investigations with inversion method (for high concentrations) have been performed with APAVER 2021 data from UBSediFlow system.

Model/GrainSizeTools.py

@@ -155,15 +155,33 @@ class demodul_granulo:
         self.cumul_interpolated = f_logsize(self.sizes_log_resampled)
         
         # Computing volume and number probability from interpolated data
-        self.resampled_proba_vol = np.diff(self.cumul_interpolated)
+        self.resampled_proba_vol = np.abs(np.diff(self.cumul_interpolated))   # np.abs is added because probability p
+                                                                              # must be positive in np.random.choice
+                                                                              # used in function below "demodulate"
         # Computing number probability,
         ss = np.sum(self.resampled_proba_vol /(np.exp(self.sizes_log_resampled[:-1]))**3)
         self.resampled_proba_num = (self.resampled_proba_vol / (np.exp(self.sizes_log_resampled[:-1]))**3) / ss
-        
+
+        # fig, ax = plt.subplots(1, 2)
+        # ax[0].plot(list(range(len(self.cumul_interpolated))), self.cumul_interpolated)
+        # ax[0].set_xlabel("n° point")
+        # ax[0].set_ylabel("cumul interpolated")
+        #
+        # ax[1].plot(list(range(len(self.cumul_interpolated)-1)), self.resampled_proba_vol)
+        # ax[1].set_xlabel("n° point")
+        # ax[1].set_ylabel("diff cumul interpolated")
+        #
+        # plt.show()
     
     def demodulate(self):
         """This function demodulates the interpolated probability distribution"""
         # Sampling the size classes. Need to normalize the probability to get sum = 1.0
+        # print("**************************")
+        # for i in range(self.resampled_proba_vol.shape[0]):
+        #     print(f"vol = {self.resampled_proba_vol[i]}, sum = {np.sum(self.resampled_proba_vol)}, "
+        #           f" p = {self.resampled_proba_vol[i] / np.sum(self.resampled_proba_vol)}")
+        # print("**************************")
+
         ech_demod = np.random.choice(self.sizes_log_resampled[:-1], (10000,1), 
                                      p=(self.resampled_proba_vol / np.sum(self.resampled_proba_vol)))
         # fitting models from 1 to 10 components
@@ -182,7 +200,7 @@ class demodul_granulo:
     def plot_interpolation(self):
         """This method plots the cumulative probability interpolation"""
         # ------- Plotting interpolation----------- #
-        fig,(ax)=plt.subplots(1, 2, figsize = (11,4), dpi=100)
+        fig, ax = plt.subplots(1, 2, figsize = (11,4), dpi=100)
 
         # Data
         ax[0].plot(np.log(self.size_classes), self.proba_cumul_vol, linestyle='', marker='+')

Model/acoustic_data_loader.py

@@ -7,7 +7,7 @@ from matplotlib.colors import LogNorm
 
 # path_BS_raw_data = "/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/" \
 #                         "Data/Acoustic_data/20180107123500.aqa"
-# path_noise_data = "/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/" \
+# path_BS_raw_data = "/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/" \
 #                        "Data/AcousticNoise_data/20180107121600.aqa"
 
 
@@ -15,7 +15,7 @@ class AcousticDataLoader:
 
     def __init__(self, path_BS_raw_data: str):
 
-        self.path_BS_raw_data = path_BS_raw_data
+        # self.path_BS_raw_data = path_BS_raw_data
 
         # --- Load Backscatter acoustic raw data with RawAquascatData class ---
 
@@ -65,11 +65,13 @@ class AcousticDataLoader:
         # print(np.where((self._time) == 155)[0][0])
 
         # fig, ax = plt.subplots(nrows=1, ncols=1)
-        # ax.pcolormesh(self._time[0, :2200], -self._r[0, :], (self._BS_raw_data[0, :, :2200]),
-        #               cmap='viridis',
-        #               norm=LogNorm(vmin=1e-5, vmax=np.max(self._BS_raw_data[0, :, :2200])))  # , shading='gouraud')
-        # ax.pcolormesh(range(self._BS_raw_data.shape[2]), range(self._BS_raw_data.shape[0]), self._BS_raw_data[:, 1, :], cmap='viridis',
+        # # ax.pcolormesh(self._time[0, :2200], -self._r[0, :], (self._BS_raw_data[0, :, :2200]),
+        # #               cmap='viridis',
+        # #               norm=LogNorm(vmin=1e-5, vmax=np.max(self._BS_raw_data[0, :, :2200])))  # , shading='gouraud')
+        # ax.pcolormesh(range(self._BS_raw_data.shape[2]), range(self._BS_raw_data.shape[1]), self._BS_raw_data[2, :, :], cmap='viridis',
         #               norm=LogNorm(vmin=1e-5, vmax=np.max(self._BS_raw_data[:, 0, :])))  # , shading='gouraud')
+        # ax.set_xticks([])
+        # ax.set_yticks([])
         # plt.show()
 
         # --- Plot vertical profile for bottom detection ---

Model/acoustic_data_loader_UBSediFlow.py

@@ -30,7 +30,10 @@ from Model.udt_extract.raw_extract import raw_extract
     #                 "Rhone_20210519/Rhone_20210519/record/")
 
 # path_BS_raw_data0 = ("/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/Data/Raw_data_udt/")
-# filename0 = "raw_20210519_130000.udt"
+# filename0 = "raw_20210519_130643.udt"
+
+# path_BS_raw_data0 = ("/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/Data/Raw_data_udt/")
+# filename0 = "raw_20210520_085958.udt"
 
 # filename = "raw_20210519_115128.udt"
     # "raw_20210526_153310.udt"
@@ -177,12 +180,15 @@ class AcousticDataLoaderUBSediFlow:
                 if config == 1:
                     BS_data = np.array([data_us_dicts[config][channel]['echo_avg_profile']['data'][0]])
                     # print("BS_data shape ", BS_data.shape)
-                    print("******************************")
-                    date_list = [np.abs(datetime.datetime(2021, 5, 19, 13, 10, 00).timestamp()
-                          - date.timestamp()) for date in data_us_dicts[config][channel]['echo_avg_profile']['time']]
-                    print(np.where(date_list == np.min(date_list)))
-                    # == datetime.datetime(2021, 5, 19, 14, 10, 2, 644000))
-                    print("******************************")
+                    # print("******************************")
+                    # date_list = [np.abs(datetime.datetime(2021, 5, 19, 14, 10, 00).timestamp()
+                    #       - date.timestamp()) for date in data_us_dicts[config][channel]['echo_avg_profile']['time']]
+                    # print(date_list)
+                    # print(np.where(date_list == np.min(date_list)))
+                    # print((data_us_dicts[config][channel]['echo_avg_profile']['time'][np.where(date_list == np.min(date_list))[0][0]] -
+                    #       data_us_dicts[config][channel]['echo_avg_profile']['time'][0]).total_seconds())
+                    # # == datetime.datetime(2021, 5, 19, 14, 10, 2, 644000))
+                    # print("******************************")
 
                     for i in range(self._time.shape[1]):
                         BS_data = np.append(BS_data,
@@ -282,6 +288,11 @@ class AcousticDataLoaderUBSediFlow:
 
         self._time = self._time[:, :self._BS_raw_data.shape[2]]
         print("self._time.shape ", self._time.shape)
+        # print(f"time : {self._time}")
+        # print("****************************")
+        # for i in range(len(self._time[0, :])-1):
+        #     print(self._time[0, i+1] - self._time[0, i])
+        # print("****************************")
 
         print("self._r.shape ", self._r.shape)
 
@@ -345,7 +356,7 @@ class AcousticDataLoaderUBSediFlow:
         # --- Plot Backscatter US data ---
 
         # fig, ax = plt.subplots(nrows=1, ncols=1, layout="constrained")
-        # pcm = ax.pcolormesh(self._time[0, :], -self._r[0, 2:], np.log(self._BS_raw_data[0, 2:, :]),
+        # pcm = ax.pcolormesh(self._time[0, :], -self._r[0, :], np.log(self._BS_raw_data[0, :, :]),
         #                            cmap='plasma')#, shading='gouraud')
         # # pcm = ax.pcolormesh(list(range(self._BS_raw_data.shape[2])), list(range(self._BS_raw_data.shape[1])),
         # #                     np.log(self._BS_raw_data[0, :, :]),

Model/acoustic_inversion_method_high_concentration.py

@@ -83,10 +83,10 @@ class AcousticInversionMethodHighConcentration():
         x = k * a
         f = (x ** 2 * (1 - 0.25 * np.exp(-((x - 1.5) / 0.35) ** 2)) * (1 + 0.6 * np.exp(-((x - 2.9) / 1.15) ** 2))) / (
                     42 + 28 * x ** 2)
-        print(f"form factor = {f}")
+        # print(f"form factor = {f}")
         return f
 
-    def ks(self, freq, pdf):
+    def ks(self, num_sample, freq, pdf):
         # --- Calcul de la fonction de form ---
         # form_factor = self.form_factor_function_MoateThorne2012(a, freq)
         # print(f"form_factor shape = {form_factor}")
@@ -95,7 +95,7 @@ class AcousticInversionMethodHighConcentration():
         #--- Particle size distribution ---
         proba_num = (
             self.compute_particle_size_distribution_in_number_of_particles(
-                num_sample=2, r_grain=stg.radius_grain, frac_vol_cumul=stg.frac_vol_sand_cumul))
+                num_sample=num_sample, r_grain=stg.radius_grain, frac_vol_cumul=stg.frac_vol_sand_cumul))
 
         print(f"proba_num : {proba_num}")
 
@@ -228,9 +228,13 @@ class AcousticInversionMethodHighConcentration():
             # print(f"range_cells = {range_cells}")
             loc = (range_cells >= r_ini) * (range_cells < r_end)
             # print(f"loc = {loc}")
+            # print(f"loc shape = {len(loc)}")
+
             # Filling the array with interpolation values
             res[loc] = range_cells[loc] * a + b
+            # print(res.shape)
             # print(f"res = {res}")
+            # print(f"1. res.shape = {res.shape}")
 
         # Filling first and last values
         i = 0
@@ -247,20 +251,24 @@ class AcousticInversionMethodHighConcentration():
         if stg.r_bottom.size != 0:
             res[np.where(range_cells > r_bottom)] = np.nan
 
-        # print(f"res.shape = {res.shape}")
+        # print(f"2. res.shape = {res.shape}")
         # print(f"res = {res}")
 
         loc_point_lin_interp0 = range_cells[np.where((range_cells > sample_depth[0]) & (range_cells < sample_depth[-1]))]
         # print(f"range_cells : {range_cells}")
+        # print(f"loc_point_lin_interp0 shape : {len(loc_point_lin_interp0)}")
         # print(f"loc_point_lin_interp0 : {loc_point_lin_interp0}")
+        res0 = res[np.where((range_cells > sample_depth[0]) & (range_cells < sample_depth[-1]))]
 
         loc_point_lin_interp = loc_point_lin_interp0[np.where(loc_point_lin_interp0 > range_cells[0])]
+        # print(f"loc_point_lin_interp shape : {len(loc_point_lin_interp)}")
         # print(f"loc_point_lin_interp : {loc_point_lin_interp}")
+        res = res0[np.where(loc_point_lin_interp0 > range_cells[0])]
 
-        # fig, ax = plt.subplots(nrows=1, ncols=1)
-        # ax.plot(loc_point_lin_interp, res[:len(loc_point_lin_interp)], c="blue")
-        # ax.plot(sample_depth, M_profile, color="k", marker="o")
-        # plt.show()
+        fig, ax = plt.subplots(nrows=1, ncols=1)
+        ax.plot(loc_point_lin_interp, res[:len(loc_point_lin_interp)], marker="*", mfc="blue")
+        ax.plot(sample_depth, M_profile, marker="o", mfc="k", mec="k")
+        plt.show()
 
         return (loc_point_lin_interp, res)
 
@@ -270,6 +278,7 @@ class AcousticInversionMethodHighConcentration():
         delta_r = r[1] - r[0]
         zeta = alpha_s / (np.sum(M_profile_fine*delta_r))
 
+        print(f"np.sum(M_profile_fine*delta_r) : {np.sum(M_profile_fine*delta_r)}")
         # zeta0 = np.array([0.021, 0.035, 0.057, 0.229])
         # zeta = zeta0[ind]
         # zeta0 = np.array([0.04341525, 0.04832906, 0.0847188, np.nan])
@@ -335,6 +344,9 @@ class AcousticInversionMethodHighConcentration():
                          (freq2 ** X))) /
                   (zeta_freq2 - zeta_freq1))
 
+        # logVBI = (freq2**2 * np.log(j_cross_section_freq1 / freq1**X) -
+        #           freq1**2 * np.log(j_cross_section_freq2 / freq2**X)) / (freq2**2 - freq1**2)
+
         # logVBI = (( np.full((stg.r.shape[1], stg.t.shape[1]), zeta_freq2) *
         #            np.log(j_cross_section_freq1 * np.exp(4 * r2D * np.full((stg.r.shape[1], stg.t.shape[1]), water_attenuation_freq1)) /
         #                   (freq1 ** X)) -
@@ -348,8 +360,8 @@ class AcousticInversionMethodHighConcentration():
         return np.exp(logVBI)
 
     # ------------- Computing SSC fine ------------- #
-    def SSC_fine(self, zeta, r2D, VBI, freq, X, j_cross_section):
-        SSC_fine = (1/zeta) * ( 1/(4 * r2D) * np.log((VBI * freq**X) / j_cross_section) )
+    def SSC_fine(self, zeta, r2D, VBI, freq, X, j_cross_section, alpha_w):
+        SSC_fine = (1/zeta) * ( 1/(4 * r2D) * np.log((VBI * freq**X) / j_cross_section) - alpha_w)
         print("compute SSC fine finished")
         return SSC_fine
 

Model/granulo_loader.py

@@ -42,9 +42,8 @@ class GranuloLoader:
         # ==============================================================================================================
         # ==============================================================================================================
 
-#         print(self._r_grain.shape)
-#         N_sample = 2
-#
+#         N_sample = 0
+# #
 #         fig, ax = plt.subplots(1, 2)
 #         ax[0].plot(self._r_grain, self._frac_vol[N_sample, :], color="k", marker='.')
 #         ax[0].set_xscale('log')
@@ -57,17 +56,20 @@ class GranuloLoader:
 #
 #         plt.show()
 #
+#         print(f"self._r_grain.shape : {self._r_grain.shape}")
 #         print(f"self._r_grain : {self._r_grain}")
+#         print(f"self._frac_vol_cumul.shape : {self._frac_vol_cumul[N_sample, :].shape}")
 #         print(f"self._frac_vol_cumul[N_sample, :] : {self._frac_vol_cumul[N_sample, :]}")
-#
+#         print(np.where(self._frac_vol_cumul[N_sample, :] > 0))
+# #
 #         min_demodul = 1e-6
 #         max_demodul = 500e-6
-#         sample_demodul = demodul_granulo(self._r_grain,
+#         sample_demodul = demodul_granulo(self._r_grain[:],
 #                                               self._frac_vol_cumul[N_sample, :],
 #                                               min_demodul, max_demodul)
 #
 #         print(f"sample_demodul : {sample_demodul.demodul_data_list}")
-#
+
 #         N_modes = 3
 #         sample_demodul.print_mode_data(N_modes)
 #         sample_demodul.plot_interpolation()
@@ -120,10 +122,14 @@ class GranuloLoader:
 
 
 # if __name__ == "__main__":
-#     # GranuloLoader("/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/Data/Granulo_data/"
-#     #               "fine_sample_file.ods")
-#     GranuloLoader("/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/Data/Granulo_data/"
-#                   "sand_sample_file.ods")
+    # GranuloLoader("/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/Data/Granulo_data/"
+    #               "fine_sample_file.ods")
+    # GranuloLoader("/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/Data/Granulo_data/"
+    #               "sand_sample_file.ods")
+    # GranuloLoader("/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/Data/"
+    #               "pt_acoused_cnr_7nov2023/echantil/fine_new_file_acoused_101RDS 3.ods")
+    # GranuloLoader("/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/Data/"
+    #               "pt_acoused_cnr_7nov2023/echantil/sand_new_file_acoused_101RDS 3.ods")