Hi Brad,
Here is my complete script:
clear all;
filename = 'temp.h5';
alpha_coeff_input = 0.25;
alpha_power = 1.75;
% simulation settings
data_cast = 'single';
run_simulation = true;
% =========================================================================
% 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
mul = 2;
Nx = mul*256 - 2*PML_X_SIZE; % [grid points]
Ny = 256 - 2*PML_Y_SIZE; % [grid points]
Nz = 64 - 2*PML_Z_SIZE; % [grid points]
% set desired grid size in the x-direction not including the PML
x = 45e-3; % [m]
% calculate the spacing between the grid points
dx = x/Nx; % [m]
dy = mul*dx; % [m]
dz = mul*dx; % [m]
% create the k-space grid
kgrid = makeGrid(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]
% create the time array
t_end = (Nx*dx)*2.2/c0; % [s]
kgrid.t_array = makeTime(kgrid, c0, [], t_end);
% =========================================================================
% DEFINE THE INPUT SIGNAL
% =========================================================================
% % define properties of the input signal
source_strength = 1e6; % [Pa]
tone_burst_freq = 3e6; % [Hz]
tone_burst_cycles = 4;
% 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./(c0*rho0)).*input_signal;
% plot(input_signal); return,
% =========================================================================
% DEFINE THE ULTRASOUND TRANSDUCER
% =========================================================================
% physical properties of the transducer
transducer.number_elements = 192; % total number of transducer elements
transducer.element_width = 1; % 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 = 10000e-3; % focus distance [m]
transducer.elevation_focus_distance = 30e-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
number_active_elements = 192;
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;
% =========================================================================
% DEFINE THE MEDIUM PROPERTIES
% =========================================================================
% define a large image size to move across
number_scan_lines = 1;
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 a random distribution of scatterers for the highly scattering
% region
scattering_map = randn([Nx_tot, Ny_tot, Nz_tot]);
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_tot, Ny_tot, Nz_tot).*background_map;
density_map = rho0*ones(Nx_tot, Ny_tot, Nz_tot).*background_map;
% define a sphere for a highly scattering region
radius = 7e-3; % [m]
scattering_region3 = makeBall(Nx_tot, Ny_tot, Nz_tot, 3*Nx_tot/4, Ny_tot/2, Nz_tot/2, round(radius/dy));
scattering_region3 = repmat(scattering_region3(:,:,end/2),[1 1 Nz_tot]);
% assign region
sound_speed_map(scattering_region3 == 1) = scattering_c0(scattering_region3 == 1);
density_map(scattering_region3 == 1) = scattering_rho0(scattering_region3 == 1);
% =========================================================================
% SET HETEREGENEOUS ATTENUATION IN THE INCLUSION
% =========================================================================
alpha_coeff = alpha_coeff_input ;
medium.alpha_power = alpha_power ;
medium.BonA = 5; % coefficient de non-linearit?
% =========================================================================
% RUN THE SIMULATION
% =========================================================================
% set the input settings
input_args = {...
'PMLInside', false, 'PlotSim', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
'DataCast', data_cast};
% run the simulation if set to true, otherwise, load previous results from
% disk
% set medium position
medium_position = 1;
% load the current section of the medium
medium.sound_speed = sound_speed_map(:, medium_position:medium_position + Ny - 1, :);
medium.density = density_map(:, medium_position:medium_position + Ny - 1, :);
medium.alpha_coeff = alpha_coeff;
% run the simulation
kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:}, 'SaveToDisk', filename);
save(filename(1:end-3));
Best
Gaoyang