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);