写一个层状介质二维波动方程数值模拟matlab代码，并解释每行代码的意义

以下是一个基于有限差分法的层状介质二维波动方程数值模拟的Matlab代码：

% 定义常数
c1 = 200; % 第一层波速度
c2 = 400; % 第二层波速度
rho1 = 1; % 第一层密度
rho2 = 2; % 第二层密度

% 计算网格参数
dx = 0.01; % x方向网格间距
dy = 0.01; % y方向网格间距
dt = 0.0005; % 时间步长
x = 0:dx:1; % x方向网格坐标
y = 0:dy:1; % y方向网格坐标
t = 0:dt:1; % 时间坐标

% 初始化网格
nx = length(x);
ny = length(y);
nt = length(t);
u = zeros(nx,ny,nt); % 三维数组表示u(x,y,t)
u(:,:,1) = 0; % 初始条件

% 计算有限差分系数
r1 = c1 * dt / dx;
r2 = c2 * dt / dx;
k1 = (rho1 * c1^2) / (dx^2);
k2 = (rho2 * c2^2) / (dx^2);

% 迭代求解
for k = 2:nt
    for i = 2:nx-1
        for j = 2:ny-1
            % 根据层状介质二维波动方程计算u(i,j,k)
            if y(j) <= 0.5
                u(i,j,k) = 2*u(i,j,k-1) - u(i,j,k-2) + k1*(u(i+1,j,k-1) + u(i-1,j,k-1) + u(i,j+1,k-1) + u(i,j-1,k-1) - 4*u(i,j,k-1));
            else
                u(i,j,k) = 2*u(i,j,k-1) - u(i,j,k-2) + k2*(u(i+1,j,k-1) + u(i-1,j,k-1) + u(i,j+1,k-1) + u(i,j-1,k-1) - 4*u(i,j,k-1));
            end
        end
    end
end

% 绘制波动图像
figure
for k = 1:nt
    surf(x,y,u(:,:,k))
    zlim([-0.5,0.5])
    pause(0.01)
end

代码解释：

第1-4行：定义了介质的两个层的波速度和密度常数。

第6-9行：计算了网格参数，包括x、y、t方向的网格坐标和时间步长。

第12-14行：初始化了三维数组u(x,y,t)，表示波动的振幅。

第16-19行：计算了有限差分系数r和k，用于后面迭代计算波动方程。

第21-28行：使用三重循环迭代求解波动方程，根据层状介质的不同区域，使用不同的波速度和密度常数计算u(x,y,t)。

第31-38行：绘制了波动随时间变化的图像，使用surf()函数。

总体来说，这个代码使用了有限差分法对层状介质二维波动方程进行了数值模拟。