三层介质射线追踪MATLAB程序

三层介质射线追踪MATLAB程序。这个程序可以模拟射线在三种不同介质中的传播路径，包括折射、反射和全反射现象。

基础物理模型

考虑三层水平分层的介质，每层具有不同的折射率：

• 第一层：折射率 n₁，厚度 d₁
• 第二层：折射率 n₂，厚度 d₂
• 第三层：折射率 n₃，半无限空间
• 射线：从第一层顶部以角度θ₁入射

%% 三层介质射线追踪 - 基础版本
clear; close all; clc;

%% 1. 参数设置
% 介质参数
n = [1.0, 1.5, 2.0];     % 三层介质的折射率 [n1, n2, n3]
d = [10, 5, inf];        % 层厚度 [d1, d2, d3] (第三层为无限厚)

% 射线参数
theta1_deg = 30;         % 入射角 (度)
source = [0, 0];         % 射线起点坐标 [x, y]

% 仿真参数
num_rays = 5;            % 射线数量（用于角度扫描）
max_reflections = 3;     % 最大反射次数

%% 2. 主射线追踪函数
function ray_paths = trace_ray_three_layers(n, d, theta1_deg, source)
    % 将角度转换为弧度
    theta1 = deg2rad(theta1_deg);
    
    % 初始化射线路径存储
    ray_paths = struct();
    ray_paths.x = [];
    ray_paths.y = [];
    ray_paths.layer = [];
    ray_paths.type = {}; % 'transmit' or 'reflect'
    
    % 当前射线状态
    current_point = source;
    current_theta = theta1;
    current_layer = 1;
    direction = 1; % 1: 向下, -1: 向上
    
    % 追踪射线直到离开系统或达到最大交点
    max_points = 20;
    point_count = 0;
    
    while point_count < max_points
        point_count = point_count + 1;
        
        % 存储当前点
        ray_paths.x(end+1) = current_point(1);
        ray_paths.y(end+1) = current_point(2);
        ray_paths.layer(end+1) = current_layer;
        
        % 计算与当前层边界的交点
        if direction == 1 % 向下传播
            if current_layer == 1
                boundary_y = d(1);
            elseif current_layer == 2
                boundary_y = d(1) + d(2);
            else % 第三层
                % 无限介质，射线继续向下
                next_x = current_point(1) + 100 * tan(current_theta);
                next_y = current_point(2) + 100;
                ray_paths.x(end+1) = next_x;
                ray_paths.y(end+1) = next_y;
                ray_paths.layer(end+1) = current_layer;
                ray_paths.type{end+1} = 'transmit';
                break;
            end
            
            % 计算交点
            dy = boundary_y - current_point(2);
            dx = dy * tan(current_theta);
            next_point = [current_point(1) + dx, boundary_y];
            
        else % 向上传播
            if current_layer == 2
                boundary_y = d(1);
            elseif current_layer == 3
                boundary_y = d(1) + d(2);
            else % 第一层顶部
                boundary_y = 0;
            end
            
            % 计算交点
            dy = boundary_y - current_point(2);
            dx = dy * tan(current_theta);
            next_point = [current_point(1) + dx, boundary_y];
        end
        
        % 存储交点
        ray_paths.x(end+1) = next_point(1);
        ray_paths.y(end+1) = next_point(2);
        ray_paths.layer(end+1) = current_layer;
        
        % 确定边界处的行为（折射或反射）
        if direction == 1 % 向下传播到边界
            if current_layer < 3
                % 计算临界角
                critical_angle = asin(n(current_layer+1)/n(current_layer));
                if abs(current_theta) >= critical_angle
                    % 全反射
                    ray_paths.type{end+1} = 'reflect';
                    direction = -1;
                    current_theta = -current_theta; % 反射角等于入射角
                else
                    % 折射到下一层
                    ray_paths.type{end+1} = 'transmit';
                    current_layer = current_layer + 1;
                    % 应用斯涅尔定律
                    current_theta = asin(n(current_layer-1)/n(current_layer) * sin(current_theta));
                end
            end
        else % 向上传播到边界
            if current_layer > 1
                % 计算临界角
                critical_angle = asin(n(current_layer-1)/n(current_layer));
                if abs(current_theta) <= critical_angle
                    % 折射到上一层
                    ray_paths.type{end+1} = 'transmit';
                    current_layer = current_layer - 1;
                    % 应用斯涅尔定律
                    current_theta = asin(n(current_layer+1)/n(current_layer) * sin(current_theta));
                else
                    % 全反射
                    ray_paths.type{end+1} = 'reflect';
                    direction = 1;
                    current_theta = -current_theta;
                end
            else % 到达顶部边界
                ray_paths.type{end+1} = 'escape';
                break;
            end
        end
        
        % 更新当前点
        current_point = next_point;
    end
