import matplotlib.pyplot as plt import numpy as np import pandas as pd from Model.GrainSizeTools import demodul_granulo, mix_gaussian_model class GranuloLoader: """ This class allows to load granulo data file """ def __init__(self, path: str): self._path = path self._data = pd.read_excel(self._path, engine="odf", header=0) self._time = self._data.iloc[:, 0].tolist() self._y = self._data.iloc[:, 1].tolist() # distance from left bank (m) self._z = self._data.iloc[:, 2].tolist() # depth (m) self._r_grain = 1e-6 * self._data.columns.values[5:] / 2 # self._r_grain = 1e-6 * np.array(self._data.columns.values)[5:].astype(float) / 2 # grain radius (m) self._Ctot = self._data.iloc[:, 3].tolist() # Total concentration (g/L) self._D50 = self._data.iloc[:, 4].tolist() # median diameter (um) self._frac_vol = np.array(self._data.iloc[:, 5:]) / 100 # Volume fraction (%) self._frac_vol_cumul = np.cumsum(self._frac_vol, axis=1) # Cumulated volume fraction (%) # print(type(self._frac_vol_cumul), self._frac_vol_cumul.shape, self._frac_vol_cumul) # # --- Load sand sediments data file --- # self.path_sand = path_sand # self._data_sand = pd.read_excel(self.path_sand, engine="odf", header=0) # # self._Ctot_sand = np.array(self._data_sand.iloc[:, 2]) # Total concentration (g/L) # self._D50_sand = np.array(self._data_sand.iloc[:, 3]) # median diameter (um) # self._frac_vol_sand = np.array(self._data_sand.iloc[:, 4:]) # Volume fraction (%) # # self._frac_vol_sand_cumul = np.cumsum(self._frac_vol_sand, axis=1) # Cumulated volume fraction (%) # # # --- Compute % of fine and % of sand sediment in total concentration --- # # self._Ctot_fine_per_cent = 100 * self._Ctot_fine / (self._Ctot_fine + self._Ctot_sand) # self._Ctot_sand_per_cent = 100 * self._Ctot_sand / (self._Ctot_fine + self._Ctot_sand) # ============================================================================================================== # ============================================================================================================== # 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') # ax[0].set_xlabel('Radius ($\mu m$)') # ax[0].set_ylabel('Class size volume fraction') # # ax[1].plot([self._r_grain[i+1]-self._r_grain[i] for i in range(self._r_grain.shape[0]-1)], list(range(self._r_grain.shape[0]-1)), color="k", marker="x") # ax[1].set_xlabel('Ecart inter-class') # ax[1].set_ylabel('n° échantillon') # # 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[:], # 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() # sample_demodul.plot_modes(N_modes) # # print(f"mu_list : {sample_demodul.demodul_data_list[3 - 1].mu_list}") # print(f"sigma_list : {sample_demodul.demodul_data_list[3 - 1].sigma_list}") # print(f"w_list : {sample_demodul.demodul_data_list[3 - 1].w_list}") # # resampled_log_array = np.log(np.logspace(-10, -2, 3000)) # proba_vol_demodul = mix_gaussian_model(resampled_log_array, # sample_demodul.demodul_data_list[2].mu_list, # sample_demodul.demodul_data_list[2].sigma_list, # sample_demodul.demodul_data_list[2].w_list) # # # # proba_vol_demodul = proba_vol_demodul / np.sum(proba_vol_demodul) # ss = np.sum(proba_vol_demodul / np.exp(resampled_log_array) ** 3) # proba_num = proba_vol_demodul / np.exp(resampled_log_array) ** 3 / ss # # print(f"proba_num : {proba_num}") # freq = 5e6 # # a2f2pdf = 0 # a3pdf = 0 # for i in range(len(resampled_log_array)): # a = np.exp(resampled_log_array)[i] # a2f2pdf += a**2 * form_factor_function_MoateThorne2012(a, freq)**2 * proba_num[i] # a3pdf += a**3 * proba_num[i] # # print(f"a2f2pdf = {a2f2pdf}") # print(f"a3pdf = {a3pdf}") # # ks = (a2f2pdf / a3pdf) # # print(f"ks = {ks}") # # # def form_factor_function_MoateThorne2012(a_s, freq, C=1500): # """This function computes the form factor based on the equation of # Moate and Thorne (2012)""" # # computing the wave number # k = 2 * np.pi * freq / C # x = k * a_s # 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) # return f # if __name__ == "__main__": # GranuloLoader("/home/bmoudjed/Documents/2 Data/Confluence_Rhône_Isere_2018/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") # # form_factor = form_factor_function_MoateThorne2012(GranuloLoader.) # # def ks(a_s, rho_s, freq, pdf): # # --- Calcul de la fonction de form --- # form_factor = form_factor_function_MoateThorne2012(a_s, freq) # print(f"form_factor shape = {len(form_factor)}") # # print(f"form_factor = {form_factor}") # # # --- Gaussian mixture --- # # sample_demodul = demodul_granulo(a_s.astype(float), pdf, 0.17e-6, 200e-6) # # sample_demodul.plot_interpolation() # # # ss = np.sum(pdf / a_s ** 3) # # proba_num = (pdf / a_s ** 3) / ss # # # --- Compute k_s by dividing two integrals --- # a2f2pdf = 0 # a3pdf = 0 # for i in range(len(pdf)): # a2f2pdf += a_s[i] ** 2 * form_factor[i] * proba_num[i] # a3pdf += a_s[i] ** 3 * proba_num[i] # # ks = np.sqrt(a2f2pdf / (rho_s * a3pdf)) # # # ks = np.array([0.04452077, 0.11415143, 0.35533713, 2.47960051]) # # ks = ks0[ind] # return ks