你能给出二维面波正演的matlab代码吗

以下是一个简单的二维面波正演的 matlab 代码：

% 定义模型参数和波形参数
vp = 2000;      % P波速度
vs = 1000;      % S波速度
rho = 2000;     % 密度
f = 20;         % 主频
dx = 5;         % 网格间距
dt = 0.001;     % 时间间隔
nt = 1000;      % 时间步数
nx = 100;       % 网格数

% 初始化波场和速度模型
u = zeros(nx, nx, nt);   % 波场
v = vp * ones(nx, nx);   % P波速度

% 计算S波速度模型
vs_model = vs * ones(nx, nx);   % S波速度模型

% 计算P波和S波的传播时间
tp = dx / vp;   % P波传播时间
ts = dx / vs;   % S波传播时间

% 计算初始的震源函数
t = (0 : nt-1) * dt;    % 时间轴
src = 2 * pi * f * exp(-pi^2 * f^2 * (t - 1/f).^2);   % 震源函数

% 进行时间步进计算
for it = 2 : nt
    % 计算P波和S波的时间差
    if mod(it, round(ts/dt)) == 0
        v = vs_model;
    end
    
    % 计算波场
    for i = 2 : nx-1
        for j = 2 : nx-1
            u(i,j,it) = 2*u(i,j,it-1) - u(i,j,it-2) + ...
                        v(i,j)^2 * dt^2 * ( ...
                            (u(i+1,j,it-1) - 2*u(i,j,it-1) + u(i-1,j,it-1)) / dx^2 + ...
                            (u(i,j+1,it-1) - 2*u(i,j,it-1) + u(i,j-1,it-1)) / dx^2 ...
                        );
        end
    end
    
    % 加上震源
    u(50,50,it) = u(50,50,it) + src(it);
    
    % 绘制波场图像
    surf(u(:,:,it)), shading interp, axis equal, view(0,90), caxis([-0.1,0.1]), drawnow
end

注：此代码仅用于学习参考，实际应用中需要进行更严谨的模型参数设置和波形参数设置。