end

%% 3. 执行射线追踪并可视化
figure('Position', [100, 100, 1200, 800]);

% 绘制介质层
subplot(2, 3, [1, 2, 4, 5]);
hold on;

% 绘制介质层边界
y_boundaries = [0, d(1), d(1)+d(2), d(1)+d(2)+10];
for i = 1:length(y_boundaries)-1
    if i <= length(n)
        color = [0.8, 0.9, 1.0]; % 浅蓝色表示介质
        if i == 1
            color = [1.0, 0.95, 0.9]; % 第一层用浅橙色
        elseif i == 2
            color = [0.9, 1.0, 0.95]; % 第二层用浅绿色
        end
        
        % 填充介质区域
        fill([-20, 40, 40, -20], ...
             [y_boundaries(i), y_boundaries(i), y_boundaries(i+1), y_boundaries(i+1)], ...
             color, 'FaceAlpha', 0.3, 'EdgeColor', 'k');
    end
end

% 标注介质层
text(-18, d(1)/2, sprintf('层1: n=%.1f', n(1)), 'FontSize', 10, 'FontWeight', 'bold');
text(-18, d(1) + d(2)/2, sprintf('层2: n=%.1f', n(2)), 'FontSize', 10, 'FontWeight', 'bold');
text(-18, d(1) + d(2) + 5, sprintf('层3: n=%.1f', n(3)), 'FontSize', 10, 'FontWeight', 'bold');

% 追踪多条射线（不同入射角）
angles = linspace(10, 60, num_rays);
colors = jet(num_rays);

for i = 1:num_rays
    ray = trace_ray_three_layers(n, d, angles(i), source);
    
    % 绘制射线路径
    plot(ray.x, ray.y, 'Color', colors(i, :), 'LineWidth', 1.5, ...
        'DisplayName', sprintf('θ=%.1f°', angles(i)));
    
    % 标记转折点
    for j = 1:length(ray.type)
        idx = j*2; % 转折点的索引
        if ~isempty(ray.type{j})
            if strcmp(ray.type{j}, 'reflect')
                plot(ray.x(idx), ray.y(idx), 'ro', 'MarkerSize', 6, 'MarkerFaceColor', 'r');
            elseif strcmp(ray.type{j}, 'transmit')
                plot(ray.x(idx), ray.y(idx), 'go', 'MarkerSize', 6, 'MarkerFaceColor', 'g');
            end
        end
    end
end

% 标注源点
plot(source(1), source(2), 'kp', 'MarkerSize', 12, 'MarkerFaceColor', 'y', ...
    'DisplayName', '射线源');

xlabel('水平距离 (x)');
ylabel('深度 (y)');
title('三层介质射线追踪');
legend('Location', 'best');
grid on;
axis equal;
xlim([-20, 40]);
ylim([-2, d(1)+d(2)+10]);

%% 4. 计算传播参数
fprintf('=== 三层介质射线传播参数 ===\n');
fprintf('介质折射率: n1=%.2f, n2=%.2f, n3=%.2f\n', n(1), n(2), n(3));
fprintf('层厚度: d1=%.1f, d2=%.1f\n', d(1), d(2));
fprintf('\n');

% 计算临界角
critical_12 = rad2deg(asin(n(2)/n(1)));
critical_23 = rad2deg(asin(n(3)/n(2)));
critical_13 = rad2deg(asin(n(3)/n(1)));

fprintf('临界角:\n');
fprintf('  层1->层2: %.2f°\n', critical_12);
fprintf('  层2->层3: %.2f°\n', critical_23);
fprintf('  层1->层3: %.2f°\n', critical_13);

%% 5. 分析不同入射角的行为
subplot(2, 3, 3);
theta_test = linspace(0, 90, 100);
transmission_flag = zeros(size(theta_test));

for i = 1:length(theta_test)
    ray = trace_ray_three_layers(n, d, theta_test(i), source);
    % 检查射线是否到达第三层
    transmission_flag(i) = any(ray.layer == 3);
end

