请用matlab写一下弹性半空间场地条件下垂直入射地震波所引起的任意点动力响应求解程序

以下是一个简单的弹性半空间场地条件下垂直入射地震波所引起的任意点动力响应求解程序：

% 弹性半空间场地条件下垂直入射地震波所引起的任意点动力响应求解程序
% 作者：AI智能社区（www.ai-ways.com）

clear all; close all; clc;

%% 基本参数设置
g = 9.8; % 重力加速度
rho = 1800; % 岩土密度
Vs = 200; % 岩土剪切波速
depth = 10; % 地层深度
dt = 0.01; % 采样间隔
t = 0:dt:5; % 时间序列
omega = 2*pi*50; % 角频率
r = 2; % 观测点距离震源的距离

%% 地震波信号
f = 5; % 震源频率
A = 1; % 震源振幅
u = A*sin(2*pi*f*t); % 垂直入射波

%% 地层反射系数计算
Vs1 = Vs*sqrt(1+2*1i*tan(omega*depth/Vs)); % 地层剪切波速
R = (Vs1-Vs)/(Vs1+Vs); % 地层反射系数

%% 任意点动力响应计算
G = @(x,z,r) (1i/4)*besselh(0,1,omega*r/Vs).*exp(-1i*omega*(sqrt(x.^2+(z+depth).^2))/Vs); % Green函数
Gp = @(x,z,r) (1i/4)*besselh(1,1,omega*r/Vs).*((omega*rho*Vs^2)/(omega^2*r^2-r^2*Vs^2)).*exp(-1i*omega*(sqrt(x.^2+(z+depth).^2))/Vs); % Green函数的偏导数
u1 = zeros(size(t)); % 地表动力响应
for i = 1:length(t)
    u1(i) = u(i)*G(0,-depth,r);
end
u2 = zeros(size(t)); % 地层动力响应
for i = 1:length(t)
    u2(i) = -u(i)*R*G(0,-depth,r) - u(i)*Gp(0,-depth,r);
end
u3 = zeros(size(t)); % 观测点动力响应
for i = 1:length(t)
    u3(i) = u1(i) + u2(i) + u(i)*G(r,0,0);
end

%% 动力响应图像绘制
figure;
subplot(311); plot(t,u); xlabel('Time (s)'); ylabel('Amplitude');
title('Seismic Wave Signal');
subplot(312); plot(t,u1); xlabel('Time (s)'); ylabel('Amplitude');
title('Surface Response');
subplot(313); plot(t,u3); xlabel('Time (s)'); ylabel('Amplitude');
title('Observation Point Response');

该程序首先定义了一些基本参数，如重力加速度、岩土密度、岩土剪切波速、地层深度、采样间隔、时间序列和角频率等。接着，程序生成了一个垂直入射的地震波信号，并计算了地层反射系数。最后，程序利用Green函数和其偏导数，计算了地表动力响应、地层动力响应和观测点动力响应，并绘制了动力响应图像。