%% 環境參數設置
h = 100; % 水深(m)
c_profile = [0,1500; 50,1520; 100,1540]; % 聲速剖面(c(z)=1500+20z/100)
bottom_type = 'mud'; % 海底類型

%% 聲源與接收器
src_depth = 50; % 聲源深度(m)
rx_depth = 50;  % 接收器深度(m)
range = 500;    % 傳播距離(m)

%% 射線追蹤
num_rays = 50;  % 射線數量
angles = linspace(-30,30,num_rays); % 發射角(°)
traces = cell(num_rays,1);

for i = 1:num_rays
    [path, loss, delay] = ray_trace(src_depth, rx_depth, range, angles(i), c_profile, bottom_type);
    traces{i} = struct('path',path, 'loss',loss, 'delay',delay);
end

%% 結果可視化
figure;
plot(range,path(:,1),'r',range,path(:,2),'b'); % 繪製聲線路徑
xlabel('傳播距離(m)'); ylabel('深度(m)');
title('聲線傳播路徑');

figure;
plot(0:0.1:range,10*log10(cumsum(loss))); % 繪製傳播損失曲線
xlabel('距離(m)'); ylabel('傳播損失(dB)');

function [path, loss, delay] = ray_trace(src_z, rx_z, R, theta, c_profile, bottom)
    % 初始化參數
    g = 9.81; rho = 1025; % 重力加速度、海水密度
    c0 = interp1(c_profile(:,1),c_profile(:,2),src_z); % 起始聲速
    dz = 0.1; % 垂直步長(m)
    
    % 聲線追蹤循環
    path = [src_z; 0]; % 起始點
    total_loss = 0;
    time = 0;
    z = src_z;
    
    while z > 0
        % 計算聲速梯度
        c = interp1(c_profile(:,1),c_profile(:,2),z);
        dz_step = dz * (c0/c)^(1/2); % 聲速梯度修正步長
        
        % 更新位置
        z = z - dz_step;
        path = [path; z];
        
        % 計算吸收損失
        alpha = 0.1 * (1 + 0.01*(c-1500)^2); % 吸收係數模型
        total_loss = total_loss + alpha * dz_step;
        
        % 時間累積
        time = time + dz_step / c;
    end
    
    % 邊界反射損失計算
    if z < 0
        reflection_loss = 10*log10(0.5*(1 + (c0/c_bottom)^2)); % 海底反射模型
        total_loss = total_loss + reflection_loss;
    end
    
    loss = total_loss;
    delay = time;
end