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);