Hi Brad,
is it ok, just post a script here?
----------------------------------------------------------------------------
clear all; close all;
clearvars;
% simulation settings
DATA_CAST = 'single';
% =========================================================================
% DEFINE THE K-WAVE GRID
% =========================================================================
% set the size of the perfectly matched layer (PML)
PML_X_SIZE = 10; % [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 = 160 * 2 - 2*PML_X_SIZE; % [grid points]
Ny = 190 * 2 - 2*PML_Y_SIZE; % [grid points]
Nz = 190 * 2 - 2*PML_Z_SIZE; % [grid points]
% calculate the spacing between the grid points
dx = 0.2e-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
c0 = 1540; % [m/s]
rho0 = 1000; % [kg/m^3]
medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;
medium.BonA = 6;
% create the time array
t_end = 40e-6;
kgrid.makeTime(c0, [], t_end);
% =========================================================================
% DEFINE THE SOURCE
% =========================================================================
% create a concave sensor
sphere_offset = 0; % [grid points]
diameter = round(64e-3 / dx) +1 ; % [grid points]
radius = round(45e-3 / dx); % [grid points]
bowl_pos = [1 + sphere_offset, Ny/2, Nz/2]; % [grid points]
focus_pos = [Nx/2, Ny/2, Nz/2]; % [grid points]
source.p_mask = makeBowl([Nx, Ny, Nz], bowl_pos, radius, diameter, focus_pos) ...
- makeBowl([Nx, Ny, Nz], bowl_pos, radius, diameter/2, focus_pos);
% define properties of the input signal
source_strength = 1e6; % [Pa]
tone_burst_freq = 1e6; % [Hz]
tone_burst_cycles = 5;
% create the input signal using toneBurst
source.p = source_strength .* toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);
% =========================================================================
% DEFINE THE ULTRASOUND TRANSDUCER
% =========================================================================
% physical properties of the transducer
transducer.number_elements = 128; % total number of transducer elements
transducer.element_width = round(0.2e-3 / dx); % width of each element [grid points/voxels]
transducer.element_length = round(3.5e-3 / dx); % length of each element [grid points/voxels]
transducer.element_spacing = 0; % spacing (kerf width) between the elements [grid points/voxels]
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([80, 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 = 0e-3; % focus distance [m]
transducer.elevation_focus_distance = 0e-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);
% create the transducer using the defined settings
transducer = kWaveTransducer(kgrid, transducer);
% print out transducer properties
transducer.properties;
% voxelPlot(double(transducer.mask)+double(source.p_mask));
% view(127, 18);
% =========================================================================
% DEFINE THE MEDIUM PROPERTIES
% =========================================================================
% 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, Ny, Nz]);
% define a random distribution of scatterers for the highly scattering
% region
scattering_map = randn([Nx, Ny, Nz]);
scattering_c0 = c0 + 25 + 75 * scattering_map;
scattering_c0(scattering_c0 > 1600) = 1600;
scattering_c0(scattering_c0 < 1400) = 1400;
scattering_rho0 = scattering_c0 / 1.5;
% define properties
sound_speed_map = c0 * ones(Nx, Ny, Nz) .* background_map;
density_map = rho0 * ones(Nx, Ny, Nz) .* background_map;
% define a sphere for a highly scattering region
radius = 8e-3;
x_pos = 32e-3;
y_pos = dy * Ny/2;
scattering_region1 = makeBall(Nx, Ny, Nz, round(x_pos / dx), round(y_pos / dx), Nz/2, round(radius / dx));
% assign region
sound_speed_map(scattering_region1 == 1) = scattering_c0(scattering_region1 == 1);
density_map(scattering_region1 == 1) = scattering_rho0(scattering_region1 == 1);
% assign to the medium inputs
medium.sound_speed = sound_speed_map;
medium.density = density_map;
% =========================================================================
% RUN THE SIMULATION
% =========================================================================
% set the input settings
input_args = {'DisplayMask', transducer.active_elements_mask, ...
'PMLInside', false, 'PlotPML', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
'DataCast', DATA_CAST, 'PlotScale', [-1/4, 1/4] * source_strength};
% run the simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, transducer, input_args{:});
% extract a single scan line from the sensor data using the current
% beamforming settings
scan_line = transducer.scan_line(sensor_data);
% =========================================================================
% VISUALISATION
% =========================================================================
% plot the recorded time series
figure;
stackedPlot(kgrid.t_array * 1e6, sensor_data);
xlabel('Time [\mus]');
ylabel('Transducer Element');
title('Recorded Pressure');
scaleFig(1, 1.5);
------------------------------------------------------------------------------
Thanks
Gwansuk