From 65a9422f711a9bc8622421d0f8f10b5bc4a9c10d Mon Sep 17 00:00:00 2001 From: brahim Date: Thu, 6 Mar 2025 10:23:27 +0100 Subject: [PATCH] acoustic inversion method high concentration file is cleaned from useless commented lines and useless print --- ...tic_inversion_method_high_concentration.py | 247 +----------------- 1 file changed, 3 insertions(+), 244 deletions(-) diff --git a/Model/acoustic_inversion_method_high_concentration.py b/Model/acoustic_inversion_method_high_concentration.py index 83482b3..b2a2cbd 100644 --- a/Model/acoustic_inversion_method_high_concentration.py +++ b/Model/acoustic_inversion_method_high_concentration.py @@ -21,7 +21,6 @@ # -*- coding: utf-8 -*- -import matplotlib.pyplot as plt import numpy as np import settings as stg from Model.GrainSizeTools import demodul_granulo, mix_gaussian_model @@ -58,17 +57,6 @@ class AcousticInversionMethodHighConcentration(): (np.log(10) / 20) * (freq * 1e-3) ** 2 return alpha - # ---------- Conmpute FBC ---------- - # def compute_FCB(self): - # # print(self.BS_averaged_cross_section_corr.V.shape) - # # print(self.r_2D.shape) - # FCB = np.zeros((256, 4, 1912)) - # for f in range(4): - # # print(self.alpha_w_function(self.Freq[f], self.temperature)) - # FCB[:, f, :] = np.log(self.BS_averaged_cross_section_corr.V[:, f, :]) + np.log(self.r_3D[:, f, :]) + \ - # np.log(2 * self.alpha_w_function(self.Freq[f], self.temperature) * self.r_3D[:, f, :]) - # return FCB - # --- Gaussian mixture --- def compute_particle_size_distribution_in_number_of_particles(self, num_sample, r_grain, frac_vol_cumul): min_demodul = 1e-6 @@ -82,15 +70,6 @@ class AcousticInversionMethodHighConcentration(): sample_demodul.demodul_data_list[2].sigma_list, sample_demodul.demodul_data_list[2].w_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}") - 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 @@ -106,23 +85,9 @@ 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}") return f - # def ks(self, num_sample_sand, radius_grain_sand, frac_vol_sand_cumul, freq, C): def ks(self, proba_num, freq, C): - # --- Calcul de la fonction de form --- - # form_factor = self.form_factor_function_MoateThorne2012(a, freq) - # print(f"form_factor shape = {form_factor}") - # print(f"form_factor = {form_factor}") - - #--- Particle size distribution --- - # proba_num = ( - # self.compute_particle_size_distribution_in_number_of_particles( - # num_sample=num_sample_sand, r_grain=radius_grain_sand, frac_vol_cumul=frac_vol_sand_cumul[num_sample_sand])) - - # print(f"proba_num : {proba_num}") - # --- Compute k_s by dividing two integrals --- resampled_log_array = np.log(np.logspace(-10, -2, 3000)) a2f2pdf = 0 @@ -132,28 +97,17 @@ class AcousticInversionMethodHighConcentration(): a2f2pdf += a**2 * self.form_factor_function_MoateThorne2012(a, freq, C)**2 * proba_num[i] a3pdf += a**3 * proba_num[i] - # print("form factor ", self.form_factor_function_MoateThorne2012(a, freq, C)) - # print(f"a2f2pdf = {a2f2pdf}") - # print(f"a3pdf = {a3pdf}") - ks = np.sqrt(a2f2pdf / a3pdf) - # ks = np.array([0.04452077, 0.11415143, 0.35533713, 2.47960051]) - # ks = ks0[ind] return ks # ------------- Computing sv ------------- # def sv(self, ks, M_sand): - # print(f"ks = {ks}") - # print(f"M_sand = {M_sand}") sv = (3 / (16 * np.pi)) * (ks ** 2) * M_sand - # sv = np.full((stg.r.shape[1], stg.t.shape[1]), sv0) return sv # ------------- Computing X ------------- # def X_exponent(self, freq1, freq2, sv_freq1, sv_freq2): - # X0 = [3.450428714146802, 3.276478927777019, 3.6864638665972893, 0] - # X = X0[ind] X = np.log(sv_freq1 / sv_freq2) / np.log(freq1 / freq2) return X @@ -174,165 +128,43 @@ class AcousticInversionMethodHighConcentration(): gain = 10 ** ((RxGain + TxGain) / 20) # Computing Kt kt = kt_ref * gain * np.sqrt(tau * cel / (tau_ref * c_ref)) # 1D numpy array - # kt = np.reshape(kt0, (1, 2)) # convert to 2d numpy array to compute J_cross_section - # print(f"kt = {kt}") - # kt_2D = np.repeat(np.array([kt]), stg.r.shape[1], axis=0) - # print("kt 2D ", kt_2D) - # print("kt 2D shape ", kt_2D.shape) - # # kt_3D = np.zeros((kt_2D.shape[1], kt_2D.shape[0], stg.t.shape[1])) - # # for k in range(kt_2D.shape[1]): - # # kt_3D[k, :, :] = np.repeat(kt_2D, stg.t.shape[1], axis=1)[:, k * stg.t.shape[1]:(k + 1) * stg.t.shape[1]] - # kt_3D = np.repeat(kt_2D.transpose()[:, :, np.newaxis], stg.t.shape[1], axis=2) - # # print("kt 3D ", kt_3D) - # print("kt 3D shape ", kt_3D.shape) - return kt # ------------- Computing J_cross_section ------------- # def j_cross_section(self, BS, r2D, kt): - # J_cross_section = np.zeros((1, BS.shape[1], BS.shape[2])) # 2 because it's a pair of frequencies - # print("BS.shape", BS.shape) - # print("r2D.shape", r2D.shape) - # print("kt.shape", kt.shape) - - # if stg.ABS_name == "Aquascat 1000R": - # print("--------------------------------") - # print("BS : ", BS) - # print("BS min : ", np.nanmin(BS)) - # print("BS max : ", np.nanmax(BS)) - # print("r2D : ", r2D) - # print("kt shape : ", kt.shape) - # print("kt : ", kt) - # print("--------------------------------") - # for k in range(1): - # J_cross_section[k, :, :] = (3 / (16 * np.pi)) * ((BS[k, :, :]**2 * r2D[k, :, :]**2) / kt[k, :, :]**2) J_cross_section = (3 / (16 * np.pi)) * ((BS**2 * r2D**2) / kt**2) - # J_cross_section[J_cross_section == 0] = np.nan - # print("J_cross_section.shape", J_cross_section.shape) - # elif stg.ABS_name == "UB-SediFlow": - # for k in range(1): - # J_cross_section[k, :, :] = (3 / (16 * np.pi)) * ((BS[k, :, :]**2 * r2D[0, :, :]**2) / kt[k, :, :]**2) - # print("compute j_cross_section finished") return J_cross_section # ------------- Computing alpha_s ------------- # def alpha_s(self, sv, j_cross_section, depth, alpha_w): alpha_s = (np.log(sv / j_cross_section) / (4 * depth)) - alpha_w - print("----------------------------") - print(f"sv = {sv}") - print(f"j_cross_section = {j_cross_section}") - print(f"depth = {depth}") - print(f"alpha_w = {alpha_w}") - print(f"(np.log(sv / j_cross_section) / (4 * depth)) = {(np.log(sv / j_cross_section) / (4 * depth))}") - print(f"alpha_s {alpha_s}") - return alpha_s - # ------------- Computing interpolation of fine SSC data obtained from water sampling ------------- - # ------------- collected at various depth in the vertical sample ------------- - # def M_profile_SCC_fine_interpolated(self, sample_depth, M_profile, range_cells, r_bottom): - # res = np.zeros((len(range_cells),)) * np.nan - # for i in range(len(M_profile) - 1): - # # print(f"i = {i}") - # r_ini = sample_depth[i] - # # print(f"r_ini = {r_ini}") - # c_ini = M_profile[i] - # # print(f"c_ini = {c_ini}") - # r_end = sample_depth[i + 1] - # # print(f"r_end = {r_end}") - # c_end = M_profile[i + 1] - # # print(f"c_end = {c_end}") - # - # # Computing the linear equation - # a = (c_end - c_ini) / (r_end - r_ini) - # # print(f"a = {a}") - # b = c_ini - a * r_ini - # # print(f"b = {b}") - # - # # Finding the indices of r_ini and r_end in the interpolated array - # # 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 - # while np.isnan(res[i]): - # res[i] = M_profile[0] - # i += 1 - # - # # Filling the last values - # i = -1 - # while np.isnan(res[i]): - # res[i] = M_profile[-1] - # i += -1 - # # print(f"res.shape = {res.shape}") - # # print(f"res = {res}") - # # print(f"r_bottom.shape = {r_bottom.shape}") - # # print(f" = {res}") - # - # if r_bottom.shape != (0,): - # res[np.where(range_cells > r_bottom)] = np.nan - # - # 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)], marker="*", mfc="blue") - # # ax.plot(sample_depth, M_profile, marker="o", mfc="k", mec="k") - # # plt.show() - # - # return (loc_point_lin_interp, res) + # ------------- Computing interpolation of fine SSC ------------- def M_profile_SCC_fine_interpolated(self, sample_depth, M_profile, range_cells, r_bottom): + '''Computing interpolation of fine SSC data obtained from water sampling + collected at various depth in the vertical sample''' res = np.zeros((len(range_cells),)) * np.nan - print("range_cells ", range_cells.shape) l0 = sample_depth - print("l0 = ", l0) l1 = [l0.index(x) for x in sorted(l0)] - print("l1 = ", l1) l2 = [l0[k] for k in l1] - print("l2 = ", l2) c1 = [list(M_profile)[j] for j in l1] - print("c1 = ", c1) for i in range(len(c1) - 1): # print("i = ", i) r_ini = l2[i] c_ini = c1[i] r_end = l2[i + 1] c_end = c1[i + 1] - print("r_ini ", r_ini, "c_ini ", c_ini, "r_end ", r_end, "c_end ", c_end) # Computing the linear equation a = (c_end - c_ini) / (r_end - r_ini) b = c_ini - a * r_ini - print("range_cells ", (range_cells)) # Finding the indices of r_ini and r_end in the interpolated array loc = (range_cells >= r_ini) * (range_cells < r_end) - print("range_cells >= r_ini ", range_cells >= r_ini) - print("range_cells < r_end ", range_cells < r_end) - print("loc ", loc) # Filling the array with interpolation values res[loc] = range_cells[loc] * a + b - print("a = ", a, "b = ", b) - - print("res ", res) - # Filling first and last values i = 0 while np.isnan(res[i]): @@ -346,9 +178,6 @@ class AcousticInversionMethodHighConcentration(): i += -1 if r_bottom.size != 0: - print("res ", res.shape) - print("range_cells ", len(range_cells)) - # print("r_bottom ", len(r_bottom)) res[np.where(range_cells > r_bottom)] = np.nan loc_point_lin_interp0 = range_cells[np.where((range_cells > l2[0]) & (range_cells < l2[-1]))] @@ -357,13 +186,6 @@ class AcousticInversionMethodHighConcentration(): loc_point_lin_interp = loc_point_lin_interp0[np.where(loc_point_lin_interp0 > l2[0])] res = res0[np.where(loc_point_lin_interp0 > l2[0])] - # fig, ax = plt.subplots(nrows=1, ncols=1) - # ax.plot(res[:len(loc_point_lin_interp)], -loc_point_lin_interp, marker="*", mfc="blue") - # ax.plot(c1, [-x for x in l2], marker="o", mfc="k", mec="k", ls="None") - # ax.set_xlabel("Concentration (g/L)") - # ax.set_ylabel("Depth (m)") - # plt.show() - return (loc_point_lin_interp, res) # ------------- Computing zeta ------------- # @@ -372,39 +194,6 @@ class AcousticInversionMethodHighConcentration(): delta_r = r[1] - r[0] zeta = alpha_s / (np.sum(np.array(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]) - # zeta = zeta0[[ind1, ind2]] - - # for k in range(3): - # for p in range(3): - # if np.isnan(ind_X_min_around_sample[p, k]): - # zeta_list_exp.append(np.nan) - # else: - # ind_X_min = int(ind_X_min_around_sample[p, k]) - # ind_X_max = int(ind_X_max_around_sample[p, k]) - # ind_r_min = int(ind_r_min_around_sample[p, k]) - # ind_r_max = int(ind_r_max_around_sample[p, k]) - # - # R_temp = R_cross_section[ind_r_min:ind_r_max, :, ind_X_min:ind_X_max] - # J_temp = J_cross_section[ind_r_min:ind_r_max, :, ind_X_min:ind_X_max] - # aw_temp = aw_cross_section[ind_r_min:ind_r_max, :, ind_X_min:ind_X_max] - # sv_temp_1 = np.repeat([sv_list_temp[3 * k + p]], np.shape(R_temp)[0], axis=0) - # sv_temp = np.swapaxes(np.swapaxes(np.repeat([sv_temp_1], np.shape(R_temp)[2], axis=0), 1, 0), 2, 1) - # ind_depth = np.where(R_cross_section[:, 0, 0] >= M_list_temp[k][0, p + 1])[0][0] - # # Using concentration profile - # zeta_temp = alpha_s / ((1 / M_list_temp[k][0, p + 1]) * (R_cross_section[0, 0, 0] * M_list_temp[k][1, 0] + - # delta_r * np.sum(M_interpolate_list[k][:ind_depth]))) - # zeta_temp = (1 / (4 * R_temp) * - # np.log(sv_temp / J_temp) - aw_temp) / ((1 / M_list_temp[k][0, p + 1]) * - # (R_cross_section[0, 0, 0] * M_list_temp[k][ - # 1, 0] + - # delta_r * np.sum( - # M_interpolate_list[k][:ind_depth]))) - # zeta_list_exp.append(np.mean(np.mean(zeta_temp, axis=0), axis=1)) - return zeta # ------------- Computing VBI ------------- # @@ -415,21 +204,6 @@ class AcousticInversionMethodHighConcentration(): water_attenuation_freq1, water_attenuation_freq2, X): - # print('self.zeta_exp[ind_j].shape', self.zeta_exp[ind_j]) - # print('np.log(self.j_cross_section[:, ind_i, :]).shape', np.log(self.j_cross_section[:, ind_i, :]).shape) - # print('self.r_3D[:, ind_i, :]', self.r_3D[:, ind_i, :].shape) - # print('self.water_attenuation[ind_i]', self.water_attenuation[ind_i]) - # print('self.x_exp[0.3-1 MHz]', self.x_exp['0.3-1 MHz'].values[0]) - # print("start computing VBI") - # print("================================") - # print(f"zeta_freq2 : {zeta_freq2}") - # print(f"j_cross_section_freq1 : {j_cross_section_freq1.shape}") - # print(f"r2D : {r2D.shape}") - # print(f"water_attenuation_freq1 : {water_attenuation_freq1}") - # print(f"freq1 : {freq1}") - # print(f"X : {X}") - # print("================================") - logVBI = ((zeta_freq2 * np.log(j_cross_section_freq1 * np.exp(4 * r2D * water_attenuation_freq1) / (freq1 ** X)) - @@ -438,31 +212,16 @@ 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)) - - # np.full((stg.r.shape[1], stg.t.shape[1]), zeta_freq1) * - # np.log(j_cross_section_freq2 * np.exp(4 * r2D * np.full((stg.r.shape[1], stg.t.shape[1]), water_attenuation_freq2)) / - # (freq2 ** X))) / - # (zeta_freq2 - zeta_freq1)) - - print("compute VBI finished") - return np.exp(logVBI) # ------------- Computing SSC fine ------------- # 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 # ------------- Computing SSC sand ------------- # def SSC_sand(self, VBI, freq, X, ks): SSC_sand = (16 * np.pi * VBI * freq ** X) / (3 * ks**2) - print("compute SSC sand finished") return SSC_sand