【光学】弹性FDTD二维波传播【含Matlab源码 3704期】

Matlab领域

已于 2023-12-30 15:30:52 修改

阅读量856

点赞数 17

分类专栏： Matlab图像处理（初级版）文章标签： matlab

于 2023-12-20 22:25:26 首次发布

本文链接：https://blog.csdn.net/weixin_63266434/article/details/135119092

版权

Matlab图像处理（初级版）专栏收录该内容

369 篇文章 8 订阅

订阅专栏

在这里插入图片描述

⛄一、获取代码方式

获取代码方式1：
完整代码已上传我的资源：【光学】基于matlab弹性FDTD二维波传播【含Matlab源码 3704期】
点击上面蓝色字体，直接付费下载，即可。

获取代码方式2：
付费专栏Matlab物理应用（初级版）

备注：
点击上面蓝色字体付费专栏Matlab物理应用（初级版），扫描上面二维码，付费29.9元订阅海神之光博客付费专栏Matlab物理应用（初级版），凭支付凭证，私信博主，可免费获得1份本博客上传CSDN资源代码（有效期为订阅日起，三天内有效）；
点击CSDN资源下载链接：1份本博客上传CSDN资源代码

⛄二、部分源代码

close all;
clear all;

% Output periodicity in time steps
IT_DISPLAY = 10;

%% MODEL
% Model dimensions, [m]
nx = 401;
nz = 401;
dx = 10;
dz = dx;

% Elastic parameters
vp = 3300.0 * ones(nz, nx); % velocity of compressional waves, [m/s]
vs = vp / 1.732; % velocity of shear waves, [m/s]
rho = 2800.0 * ones(size(vp)); % density, [kg/m3]

% Uncomment below to load Marmousi II benchmark model
% load(‘./data/marmousi_ii_resampled.mat’)

% Lame parameters
lam = rho.(vp.^2 - 2vs.^2); % first Lame parameter
mu = rho.*vs.^2; % shear modulus, [N/m2]

%% TIME STEPPING
t_total = .7; % [sec] recording duration
dt = 1/(max(vp(😃) * sqrt(1.0/dx^2 + 1.0/dz^2));
nt = round(t_total/dt); % number of time steps
t = [0:nt]*dt;

%% SOURCE
f0 = 10.0; % dominant frequency of the wavelet
t0 = 1.20 / f0; % excitation time
factor = 1e10; % amplitude coefficient
angle_force = 90.0; % spatial orientation

jsrc = round(nz/2); % source location along OZ
isrc = round(nx/2); % source location along OX

deg2rad = pi / 180.d0; % convet degrees to radians
a = pipif0f0;
dt2rho_src = dt^2/rho(jsrc, isrc);
source_term = factor * exp(-a(t-t0).^2); % Gaussian
% source_term = -factor2.0a*(t-t0)exp(-a(t-t0)^2); % First derivative of a Gaussian:
% source_term = -factor * (1.0 - 2.0a(t-t0).^{2).*exp(-a*(t-t0).}2); % Ricker source time function (second derivative of a Gaussian):

force_x = sin(angle_force * deg2rad) * source_term * dt2rho_src / (dx * dz);
force_z = cos(angle_force * deg2rad) * source_term * dt2rho_src / (dx * dz);
min_wavelengh = 0.5*min(vs(vs>330))/f0; % shortest wavelength bounded by velocity in the air

%% ABSORBING BOUNDARY (ABS)
abs_thick = min(floor(0.20nx), floor(0.20nz)); % thicknes of the layer
abs_rate = 0.3/abs_thick; % decay rate

lmargin = [abs_thick abs_thick];
rmargin = lmargin;
weights = ones(nz+2,nx+2);
for iz = 1:nz+2
for ix = 1:nx+2
i = 0;
j = 0;
k = 0;
if (ix < lmargin(1) + 1)
i = lmargin(1) + 1 - ix;
end
if (iz < lmargin(2) + 1)
k = lmargin(2) + 1 - iz;
end
if (nx - rmargin(1) < ix)
i = ix - nx + rmargin(1);
end
if (nz - rmargin(2) < iz)
k = iz - nz + rmargin(2);
end
if (i == 0 && j == 0 && k == 0)
continue
end
rr = abs_rate * abs_rate * double(ii + jj + k*k );
weights(iz,ix) = exp(-rr);
end
end