用matlab写s波、p波的反射率和折射率随入射角的变化的代码

在MATLAB中，你可以使用数值计算库来模拟地震波的传播特性。以下是一个简单的示例，演示如何编写一个函数来计算S波（纵波）和P波（横波）的反射率和折射率随着入射角变化的情况。请注意，这只是一个基本框架，并未包括所有细节，如岩石模型的具体设置和复杂计算。

function [R_s, R_p, n_s, n_p] = wave_propagation(theta)
    % 假设你有一个岩体模型（比如各向同性的弹性常数）
    Young's modulus = E;
    Poisson's ratio = nu;

    % S波速度 (纵波) 和 P波速度 (横波)
    vs = sqrt(E / (2 * (1 + nu)));
    vp = sqrt((E + nu * E) / (2 * (1 - nu)));

    % 入射角
    theta_deg = deg2rad(theta);

    % 反射率公式（假设界面是垂直的）
    R_s = (vp^2 - vs^2) / (vp^2 + vs^2);
    R_p = ((vp^2 - vs^2) * sin^2(theta_deg)) / (vp^2 + vs^2);

    % 折射率公式
    sin_theta_prime = vp ./ vs * sin(theta_deg);
    n_s = vp ./ vs;
    n_p = sin_theta_prime; % 横波折射率等于sin(θ')（斯涅尔定律）

    % 返回结果
    R_s = R_s .* (theta <= 90); % 只考虑第一波反射（小于90度）
    R_p = R_p .* (theta <= 90);
end

% 使用示例：
theta_range = linspace(0, 90, 100); % 角度范围，从0到90度，步长为1度
[R_s_values, R_p_values, n_s_values, n_p_values] = meshgrid(theta_range, [], []);

% 计算并存储各个角度下的反射率和折射率
R_s_matrix = reshape(wave_propagation(theta_range), [], 1);
R_p_matrix = reshape(wave_propagation(theta_range), [], 1);
n_s_matrix = reshape(n_s_values, [], 1);
n_p_matrix = reshape(n_p_values, [], 1);

% 现在你可以画出反射率和折射率随入射角变化的图像
% plot(R_s_matrix, 'b', R_p_matrix, 'r')
% plot(n_s_matrix, 'g', n_p_matrix, 'm')