Hi!
I just started learning KWAVE since last week. I am trying to build a spherical source to model our single element, point focused transducer. I did that. But the problem is when I bump up the peak amplitude of the source sine wave to 2MPa, which is what we expect the transducer to achieve, the simulation is just pure black... like it's been overwhelmed...
How to solve this problem?
I attached my code here:
---------------------------------------
% Simulating Transducer Field Patterns Example
%
% This example demonstrates the use of k-Wave to compute the field pattern
% generated by a curved single element transducer in two dimensions. It
% builds on the Monopole Point Source In A Homogeneous Propagation Medium
% Example.
%
% author: Bradley Treeby
% date: 10th December 2009
% last update: 4th May 2017
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2017 Bradley Treeby
% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.
clearvars;
% =========================================================================
% SIMULATION
% =========================================================================
a=1;
% create the computational grid
Nx = 216/2*a; % number of grid points in the x (row) direction
Ny = 216*a; % number of grid points in the y (column) direction
dx = 50e-3/Nx; % grid point spacing in the x direction [m]
dy = dx; % grid point spacing in the y direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy);
% define the properties of the propagation medium
% medium.sound_speed = 1500; % [m/s]
% medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
% medium.alpha_power = 1.5;
medium.sound_speed = 1500 * ones(Nx, Ny); % [m/s]
medium.sound_speed(:, Ny/4:Ny/4+10) = 1644; % [m/s]
medium.density = 1000 * ones(Nx, Ny); % [kg/m^3]
medium.density(:, Ny/4:Ny/4+10) = 1038; % [kg/m^3]
medium.alpha_coeff = 0.75*ones(Nx,Ny); % [dB/(MHz^y cm)]
medium.alpha_coeff(:,Ny/4:Ny/4+10) = 2.24;
medium.alpha_power = 1.5;
% create the time array
kgrid.makeTime(medium.sound_speed);
% define a curved transducer element
arc_pos = [Nx/2, 1]; % [grid points]
radius = 56*a; % [grid points]
diameter = 111; % [grid points]
focus_pos = [Nx/2, Nx/2]; % [grid points]
source.p_mask = makeArc([Nx, Ny], arc_pos, radius, diameter, focus_pos);
% define a time varying sinusoidal source
source_freq = 0.2e6; % [Hz]
source_mag = 0.5; % [Pa]
source.p = source_mag * sin(2 * pi * source_freq * 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].';
% % define a single sensor point
% sensor.mask = zeros(Nx, Ny);
% sensor.mask(Nx/2, Ny/4) = 1;
% set the record mode capture the final wave-field and the statistics at
% each sensor point
sensor.record = {'p','p_final', 'p_max', 'p_rms'};
% assign the input options
input_args = {'DisplayMask', display_mask, 'PMLInside', false, 'PlotPML', false};
% 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_max(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, [-1 1]);
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_max, [-1 1]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('Maximum Pressure');
% plot the rms recorded pressure
subplot(1, 3, 3);
imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.p_rms, [-1 1]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('RMS Pressure');
scaleFig(2, 1);
[t_sc, scale, prefix] = scaleSI(max(kgrid.t_array(:)));
figure
%source signal amplitude
% subplot(2, 1, 1);
% plot(kgrid.t_array * scale, source.p, 'k-');
% xlabel(['Time [' prefix 's]']);
% ylabel('Signal Amplitude');
% axis tight;
% title('Input Pressure Signal');
%
% subplot(2, 1, 2);
% plot(kgrid.t_array * scale, sensor_data.p, 'r-');
% xlabel(['Time [' prefix 's]']);
% ylabel('Signal Amplitude');
% axis tight;
% title('Sensor Pressure Signal');
%