考虑P波入射到固-固界面上，在MATLAB中编写程序生成不同入射角下射线路径（如下图所示）。

在MATLAB中模拟P波入射到固-固界面的情况，通常需要使用偏微分方程求解器（如pdepe函数）以及一些图形绘制库（如meshgrid、plot或surf等）。以下是大致步骤：

1. **设置物理模型**:
定义固体介质的属性（比如弹性常数），并确定界面的位置。

2. **网格创建**:
使用`meshgrid`生成二维网格，x轴代表水平方向，y轴代表深度方向，入射角θ作为额外维度。

3. **方程设定**:
P波在各向同性的均匀介质中的波动方程可以表示为平面波边界条件下的拉普拉斯方程。在MATLAB中，这可以用`pdepe`函数来解决。

4. **边界条件**:
- 对于入射面，应用P波的入射条件（例如，垂直入射为`u(x,0,theta) = sin(omega*t)`）。
- 对于界面，可能是反射和透射条件，需要根据物质特性计算相应的边界值。

5. **求解方程**:
调用`pdepe`函数求解偏微分方程，并存储结果在适当的变量中（比如`U`矩阵，其中`U(:,:,i)`对应每个时间步的解）。

6. **绘图显示**:
利用`surf`或`quiver`函数绘制出不同入射角度下P波的传播轨迹。

% 参数设定
L = [0 1]; % 界面区域
d = [0 1]; % 深度范围
tspan = [0 0.1]; % 时间区间
theta = linspace(0, pi/2, num_angles); % 入射角数组

% 创建网格
[X,Y] = meshgrid(L,d);
Z = zeros(size(X)); % 初始化深度数据
Z(:,1) = theta; % 将入射角作为额外维度

% 方程定义
fun = @(~,~,state) spatial_derivative(state, X, Y, Z); % 自定义函数处理空间导数
ic = @(x,y) initial_conditions(x, y, Z(y,1)); % 初始条件函数

% 解方程
[u,t] = pdepe(@stiffness, @flux, ic, tspan, [X(:),Y(:),Z(:)]);

% 可视化结果
surf(X,Y,u)
xlabel('Horizontal Distance')
ylabel('Depth')
zlabel('Amplitude')
title(['P-Wave Propagation for Different Incidence Angles'])
hold on
for i = 1:length(theta)
    quiver(X(:,1), Y(1,:), u(:,1,i), [], 'LineWidth', 1.5);
    line([0 L(1)], [Z(1,i) Z(1,i)], 'Color', 'k'); % 入射线
end
hold off