splitted beam pattern of transducer

clearvars;

% simulation settings
DATA_CAST = 'gpuArray-single'; % set to 'single' or 'gpuArray-single' to speed up computations
MASK_PLANE = 'xz'; % set to 'xy' or 'xz' to generate the beam pattern in different planes
USE_STATISTICS = true; % set to true to compute the rms or peak beam patterns, set to false to compute the harmonic beam patterns

% =========================================================================
% DEFINE THE K-WAVE GRID
% =========================================================================

% set the size of the perfectly matched layer (PML)
PML_X_SIZE = 20; % [grid points]
PML_Y_SIZE = 10; % [grid points]
PML_Z_SIZE = 10; % [grid points]

% set total number of grid points not including the PML
Nx = 540 - 2*PML_X_SIZE; % [grid points]
Ny = 780 - 2*PML_Y_SIZE; % [grid points]
Nz = 380 - 2*PML_Z_SIZE; % [grid points]

% set desired grid size in the x-direction not including the PML
x = 40e-3; % [m]

% calculate the spacing between the grid points
dx = 0.05e-3; % [m]
dy = dx; % [m]
dz = dx; % [m]

% create the k-space grid
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);

% =========================================================================
% DEFINE THE MEDIUM PARAMETERS
% =========================================================================

% define the properties of the propagation medium
medium.sound_speed = 1540; % [m/s]
medium.density = 1000; % [kg/m^3]
medium.alpha_coeff = 0.0;
medium.alpha_power = 0.0;

medium.BonA = 6;

% create the time array
t_end = 20e-6; % [s]
kgrid.makeTime(medium.sound_speed, [], t_end);

% =========================================================================
% DEFINE THE INPUT SIGNAL
% =========================================================================

% define properties of the input signal
source_strength = 1e6; % [Pa]
tone_burst_freq = 7e6; % [Hz]
tone_burst_cycles = 5;

% create the input signal using toneBurst
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);

% scale the source magnitude by the source_strength divided by the
% impedance (the source is assigned to the particle velocity)
input_signal = (source_strength ./ (medium.sound_speed * medium.density)) .* input_signal;

% =========================================================================
% DEFINE THE ULTRASOUND TRANSDUCER
% =========================================================================

% physical properties of the transducer
transducer.number_elements = 1; % total number of transducer elements
transducer.element_width = round(35e-3/dx); % width of each element [grid points]
transducer.element_length = round(15e-3/dx); % length of each element [grid points]
transducer.element_spacing = 0; % spacing (kerf width) between the elements [grid points]
transducer.radius = inf; % radius of curvature of the transducer [m]

% calculate the width of the transducer in grid points
transducer_width = transducer.number_elements * transducer.element_width ...
+ (transducer.number_elements - 1) * transducer.element_spacing;

% use this to position the transducer in the middle of the computational grid
transducer.position = round([1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2]);

% properties used to derive the beamforming delays
transducer.sound_speed = 1540; % sound speed [m/s]
transducer.focus_distance = inf; % focus distance [m]
transducer.elevation_focus_distance = 16.5e-3; % focus distance in the elevation plane [m]
transducer.steering_angle = 0; % steering angle [degrees]

% apodization
transducer.transmit_apodization = 'Rectangular';
transducer.receive_apodization = 'Rectangular';

% define the transducer elements that are currently active
transducer.active_elements = ones(transducer.number_elements, 1);

% append input signal used to drive the transducer
transducer.input_signal = input_signal;

% create the transducer using the defined settings
transducer = kWaveTransducer(kgrid, transducer);

% print out transducer properties
transducer.properties;

%voxelPlot(single(transducer.mask));

% =========================================================================
% DEFINE SENSOR MASK
% =========================================================================

% define a sensor mask through the central plane
sensor.mask = zeros(Nx, Ny, Nz);
switch MASK_PLANE
case 'xy'

% define mask
sensor.mask(:, :, Nz/2) = 1;

% store y axis properties
Nj = Ny;
j_vec = kgrid.y_vec;
j_label = 'y';

case 'xz'

% define mask
sensor.mask(:, Ny/2, :) = 1;

% store z axis properties
Nj = Nz;
j_vec = kgrid.z_vec;
j_label = 'z';

case 'xy-xz'
sensor.mask(:, :, Nz/2) = 1;
sensor.mask(:, Ny/2, :) = 1;

end

% set the record mode such that only the rms and peak values are stored
if USE_STATISTICS
sensor.record = {'p_rms', 'p_max'};
end

% =========================================================================
% RUN THE SIMULATION
% =========================================================================

% set the input settings
input_args = {'DisplayMask', transducer.all_elements_mask, ...
'PMLInside', false, 'PlotPML', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
'DataCast', DATA_CAST, 'DataRecast', true, 'PlotScale', [-1/2, 1/2] * source_strength};

% stream the data to disk in blocks of 100 if storing the complete time
% history
if ~USE_STATISTICS
input_args = [input_args {'StreamToDisk', 100}];
end

% run the simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, sensor, input_args{:});

% =========================================================================
% COMPUTE THE BEAM PATTERN USING SIMULATION STATISTICS
% =========================================================================

if USE_STATISTICS

% reshape the returned rms and max fields to their original position
sensor_data.p_rms = reshape(sensor_data.p_rms, [Nx, Nj]);
sensor_data.p_max = reshape(sensor_data.p_max, [Nx, Nj]);

save("WS_LF_7MHz_xz");

% plot the beam pattern using the pressure maximum
figure;
imagesc(j_vec * 1e3, (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3, sensor_data.p_max * 1e-6);
xlabel([j_label '-position [mm]']);
ylabel('x-position [mm]');
title('Total Beam Pattern Using Maximum Of Recorded Pressure');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Pressure [MPa]');
axis image;

% plot the beam pattern using the pressure rms
figure;
imagesc(j_vec * 1e3, (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3, sensor_data.p_rms * 1e-6);
xlabel([j_label '-position [mm]']);
ylabel('x-position [mm]');
title('Total Beam Pattern Using RMS Of Recorded Pressure');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Pressure [MPa]');
axis image;

% end the example
return

end

k-Wave

A MATLAB toolbox for the time-domain
simulation of acoustic wave fields