plot(theta_test, transmission_flag, 'b-', 'LineWidth', 2);
xlabel('入射角 (度)');
ylabel('到达第三层 (1=是, 0=否)');
title('透射到第三层的角度范围');
grid on;
ylim([-0.1, 1.1]);

%% 6. 计算传播时间
subplot(2, 3, 6);
% 假设光速为1（归一化）
c = 1;
v = c ./ n; % 各层中的传播速度

% 计算示例射线的传播时间
example_angle = 30;
ray_example = trace_ray_three_layers(n, d, example_angle, source);

% 计算路径长度和传播时间
total_time = 0;
segment_times = [];

for i = 1:length(ray_example.x)-1
    dx = ray_example.x(i+1) - ray_example.x(i);
    dy = ray_example.y(i+1) - ray_example.y(i);
    segment_length = sqrt(dx^2 + dy^2);
    
    % 确定当前段的层
    current_layer = ray_example.layer(i);
    segment_time = segment_length / v(current_layer);
    
    total_time = total_time + segment_time;
    segment_times(end+1) = segment_time;
end

bar(1:length(segment_times), segment_times);
xlabel('路径段');
ylabel('传播时间 (归一化)');
title(sprintf('传播时间分布 (θ=%.1f°)', example_angle));
grid on;

fprintf('\n示例射线 (θ=%.1f°) 传播分析:\n', example_angle);
fprintf('  总路径长度: %.2f\n', sum(sqrt(diff(ray_example.x).^2 + diff(ray_example.y).^2)));
fprintf('  总传播时间: %.4f (归一化)\n', total_time);
fprintf('  经过的层: %s\n', mat2str(unique(ray_example.layer)));

sgtitle('三层介质射线追踪综合分析', 'FontSize', 14, 'FontWeight', 'bold');

高级版本：包含能量衰减和参数分析

%% 三层介质射线追踪 - 高级版本（含能量计算）
clear; close all; clc;

%% 1. 高级参数设置
% 介质参数（折射率 + 衰减系数）
n = [1.0, 1.5, 2.0];          % 折射率
alpha = [0.01, 0.05, 0.02];   % 衰减系数 (dB/m)
d = [10, 5, 20];              % 层厚度

% 射线参数
theta_range = [10, 70];       % 入射角范围
num_rays = 50;                % 射线数量

% 源参数
source_power = 1.0;           % 初始功率 (归一化)

%% 2. 高级射线追踪函数（含能量计算）
function [ray, energy] = trace_ray_advanced(n, alpha, d, theta_deg, source, initial_power)
    % 初始化
    theta = deg2rad(theta_deg);
    current_point = source;
    current_theta = theta;
    current_layer = 1;
    direction = 1;
    
    % 存储结构
    ray.x = source(1);
    ray.y = source(2);
    ray.layer = 1;
    ray.event = {};
    ray.theta = theta;
    
    energy.power = initial_power;
    energy.current = initial_power;
    energy.attenuation = 0;
    
    max_segments = 50;
    
    for seg = 1:max_segments
        % 确定当前边界
        if direction == 1 % 向下
            if current_layer == 1
                boundary_y = d(1);
            elseif current_layer == 2
                boundary_y = d(1) + d(2);
            else
                boundary_y = d(1) + d(2) + d(3);
            end
        else % 向上
            if current_layer == 2
                boundary_y = d(1);
            elseif current_layer == 3
                boundary_y = d(1) + d(2);
            else
                boundary_y = 0;
                break; % 到达顶部
            end
        end
        
        % 计算交点
        dy = boundary_y - current_point(2);
        dx = dy * tan(current_theta);
        next_point = [current_point(1) + dx, boundary_y];
        
        % 计算路径衰减
        segment_length = sqrt(dx^2 + dy^2);
        attenuation = alpha(current_layer) * segment_length;
        energy.current = energy.current * exp(-attenuation);
        energy.attenuation = energy.attenuation + attenuation;
        
        % 存储路径点
        ray.x(end+1) = next_point(1);
        ray.y(end+1) = next_point(2);
        ray.layer(end+1) = current_layer;
        ray.theta(end+1) = current_theta;
        
        % 边界处理
        if direction == 1 && current_layer < 3
            % 计算反射/透射系数（简化Fresnel公式）
            n1 = n(current_layer);
            n2 = n(current_layer+1);
            
            if abs(current_theta) >= asin(n2/n1)
                % 全反射
                ray.event{end} = 'total_reflection';
                direction = -1;
                current_theta = -current_theta;
            else
                % 计算透射系数
                theta_t = asin(n1/n2 * sin(current_theta));
                R = ((n1*cos(current_theta) - n2*cos(theta_t)) / ...
                     (n1*cos(current_theta) + n2*cos(theta_t)))^2;
                T = 1 - R;
                
                % 能量分配
                energy.reflected = energy.current * R;
                energy.current = energy.current * T;
                
                % 进入下一层
                ray.event{end} = 'transmission';
                current_layer = current_layer + 1;
                current_theta = theta_t;
            end
            
        elseif direction == -1 && current_layer > 1
            % 类似处理向上传播
            n1 = n(current_layer);
            n2 = n(current_layer-1);
            
            if abs(current_theta) <= asin(n2/n1)
                theta_t = asin(n1/n2 * sin(current_theta));
                R = ((n1*cos(current_theta) - n2*cos(theta_t)) / ...
                     (n1*cos(current_theta) + n2*cos(theta_t)))^2;
                
                energy.current = energy.current * (1 - R);
                ray.event{end} = 'transmission_up';
                current_layer = current_layer - 1;
                current_theta = theta_t;
                direction = 1;
            else
                ray.event{end} = 'reflection';
                direction = 1;
                current_theta = -current_theta;
            end
        else
            break;
        end
        
        current_point = next_point;
        
        % 检查能量是否过低
        if energy.current < 0.001
            ray.event{end} = 'energy_depleted';
            break;
        end
        
        % 检查是否离开系统
        if current_point(2) > sum(d(1:3)) + 5
            ray.event{end} = 'exit_bottom';
            break;
        end
    end
    
    energy.final = energy.current;
