Here is simulation part of my code:
%% DEFINE THE K-WAVE GRID
pml_x_size = 20; % [grid points]
pml_y_size = 10; % [grid points]
pml_z_size = 10; % [grid points]
x_size = 60e-3;
y_size = 30e-3;
z_size = 10e-3;
f_max = 2.7e6;
points_per_wavelength = 3;
c0_min = 1540;
dx = c0_min/(points_per_wavelength*f_max);
dy = dx;
dz = dx;
Nx = round(x_size/dx);
Ny = round(y_size/dy);
Nz = round(z_size/dz);
%% DEFINE THE MEDIUM PARAMETERS
% define the properties of the propagation medium
c0 = 1540; % [m/s]
rho0 = 1000; % [kg/m^3]
medium.alpha_coeff = 0.5*ones(kgrid.Nx, kgrid.Ny, kgrid.Nz); % [dB/(MHz^y cm)]
medium.alpha_coeff(round(0.5*kgrid.Nx):round(0.65*kgrid.Nx),:,:) = 1;
medium.alpha_power = 1.5;
medium.BonA = 6;
% create the time array
t_end = (Nx*dx)*1.9/c0; % [s]
CFL = 0.2;
kgrid.t_array = makeTime(kgrid, c0, CFL, t_end);
%% DEFINE THE INPUT SIGNAL
% define properties of the input signal
source_strength = 1e6; % [Pa]
tone_burst_freq = 3.5e6; % [Hz]
tone_burst_cycles = 4;
signal_burst_freeq = 1/kgrid.dt;
% create the input signal using toneBurst
input_signal = toneBurst(signal_burst_freeq, 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./(c0*rho0)).*input_signal;
%% DEFINE THE ULTRASOUND TRANSDUCER
% physical properties of the transducer
transducer.number_elements = 36; % total number of transducer elements
transducer.element_width = 2; % width of each element [grid points]
transducer.element_length = 24; % 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 = c0; % sound speed [m/s]
transducer.focus_distance = 20e-3; % focus distance [m]
transducer.steering_angle = 0; % steering angle [degrees]
% apodization
transducer.transmit_apodization = 'Hanning';
transducer.receive_apodization = 'Rectangular';
% define the transducer elements that are currently active
%number_active_elements = 32;
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 = makeTransducer(kgrid, transducer);
% print out transducer properties
transducer.properties;
transducer.plot;
%% DEFINE THE MEDIUM PROPERTIES
% define a large image size to move across
number_scan_lines = 9;
Nx_tot = Nx;
Ny_tot = Ny + number_scan_lines*transducer.element_width;
Nz_tot = Nz;
% define a random distribution of scatterers for the medium
background_map_mean = 1;
background_map_std = 0.008;
background_map = background_map_mean + background_map_std*randn([Nx_tot, Ny_tot, Nz_tot]);
% define properties
sound_speed_map = c0*ones(Nx_tot, Ny_tot, Nz_tot).*background_map;
density_map = rho0*ones(Nx_tot, Ny_tot, Nz_tot).*background_map;
%% RUN THE SIMULATION
% preallocate the storage
scan_lines = zeros(number_scan_lines, kgrid.Nt);
% set the input settings
input_args = {...
'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size, pml_z_size], ...
'DataCast', DATA_CAST, 'DataRecast', true, 'PlotSim', false};
if RUN_SIMULATION
% set medium position
% medium_position = 1;
tr_width = transducer.element_width;
for scan_line_index = 1:number_scan_lines
input_args = {'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size, pml_z_size], ...
'DataCast', DATA_CAST, 'DataRecast', true, 'PlotSim', false};
n = Ny;
tr_wid = tr_width;
med = struct();
med = medium;
if (scan_line_index==1)
medium_position = 1
else
medium_position = scan_line_index*tr_width;
end
disp('');
disp(['Computing scan line ' num2str(scan_line_index) ' of ' num2str(number_scan_lines)]);
% load the current section of the medium
med.sound_speed = sound_speed_map(:, medium_position:medium_position + n - 1, :);
med.density = density_map(:, medium_position:medium_position + n - 1, :);
% run the simulation
sensor_data = kspaceFirstOrder3D(kgrid, med, transducer, transducer, input_args{:});
% extract the scan line from the sensor data
scan_lines(scan_line_index, :) = transducer.scan_line(sensor_data);
% update medium position
%medium_position = medium_position + transducer.element_width;
end
% save the scan lines to disk
save example_us_bmode_scan_lines17 scan_lines;
else
% load the scan lines from disk
load example_us_bmode_scan_lines17
end
scan_lines_orig=scan_lines;