Hi, there! First of all, thanks for the new version of k-wave! The example equivalent_source_holography helps me a lot!
I use the equivalent_source_holography example to simulate the nonlinear field of the phased array producing multi-foci field pattern. The rayleigh integral is used to calculate the linear field in the depth of one wavelength （as the value of input_plane_complex).The fundamental wave, 2nd and 3rd harmonic wave can be simulated.
My question is: the fundamental wave in the simulation result is quite different with the linear field of the Rayleigh integral.The main lobe is quite similar, but the side lobes is very different.
I think it is because the grid point spacing is not small enough. However, when I set dx =Lamda/14- Lamda/20 ( dx is he grid point spacing; Lamda is the wavelength of the fundemental wave),the L2 Error does not decline apparently and the harmonic fields become strange(the focus area has nearly no pressure and the boundary area has large pressure instead).
How can the fundamental wave of k-wave coincides with the linear field of the Rayleigh integral?
==========================the matlab script ============================
% Equivalent Source Holography Example
tic
close all
clearvars;
% =========================================================================
% DEFINE SIMULATION SETTINGS
% =========================================================================
%声源属性
Fre=1e6; %[Hz]
Sound_speed=1500; %[m/s]
Lamda=Sound_speed/Fre;%波长
% define grid properties
Nx = 309; % 去除一部分x方向数据，以便快速计算
Ny = 150; % number of grid points in the y direction
Nz = 80; % number of grid points in the z direction
dx = Lamda/10; % grid point spacing [m]
c0 = 1500; % sound speed [m/s]
%阵列参数设置
num=20;%对应曲阵元数量
a = Lamda/2;
b = 4*Lamda;
%生成阵元位置规则
dis= Lamda/2;%规则分布的间距
% define source properties
source_freq = Fre; % [Hz]
rect_x_radius = round((a*num+dis*(num-1))/dx); % [grid points]
rect_y_radius = round(b/dx); % [grid points]
input_plane_z_index = round(Lamda/dx); % [grid points] 大概一个波长
% settings for calculating the equivalent source
grid_expansion = 0;
number_optimisation_steps = 20;
%the data calculated by Rayleigh integral
load('双焦点5mm横切面-深度λ-复数.mat')
input_plane_complex=pr;
figure
imagesc(abs(input_plane_complex))
title('输入声源数据')
axis equal
% create the computational grid
kgrid = kWaveGrid(Nx, dx, Ny, dx, Nz, dx);
cfl = 0.3;
t_end = 1.2*sqrt(kgrid.x_size.^2 + kgrid.y_size.^2 + kgrid.z_size.^2) ./ c0;
% assign medium properties
medium.sound_speed = c0;
medium.density = 1000; % [kg/m^3]
medium.BonA = 5;%5
% compute sampling rates, forcing points-per-period to be an integer
points_per_wavelength = c0 / (source_freq * dx);
points_per_period = round(points_per_wavelength / cfl);
% compute corresponding time spacing
dt = 1 / (points_per_period * source_freq);
% create the time array
Nt = round(t_end / dt);
kgrid.setTime(Nt, dt);
% =========================================================================
% CALCULATE EQUIVALENT SOURCE
% =========================================================================
% offset between the input plane and source plane [grid points]
source_offset = input_plane_z_index - 1;
% calculate equivalent source
[source_estimate, optim_params] = calculateMassSourceCW(input_plane_complex, dx, ...
source_freq, c0, source_offset, grid_expansion, ...
'NumSteps', number_optimisation_steps);
% =========================================================================
% PROJECT MEASURED DATA USING EQUIVALENT SOURCE
% =========================================================================
% assign source mask at the beginning of the grid
clear source;
source.p_mask = zeros(Nx, Ny, Nz);
source.p_mask(:, 1, :) = 1;
% assign the measured data as a dirichlet boundary condition
source.p = createCWSignals(kgrid.t_array, source_freq, abs(source_estimate(:)), angle(source_estimate(:)));
% define a sensor mask through the central plane
sensor.mask = zeros(Nx, Ny, Nz);
% define mask
sensor.mask(:, :, Nz/2) = 1;
% set the start time to only record the last three periods
sensor.record_start_index = kgrid.Nt - 3 * points_per_period + 1;
% assign input arguments
input_args = {...
'PMLInside', false', ...
'PMLSize', 'auto', ..., Nx, Ny
'DataCast', 'single', ...
'DataRecast', true, ...
};
% run k-Wave simulation to project measured data
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
save('sensor_data_0508_nonlinear_Lamda10.mat','sensor_data')
Nt=size(sensor_data,2);
sensor_data = reshape(sensor_data, [Nx, Ny, Nt]);
%f0
[proj_eqs,phase,fre] = extractAmpPhase(sensor_data, 1/kgrid.dt, source_freq, ...
'Dim', 3, ...
'FFTPadding', 1, ...
'Window', 'Rectangular');
x_vec=dx*((1:Nx)-Nx/2);%[m]
y_vec=dx*(1:Ny);
figure
imagesc(1000*x_vec,1000*y_vec,proj_eqs')
xlabel('[mm]')
ylabel('[mm]')
%h2
[proj_eqs,phase,fre] = extractAmpPhase(sensor_data, 1/kgrid.dt,2* source_freq, ...
'Dim', 3, ...
'FFTPadding', 1, ...
'Window', 'Rectangular');
figure
imagesc(1000*x_vec,1000*y_vec,proj_eqs')
xlabel('[mm]')
ylabel('[mm]')
%h3
[proj_eqs,phase,fre] = extractAmpPhase(sensor_data, 1/kgrid.dt,3* source_freq, ...
'Dim', 3, ...
'FFTPadding', 1, ...
'Window', 'Rectangular');
figure
imagesc(1000*x_vec,1000*y_vec,proj_eqs')
xlabel('[mm]')
ylabel('[mm]')