Dear author:
I use the 'kspaceFirstOrder3D' function for 3D simulation, but the results are different when the different simulation step size in z-axis is set. I wonder if some parameters were wrong in the code.
Here is the code. P0 is the input computed through Rayleigh-sommerfeld Integral of focused transducer at 20mm before the focus. The result is different with the dz = 0.1e-3; Nz = 311.
clear;
clc;
close all;
% =========================================================================
% DEFINE THE K-WAVE GRID
% =========================================================================
% initial condition
[FileName, FileDir] = uigetfile('.mat', 'Initial Plane File');
P0 = load(fullfile(FileDir, FileName));
P0 = cell2mat(struct2cell(P0));
% set total number of grid points not including the PML
[Nx, Ny] = size(P0); % [grid points]
Nz = 161; ; % [grid points]
% calculate the spacing between the grid points
dx = 0.2e-3; % [m]
dy = dx; % [m]
dz = 0.2e-3; % [m]
% set desired grid size in the x-direction not including the PML
x = Nx * dx; % [m]
y = Ny * dy;
z = Nz * dz; % [m]
% create the k-space grid
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
% =========================================================================
% DEFINE THE MEDIUM PARAMETERS
% =========================================================================
sound_speed_Water = 1482;
medium.sound_speed = sound_speed_Water * ones(Nx, Ny, Nz); % [m/s]
medium.density = 994 * ones(Nx, Ny, Nz); % [kg/m^3]
medium.alpha_coeff = 0 * ones(Nx, Ny, Nz);
medium.alpha_power = 1;
medium.alpha_mode = 'no_dispersion';
medium.BonA = 5.2 * ones(Nx, Ny, Nz);
pml_size = 10; % PML size
pml_alpha = 2; % PML absorption coefficient [Np/grid point]
% =========================================================================
% DEFINE THE INPUT SIGNAL
% =========================================================================
% create the time array
source_freq = 1.2e6; % [Hz]
kgrid.makeTime(medium.sound_speed);
dT = 1/source_freq;
NumT = ceil(dT/kgrid.dt/10) * 10;
T = kgrid.dt * kgrid.Nt;
kgrid.dt = dT/NumT;
Nt = ceil(T/kgrid.dt);
kgrid.t_array = [0:1:Nt] * kgrid.dt;
% define properties of the input signal
source.p_mask = zeros(Nx, Ny, Nz);
source.p_mask(:, :, 11) = 1;
for i = 1:numel(P0)
source.p(i,:) = abs(P0(i)) * sin(2*pi*source_freq*kgrid.t_array + angle(P0(i)));
end
% define a series of Cartesian points to collect the data
sensor.mask = zeros(Nx, Ny, Nz);
sensor.mask(ceil(Nx/2), :, :) = 1;
sensor.record_start_index = ceil(numel(kgrid.t_array)/2);
% define the field parameters to record
sensor.record = {'p', 'p_max_all','p_final'};
input_args = {'DisplayMask', sensor.mask, 'DataCast', 'single', ...
'CartInterp', 'nearest', 'PMLSize', [10 10 10], 'PMLAlpha', pml_alpha, 'PMLInside', false};
% run the simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});Best wishes
kevin