end

%% 3. 批量射线追踪与参数分析
fprintf('=== 高级三层介质射线追踪分析 ===\n');
fprintf('追踪 %d 条射线，入射角范围 %.1f° 到 %.1f°\n', ...
    num_rays, theta_range(1), theta_range(2));

angles = linspace(theta_range(1), theta_range(2), num_rays);
results = struct();

for i = 1:num_rays
    [ray, energy] = trace_ray_advanced(n, alpha, d, angles(i), [0, 0], source_power);
    
    results(i).angle = angles(i);
    results(i).ray = ray;
    results(i).energy = energy;
    results(i).final_depth = ray.y(end);
    results(i).final_power = energy.final;
    results(i).num_segments = length(ray.x) - 1;
end

%% 4. 综合可视化
figure('Position', [100, 100, 1400, 900]);

% 4.1 射线路径图
subplot(3, 4, [1, 2, 5, 6]);
hold on;

% 绘制介质层
y_levels = [0, cumsum(d)];
for i = 1:length(y_levels)-1
    color = [0.9, 0.95, 1.0];
    if mod(i, 2) == 0, color = [0.95, 0.9, 0.95]; end
    patch([-30, 60, 60, -30], ...
          [y_levels(i), y_levels(i), y_levels(i+1), y_levels(i+1)], ...
          color, 'FaceAlpha', 0.2, 'EdgeColor', 'k');
end

% 选择几条代表性射线绘制
selected_indices = round(linspace(1, num_rays, 8));
colors = parula(length(selected_indices));

for idx = 1:length(selected_indices)
    i = selected_indices(idx);
    ray = results(i).ray;
    
    % 根据最终能量设置线宽
    linewidth = 1 + 2 * results(i).final_power;
    
    plot(ray.x, ray.y, 'Color', colors(idx, :), ...
        'LineWidth', linewidth, ...
        'DisplayName', sprintf('θ=%.1f°', results(i).angle));
end

xlabel('水平距离 (m)');
ylabel('深度 (m)');
title('代表性射线路径');
legend('Location', 'best', 'NumColumns', 2);
grid on;
xlim([-10, 50]);
ylim([-2, sum(d)+5]);

% 4.2 能量衰减随角度变化
subplot(3, 4, 3);
final_powers = [results.final_power];
plot(angles, 10*log10(final_powers), 'b-', 'LineWidth', 2);
xlabel('入射角 (度)');
ylabel('最终功率 (dB)');
title('最终接收功率 vs 入射角');
grid on;

% 4.3 穿透深度分析
subplot(3, 4, 4);
final_depths = [results.final_depth];
plot(angles, final_depths, 'r-', 'LineWidth', 2);
xlabel('入射角 (度)');
ylabel('最终深度 (m)');
title('射线穿透深度');
grid on;
ylim([0, sum(d)+2]);

