Bradley, here is my code. The image mask is not included, but essentially it is black and white skull curved strip. The masks apply well (I've checked the matrix).
Not sure why the amplitude isn't attenuating at the focus though, and I'm getting really weird behaviour within the skull (seems to be amplifying everything) instead of attenuating
Thanks for your help
clearvars;
% =========================================================================
% ACOUSTIC SIMULATION
% =========================================================================
% define the PML size
pml_size = 20; % [grid points]
% define the grid parameters
Nx = 256 - 2 * pml_size; % [grid points]
Ny = 256 - 2 * pml_size; % [grid points]
dx = 0.25e-3; % [m]
dy = 0.25e-3; % [m]
% create the computational grid
kgrid = kWaveGrid(Nx, dx, Ny, dy);
%define the properties of the propagation medium
%% ALPHA COEFF
medium.alpha_coeff = 0.05 * ones(Nx, Ny); % [dB/(MHz^y cm)] WATER
%skull
pix=loadimagecustom('skull216.png');
pix = pix .* (8/ max(pix(:))); %set all black pixels to 8, all white to 0.
medium.alpha_coeff =medium.alpha_coeff + pix; % add 8 to the 0.05 water, making the skull layer
medium.alpha_power = 1.43;
%% SPEED
medium.sound_speed = 1450 * ones(Nx, Ny); % WATER default fill! [m/s]
% %skull
pix=loadimagecustom('skull216.png');
pix = pix .* (700/ max(pix(:))); *add 700 to the water speed
medium.sound_speed =medium.sound_speed +pix; % 2500 bone! [kg/m^3]
%% DENSITY
medium.density = 1000 * ones(Nx, Ny); % WATER default fill [kg/m^3]
%skull
pix=loadimagecustom('skull216.png');
pix = pix .* (732 / max(pix(:)));
medium.density =medium.density +pix; % bone! [kg/m^3]
%% define the source parameters
diameter = 0.03; % [m]
radius = 0.03; % [m]
freq = 500000; % [Hz]
amp = 0.53e6; % 0.53[MPa]
% define a focused ultrasound transducer
source.p_mask = makeArc([Nx, Ny], [1, Ny/2], round(radius / dx), round(diameter / dx) + 1, [Nx/2, Ny/2]);
% create a display mask to display the transducer
display_mask = source.p_mask;
% calculate the time step using an integer number of points per period
ppw = 12.08; % points per wavelength
cfl = 0.3; % cfl number
ppp = 41; % points per period
T = 1 / freq; % period [s]
dt = 4.878e-8; % time step [s]
% calculate the number of time steps to reach steady state
t_end = 5.0575e-5;
Nt = round(t_end / dt);
% create the time array
kgrid.setTime(Nt, dt);
% define the input signal
source.p = createCWSignals(kgrid.t_array, freq, amp, 0);
% set the sensor mask to cover the entire grid
sensor.mask = ones(Nx, Ny);
sensor.record = {'p', 'p_max_all'};
% record the last 3 cycles in steady state
num_periods = 3;
T_points = round(num_periods * T / kgrid.dt);
sensor.record_start_index = Nt - T_points + 1;
% set the input arguements
input_args = {'PMLInside', false,'PlotLayout', true, 'PlotPML', false, 'DisplayMask', ...
display_mask, 'PlotScale', [-1, 1] * amp};
% run the acoustic simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
% =========================================================================
% CALCULATE PRESSURE
%===================
% extract the pressure amplitude at each position
p = extractAmpPhase(sensor_data.p, 1/kgrid.dt, freq);
% reshape the data,
p = reshape(p, Nx, Ny);
% =========================================================================
% VISUALISATION
% =========================================================================
% plot the thermal dose and lesion map
figure;
% plot the acoustic pressure
subplot(2, 3, 1);
imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, p * 1e-6);
h = colorbar;
xlabel(h, '[MPa]');
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('Acoustic Pressure Amplitude');
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
% set colormap and enlarge figure window
colormap(jet(256));
scaleFig(1.5, 1);