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]')