% 4.4 传播路径长度
subplot(3, 4, 7);
path_lengths = arrayfun(@(r) sum(sqrt(diff(r.ray.x).^2 + diff(r.ray.y).^2)), results);
plot(angles, path_lengths, 'g-', 'LineWidth', 2);
xlabel('入射角 (度)');
ylabel('总路径长度 (m)');
title('射线总路径长度');
grid on;

% 4.5 能量衰减分布
subplot(3, 4, 8);
attenuations = arrayfun(@(r) r.energy.attenuation, results);
plot(angles, attenuations, 'm-', 'LineWidth', 2);
xlabel('入射角 (度)');
ylabel('总衰减 (dB)');
title('能量总衰减');
grid on;

% 4.6 参数敏感性分析
subplot(3, 4, [9, 10, 11, 12]);

% 测试不同折射率对比
n_variations = {[1.0, 1.3, 1.6], [1.0, 1.5, 2.0], [1.0, 2.0, 3.0]};
styles = {'-', '--', ':'};
legend_text = {};

hold on;
for var_idx = 1:length(n_variations)
    n_test = n_variations{var_idx};
    powers_test = zeros(size(angles));
    
    for i = 1:length(angles)
        [~, energy] = trace_ray_advanced(n_test, alpha, d, angles(i), [0, 0], source_power);
        powers_test(i) = energy.final;
    end
    
    plot(angles, 10*log10(powers_test), 'LineWidth', 2, ...
        'LineStyle', styles{mod(var_idx-1, length(styles))+1}, ...
        'Color', [0, 0, var_idx/length(n_variations)]);
    
    legend_text{end+1} = sprintf('n=[%.1f,%.1f,%.1f]', n_test(1), n_test(2), n_test(3));
end

xlabel('入射角 (度)');
ylabel('接收功率 (dB)');
title('不同折射率配置下的性能比较');
legend(legend_text, 'Location', 'best');
grid on;

sgtitle('三层介质射线追踪高级分析', 'FontSize', 16, 'FontWeight', 'bold');

%% 5. 统计结果输出
fprintf('\n=== 统计结果 ===\n');

% 找到最佳入射角
[max_power, max_idx] = max([results.final_power]);
fprintf('最佳入射角: %.1f° (接收功率: %.4f, %.1f dB)\n', ...
    results(max_idx).angle, max_power, 10*log10(max_power));

% 计算平均性能
avg_power = mean([results.final_power]);
avg_depth = mean([results.final_depth]);
avg_length = mean(path_lengths);

fprintf('平均接收功率: %.4f (%.1f dB)\n', avg_power, 10*log10(avg_power));
fprintf('平均穿透深度: %.2f m\n', avg_depth);
fprintf('平均路径长度: %.2f m\n', avg_length);

% 计算到达第三层的射线比例
reached_layer3 = sum([results.final_depth] > (d(1) + d(2) - 0.1));
fprintf('到达第三层的射线比例: %.1f%%\n', reached_layer3/num_rays*100);

%% 6. 导出数据（可选）
% 将结果保存到文件以便进一步分析
output_data.angles = angles;
output_data.final_powers = [results.final_power];
output_data.final_depths = [results.final_depth];
output_data.path_lengths = path_lengths;
output_data.attenuations = attenuations;

save('ray_tracing_results.mat', 'output_data', 'n', 'alpha', 'd');
fprintf('\n结果已保存到 ray_tracing_results.mat\n');

参考代码 三层介质射线追踪MATLAB代码 www.3dddown.com/cna/96611.html

关键物理原理

1. 斯涅尔定律（折射定律）

n₁·sin(θ₁) = n₂·sin(θ₂)

• 用于计算射线穿过介质界面时的折射角

2. 临界角与全反射

θ_critical = arcsin(n₂/n₁)   （当n₁ > n₂时）

• 入射角大于临界角时发生全反射

3. 菲涅尔方程（反射/透射系数）

R = [(n₁·cosθ₁ - n₂·cosθ₂)/(n₁·cosθ₁ + n₂·cosθ₂)]²
T = 1 - R

• 计算边界处的能量分配

应用场景建议

1. 地震勘探：模拟地震波在不同岩层中的传播
2. 水下声学：研究声波在海水、沉积物中的传播路径
3. 光学薄膜：分析光线在多层光学涂层中的行为
4. 医学成像：超声波在人体不同组织中的传播模拟