本专栏(一)至(五)文章均介绍的是二维模拟，本文将介绍如何进行一维模拟和三维模拟。
同时也包含：定义网格、介质、压力分布(声源)、传感器；运行模拟；数据可视化几个部分。
有兴趣的同好可以加企鹅群交流：937503xxx（有意向的同好私聊或评论加群，新群，目前群内仅5人），任何科研相关问题都可以在群内交流，平常一起聊聊天、分享日常也是没问题的，欢迎~
目录
一、一维模拟
此模拟和检测是一维传播介质中的初始压力分布产生的压力场的简单演示。
基于此前的均匀传播介质和异构传播介质的示例。
1.1 定义网格、介质、压力分布
% create the computational grid
Nx = 512; % number of grid points in the x (row) direction
dx = 0.05e-3; % grid point spacing in the x direction [m]
kgrid = kWaveGrid(Nx, dx);
% define the properties of the propagation medium
medium.sound_speed = 1500 * ones(Nx, 1); % [m/s]
medium.sound_speed(1:round(Nx/3)) = 2000; % [m/s]
medium.density = 1000 * ones(Nx, 1); % [kg/m^3]
medium.density(round(4*Nx/5):end) = 1500; % [kg/m^3]
% create initial pressure distribution using a smoothly shaped sinusoid
x_pos = 280; % [grid points]
width = 100; % [grid points]
height = 1; % [au]
in = (0:pi/(width/2):2*pi).';
source.p0 = [zeros(x_pos, 1); ((height/2) * sin(in - pi/2) + (height/2));
zeros(Nx - x_pos - width - 1, 1)];（1）定义网格：
% 一维仿真以类似于二维仿真的方式进行。
% 介质离散化由"kWaveGrid"使用单一维度的输入来执行。
% 非均质声传播介质的特性也以一维列向量形式给出。
% 此一维网格为[-12.8mm,12.8mm]，可理解为只有x轴。
（2）定义介质：
% 与异构传播介质示例中方法一致，只是换成了一维。
% 这里的"round"为：四舍五入到最接近的小数或整数。
（3）定义初始压力分布：
% 与二维中一样，初始压力分布使用包含计算域中每个网格点的初始压力值的向量来设置。
% 这里一个平滑变化的压力函数是用正弦波的一部分来定义的。
x_pos = 280; % [grid points]
width = 100; % [grid points]
height = 1; % [au]
in = (0:pi/(width/2):2*pi).';
% 创建一个列向量"in"，其中包含从0到2π的等间隔的值，间隔为π/(width/2)=π/50
% ‘.’运算符表示转置操作，因此这里"in"为一个列向量，而不是行向量。
source.p0 = [zeros(x_pos, 1); ((height/2) * sin(in - pi/2) + (height/2));
zeros(Nx - x_pos - width - 1, 1)];
% 由三角函数诱导公式得：(height/2) * sin(in - pi/2) + (height/2) = -1/2*cos(in)+1/2
% 由上面所得区间与间隔，可知函数"-1/2*cosx+1/2"所取值为0，0.0010，0.0039 ...
% 于是可知——
% [-12.8mm,1.2mm]为0，因为：zeros(x_pos, 1)=zeros(280, 1)，-12.8+280*0.05=1.2[mm]：
% [1.2mm,6.25mm]为"-1/2*cosx+1/2"函数在[0,2π]间的取值；
% [6.25mm,12.8mm]为0，因为：zeros(Nx-x_pos-width-1,1)=zeros(131,1)，12.8-131 *0.05 = 6.25[mm]
% 参见下方初始压力分布图，可以看到是正弦函数的一部分，方便对照理解。
1.2 定义传感器及运行模拟
% create a Cartesian sensor mask
sensor.mask = [-10e-3, 10e-3]; % [mm]
% set the simulation time to capture the reflections
t_end = 2.5 * kgrid.x_size / max(medium.sound_speed(:));
% define the time array
kgrid.makeTime(medium.sound_speed, [], t_end);
% run the simulation
sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor, 'PlotLayout', true);
% 同样，传感器定义了在每个时间步长记录压力场的位置，
% 可以以笛卡尔坐标、二进制或直线两端的网格坐标给出。
% 在这里定义了一个带有两点的笛卡尔式传感器：即在-10mm与10mm处定义传感器。
sensor.mask = [-10e-3, 10e-3]; % [mm]
% 设置模拟时间并利用"makeTime"生成时间数组。
t_end = 2.5 * kgrid.x_size / max(medium.sound_speed(:));
% 单独运行可得：kgrid.x_size = 0.0256
% 单独运行可得： max(medium.sound_speed(:)) = 2000
% 所以：t_end = 3.2*10^(-5) = 0.32[us]
kgrid.makeTime(medium.sound_speed, [], t_end);
sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor, 'PlotLayout', true);
% 注意是"kspaceFirstOrder1D"：1D
% "PlotLayout"：使用"kspaceFirstOrder2D"的内置绘图功能。
% 绘制初始压力分布、传感器、介质属性（声速、密度）
（1）运行动画如下——
（2）内置绘图如下：初始压力分布、传感器、介质(声速、密度)——
1.3 数据可视化
% plot the recorded time signals
figure;
plot(sensor_data(1, :), 'b-');
hold on;
plot(sensor_data(2, :), 'r-');
axis tight;
set(gca, 'YLim', [-0.1, 0.7]);
ylabel('Pressure');
xlabel('Time Step');
legend('Sensor Position 1', 'Sensor Position 2');其中——
hold on; % 继续绘图，两条线在一张图上
axis tight; % 设置坐标轴显示范围为紧凑型
可视化结果如下：随时间压力变化——
二、三维模拟
2.1 定义网格、介质、压力分布
% create the computational grid
Nx = 64; % number of grid points in the x direction
Ny = 64; % number of grid points in the y direction
Nz = 64; % number of grid points in the z direction
dx = 0.1e-3; % grid point spacing in the x direction [m]
dy = 0.1e-3; % grid point spacing in the y direction [m]
dz = 0.1e-3; % grid point spacing in the z direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
% define the properties of the propagation medium
medium.sound_speed = 1500 * ones(Nx, Ny, Nz); % [m/s]
medium.sound_speed(1:Nx/2, :, :) = 1800; % [m/s]
medium.density = 1000 * ones(Nx, Ny, Nz); % [kg/m^3]
medium.density(:, Ny/4:end, :) = 1200; % [kg/m^3]
% create initial pressure distribution using makeBall
ball_magnitude = 10; % [Pa]
ball_x_pos = 38; % [grid points]
ball_y_pos = 32; % [grid points]
ball_z_pos = 32; % [grid points]
ball_radius = 5; % [grid points]
ball_1 = ball_magnitude * makeBall(Nx, Ny, Nz, ball_x_pos, ball_y_pos, ball_z_pos, ball_radius);
ball_magnitude = 10; % [Pa]
ball_x_pos = 20; % [grid points]
ball_y_pos = 20; % [grid points]
ball_z_pos = 20; % [grid points]
ball_radius = 3; % [grid points]
ball_2 = ball_magnitude * makeBall(Nx, Ny, Nz, ball_x_pos, ball_y_pos, ball_z_pos, ball_radius);
source.p0 = ball_1 + ball_2;（1）定义网格：定义方式类似，只不过新增了z轴；
（2）定义介质：定义方式类似，只不过新增了z轴；
（3）定义初始压力分布：
% 使用"makeBall"函数，而不是二维示例中的"makeDisc"。
% 在这里为：创建两个不同直径的球形的初始压力分布。
2.2 定义传感器及运行模拟
% define a series of Cartesian points to collect the data
x = (-22:2:22) * dx; % [m]
y = 22 * dy * ones(size(x)); % [m]
z = (-22:2:22) * dz; % [m]
sensor.mask = [x; y; z];
% input arguments
input_args = {'PlotLayout', true, 'PlotPML', false, ...
'DataCast', 'single', 'CartInterp', 'nearest'};
% run the simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
% plot the simulation layout using voxelplot
voxelPlot(double(source.p0 | cart2grid(kgrid, sensor.mask)));
view([50, 20]);% 在这里创建了一个线性传感器阵列。
x = (-22:2:22) * dx; % [m]
y = 22 * dy * ones(size(x)); % [m]
z = (-22:2:22) * dz; % [m]
sensor.mask = [x; y; z];
input_args = {'PlotLayout', true, 'PlotPML', false, ...
'DataCast', 'single', 'CartInterp', 'nearest'};
% "PlotLayout-true"：使用内置绘图功能；
% "PlotPML-false"：关闭PML(完美匹配层)的显示
% "CartInterp"——
% 当定义笛卡尔式传感器时，用于提取压力的插值模式。
% 如果"CartInterp"设置为"nearest"，将笛卡尔点在内部映射到二进制网格上，
% 如果多个笛卡尔传感器点映射到同一位置，则丢弃重复的点，并且返回的sensor_data将少于sensor指定的点。
% 因此，如果笛卡尔传感器掩模具有密集的点阵列，则返回的时间序列数可能与原始笛卡尔传感器点数不对应。
% 将丢弃的重复点的详细信息打印到命令行。
% 这种情况可以通过增加计算网格的分辨率或减少笛卡尔传感器点的填充来避免。
% "DataCast"——
% "DataCast"设置为"single"，则加快计算时间，但代价是精度损失；
% 也可以设置为“linear”(默认值)来改变插值模式。
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
% 注意是"kspaceFirstOrder3D"：3D
% 并在二维面上绘制三维图的三个平面图(xy zx yz)：
% 初始压力图、介质性质(声速图、密度图)。
% 使用"voxelplot"绘制三维图
voxelPlot(double(source.p0 | cart2grid(kgrid, sensor.mask)));
view([50, 20]);
% "double"：转换为双精度浮点数；
% "cart2grid"：将一组笛卡尔点插值到二进制网格上，传感器的坐标点绘制出来；
% "view(az,el)"：设置三维图的观察视角，单位是度。
（1）运行动画如下——
（2）内置绘图如下：初始压力分布、传感器、介质(声速、密度) + 三维视角图——
2.3 数据可视化
% plot the simulated sensor data
figure;
imagesc(sensor_data, [-1, 1]);
colormap(getColorMap);
ylabel('Sensor Position');
xlabel('Time Step');
colorbar;可视化结果如下：传感器数据图——
2834
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
