% 定义水力学计算参数
Q = 15e-3;    % Flow rate (m^3/s)
W = 0.6;      % Channel width (m)
h = 0.0863;   % Flow depth (m)
nu = 0.89e-6; % Kinematic viscosity of water (m^2/s)
ks = 0.01;    % Bed roughness (m)

% 计算流速
U = Q / (W * h); % Flow velocity (m/s)
fprintf('Flow velocity: %.6f m/s\n', U);

% 初始化参数
slope = 0.01;    % Initial channel slope
ustar = sqrt(9.81 * h * slope); % Initial shear velocity (m/s)

% 迭代求解剪切流速
max_iter = 10000;
tolerance = 1e-6;
ustar_new = 0;

for i = 1:max_iter
    delta_nu = nu / ustar;
    % 计算 y0 (yo 函数)
    y0 = 0.11 * delta_nu + 0.033 * ks * (1 - exp(-0.04 * ks / delta_nu));
    % 计算新的剪切流速 (e 在 MATLAB 中是 exp(1))
    ustar_new = U / (2.5 * log(h / y0 / exp(1)));
    
    % 检查收敛条件
    if abs(ustar_new - ustar) < tolerance
        break;
    end
    ustar = ustar_new;
end

% 输出结果
fprintf('Shear velocity: %.6f m/s\n', ustar);
% 重新计算渠道坡度
slope = (ustar^2) / (9.81 * h);
fprintf('Channel slope: %.6f * 1e-3\n', slope * 1e3);

% 计算 iarg
iarg = 11950 * (-slope) + 59.75;
fprintf('Flume input arg: %.6f\n', iarg);
