Dear k-WAVE,

Currently, I am trying to see the nonlinear acoustic field with a chirp input. I have some questions regarding the simulation, even though K-wave could give me some initial results.

1) the power law effect, as the frequency at each instantaneous moment was different(time dependent), how did k-wave consider this alteration ? and how could k-wave consider the frequency change at each moment?

2) the simulation time, as the chirp time used in the simulation is a little bit long (around 1ms) and the working frequency is 3MHz, the current simulation takes a very long, is there any ways to speed up

please find my code below

clear all;

% =========================================================================

% SIMULATION

% =========================================================================

PML_size=20;

DATA_CAST = 'single'; % set to 'single' or 'gpuArray-single' to speed up computations

% create the computational grid

Nx = 1024-2*PML_size; % number of grid points in the x (row) direction

Ny = 1024-2*PML_size; % number of grid points in the y (column) direction

dx = 1e-4; % grid point spacing in the x direction [m]

dy = dx; % grid point spacing in the y direction [m]

kgrid = makeGrid(Nx, dx, Ny, dy);

% define a curved transducer element

R=0.0626; % tranducer focal length

r1=0.032; % Aperture size-outer radius

r2=0.0113; % Aperture size-iner radius

% I=zeros(Nx,Ny);

angle1=asin(r1/R);

angle2=asin(r2/R);

L_shift_1=round(0.01/dx);

L_shift_2=round(0.02/dx);

R_size=round(R/dx);

circle1 = makeCircle(Nx, Ny,Nx/2,Ny/2+L_shift_2,R_size,angle1);

circle2 = makeCircle(Nx, Ny,Nx/2,Ny/2+L_shift_2,R_size,angle2);

circle3 = flipud(circle1);

circle4 = flipud(circle2);

circle5=circle3-circle4;

source.p_mask=circle1-circle2+circle5;

% define the properties of the propagation medium

% water

medium.sound_speed = 1500*ones(Nx,Ny); % [m/s]

medium.alpha_power = 2; % [dB/(MHz^y cm)]

medium.alpha_coeff = 0.0022*ones(Nx,Ny); % [dB/(MHz^y cm)]

medium.BonA=5*ones(Nx,Ny);

medium.density=1000*ones(Nx,Ny);

% phantom property

medium.sound_speed(:,Ny/2+L_shift_2:Ny) = 1544; % [m/s]

% medium.alpha_power(:,412:1024) = 1.1; % [dB/(MHz^y cm)]

medium.alpha_coeff(:,Ny/2+L_shift_2:Ny) = 0.015; % [dB/(MHz^y cm)]

medium.density(:,Ny/2+L_shift_2:Ny)=1044;

% create the time array

dt=4e-8;

tend=1e-4;

kgrid.t_array=0:dt:tend;

% define a time varying sinusoidal source

source_freq = 1.0e6; % [Hz]

source_mag = 600E3; % [Pa]

delta_freq=0.4e6;

source.p = source_mag*sin(2*pi*(source_freq+delta_freq*kgrid.t_array/tend).*kgrid.t_array);

% filter the source to remove high frequencies not supported by the grid

source.p = filterTimeSeries(kgrid, medium, source.p);

% create a display mask to display the transducer

display_mask = source.p_mask;

% create a sensor mask covering the entire computational domain using the

% opposing corners of a rectangle

sensor.mask = [1, 1, Nx, Ny].';

% set the record mode capture the final wave-field and the statistics at

% each sensor point

sensor.record = {'p_final', 'p_min', 'p_rms'};

% assign the input options

input_args = {'DisplayMask', display_mask, 'PMLInside', false, 'PlotPML', false,'DataCast', DATA_CAST, 'DataRecast', true};

% run the simulation

sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% =========================================================================

% VISUALISATION

% =========================================================================

% add the source mask onto the recorded wave-field

sensor_data.p_final(source.p_mask ~= 0) = 1;

sensor_data.p_min(source.p_mask ~= 0) = 1;

sensor_data.p_rms(source.p_mask ~= 0) = 1;

% plot the final wave-field

figure;

subplot(1, 3, 1), imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data.p_final, [-1e6 1e6]);

colormap(getColorMap);

ylabel('x-position [mm]');

xlabel('y-position [mm]');

axis image;

title('Final Wave Field');

% plot the maximum recorded pressure

subplot(1, 3, 2), imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data.p_min, [-1e6 1e6]);

colormap(getColorMap);

ylabel('x-position [mm]');

xlabel('y-position [mm]');

axis image;

title('Minimum Pressure');

% plot the rms recorded pressure

subplot(1, 3, 3), imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data.p_rms, [-1e6 1e6]);

colormap(getColorMap);

ylabel('x-position [mm]');

xlabel('y-position [mm]');

axis image;

title('RMS Pressure');

scaleFig(2, 1);