442 lines
20 KiB
Python
442 lines
20 KiB
Python
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
import settings as stg
|
|
from Model.GrainSizeTools import demodul_granulo, mix_gaussian_model
|
|
|
|
|
|
class AcousticInversionMethodHighConcentration():
|
|
|
|
""" Thi class compute acoustic inversion method adapted for high suspended sediments concentration
|
|
For instance, the case of dam flush downstream Rhone-Isere confluence (07/01/2018) evaluated at ~ 10g/L """
|
|
|
|
def __init__(self):
|
|
|
|
pass
|
|
|
|
# ==========================================
|
|
# Functions
|
|
# ==========================================
|
|
|
|
# ---------- Computing sound speed ------------- #
|
|
def water_velocity(self, T):
|
|
"""Computing sond speed from Bilaniuk and Wong 1993"""
|
|
C = 1.40238744 * 1e3 + 5.03836171 * T - 5.81172916 * 1e-2 * T ** 2 + 3.34638117 * 1e-4 * T ** 3 - \
|
|
1.48259672 * 1e-6 * T ** 4 + 3.16585020 * 1e-9 * T ** 5
|
|
return C
|
|
|
|
# -------- Computing water attenuation coefficient ----------- #
|
|
def water_attenuation(self, freq, T):
|
|
"""Computing attenuation from François and Garrison 1982"""
|
|
if T > 20:
|
|
alpha = (3.964 * 1e-4 - 1.146 * 1e-5 * T + 1.45 * 1e-7 * T ** 2 - 6.5 * 1e-10 * T ** 3) * 1e-3 * \
|
|
(np.log(10) / 20) * (freq * 1e-3) ** 2
|
|
else:
|
|
alpha = (4.937 * 1e-4 - 2.59 * 1e-5 * T + 9.11 * 1e-7 * T ** 2 - 1.5 * 1e-8 * T ** 3) * 1e-3 * \
|
|
(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
|
|
max_demodul = 500e-6
|
|
|
|
sample_demodul = demodul_granulo(r_grain, frac_vol_cumul[num_sample], min_demodul, max_demodul)
|
|
|
|
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)
|
|
|
|
# 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
|
|
|
|
return proba_num
|
|
|
|
# ------------- Computing ks ------------- #
|
|
def form_factor_function_MoateThorne2012(self, a, freq, C):
|
|
"""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
|
|
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
|
|
a3pdf = 0
|
|
for i in range(len(resampled_log_array)):
|
|
a = np.exp(resampled_log_array)[i]
|
|
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
|
|
|
|
# -------------- Computing Kt -------------- #
|
|
def kt_corrected(self, r, water_velocity, RxGain, TxGain, kt_ref):
|
|
"""Computing the instrument constant Kt that depends on gain and temperature"""
|
|
# Cell size
|
|
delta_r = r[1] - r[0]
|
|
# Pulse length
|
|
tau = 2 * delta_r / 1500
|
|
# Sound speed
|
|
cel = water_velocity
|
|
# Reference pulse length
|
|
tau_ref = 13.33333e-6
|
|
# Reference sound speed
|
|
c_ref = 1475
|
|
# Gain
|
|
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)
|
|
|
|
def M_profile_SCC_fine_interpolated(self, sample_depth, M_profile, range_cells, r_bottom):
|
|
res = np.zeros((len(range_cells),)) * np.nan
|
|
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]):
|
|
res[i] = c1[0]
|
|
i += 1
|
|
|
|
# Filling the last values
|
|
i = -1
|
|
while np.isnan(res[i]):
|
|
res[i] = c1[-1]
|
|
i += -1
|
|
|
|
if r_bottom.size != 0:
|
|
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]))]
|
|
res0 = res[np.where((range_cells > l2[0]) & (range_cells < l2[-1]))]
|
|
|
|
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 ------------- #
|
|
def zeta(self, alpha_s, r, M_profile_fine):
|
|
|
|
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 ------------- #
|
|
def VBI_cross_section(self, freq1, freq2,
|
|
zeta_freq1, zeta_freq2,
|
|
j_cross_section_freq1, j_cross_section_freq2,
|
|
r2D,
|
|
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)) -
|
|
zeta_freq1 *
|
|
np.log(j_cross_section_freq2 * np.exp(4 * r2D * water_attenuation_freq2) /
|
|
(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
|
|
|
|
|