% Script to define the skull geometry used in
% https://dx.doi.org/10.1121/10.0002177 
%
% author: Bradley Treeby
% date: 26th November 2020
% last update: 26th November 2020

clearvars;

% --------------------
% PROPERTIES
% --------------------

% source parameters
source_f0           = 1.1e6;    % source frequency [Hz]
source_roc          = 100e-3;   % bowl radius of curvature [m]
source_diameter     = 64e-3;    % bowl aperture diameter [m]
source_mag          = 2e6;      % source pressure [Pa]
source_cycles       = 4;

% material geometry
skull_radius      	= 80e-3;
skull_thickness  	= 6.5e-3;
skin_thickness      = 6.5e-3;
skin_tx_distance    = 30e-3;

% material properties
water_temp          = 22;
water_c0            = waterSoundSpeed(water_temp);
water_rho0          = waterDensity(water_temp);
water_BonA          = waterNonlinearity(water_temp);
water_a0            = waterAbsorption(source_f0 * 1e-6, water_temp);

skin_c0             = 1624;
skin_rho0           = 1109;
skin_BonA           = 6.7;
skin_a0             = 1.84 / (source_f0 * 1e-6).^2;

brain_c0            = 1546;      
brain_rho0          = 1045;
brain_BonA          = 6.7;
brain_a0            = (0.59 * 1.1^1.3) / (source_f0 * 1e-6).^2;

skull_c0            = 2820;
skull_rho0          = 1732;
skull_BonA          = 374;
skull_a0            = 2.7 / (source_f0 * 1e-6).^2;

% grid parameters
axial_size          = 120e-3;   % total grid size in the axial dimension [m]
lateral_size        = 80e-3/2;  % total grid size in the lateral dimension [m]

% computational parameters
ppw                 = 50;       % number of points per wavelength
source_x_offset     = 40;       % grid points to offset the source

% --------------------
% GRID
% --------------------

% calculate the grid spacing based on the PPW and F0
dx = water_c0 / (ppw * source_f0);   % [m]

% compute the size of the grid
Nx = roundEven(axial_size / dx) + source_x_offset;
Ny = roundEven(lateral_size / dx);

% create the computational grid
kgrid = kWaveGrid(Nx, dx, Ny, dx);
    
% --------------------
% MEDIUM
% --------------------

% convert geometry to gp
skull_radius_gp             = skull_radius / dx;
skull_thickness_gp          = skull_thickness / dx;
skin_thickness_gp           = skin_thickness / dx;
skin_transducer_distance_gp = skin_tx_distance / dx;

% calculate centers for respective
skin_center = source_x_offset + skin_transducer_distance_gp + skin_thickness_gp + skull_radius_gp;

% create a skin map
skin_mask = makeDisc(Nx * 2, Ny * 6, skin_center, Ny + 1, skull_radius_gp + skin_thickness_gp);

% create a skull map
skull_outer_mask = makeDisc(Nx * 2, Ny * 6, skin_center, Ny + 1, skull_radius_gp);
skull_inner_mask = makeDisc(Nx * 2, Ny * 6, skin_center, Ny + 1, skull_radius_gp - skull_thickness_gp);

% trim masks
skin_mask = skin_mask(1:Nx, Ny + 1:2*Ny);
skull_outer_mask = skull_outer_mask(1:Nx, Ny + 1:2*Ny);
skull_inner_mask = skull_inner_mask(1:Nx, Ny + 1:2*Ny);

% reference sound speed
medium.sound_speed_ref                      = brain_c0;

% sound speed
medium.sound_speed                          = water_c0 * ones(Nx, Ny);
medium.sound_speed(skin_mask == 1)          = skin_c0;
medium.sound_speed(skull_outer_mask == 1)   = skull_c0;
medium.sound_speed(skull_inner_mask == 1)   = brain_c0;

% sound speed
medium.density                              = water_rho0 * ones(Nx, Ny);
medium.density(skin_mask == 1)              = skin_rho0;
medium.density(skull_outer_mask == 1)       = skull_rho0;
medium.density(skull_inner_mask == 1)       = brain_rho0;

% nonlinearity
medium.BonA                                 = water_BonA * ones(Nx, Ny);
medium.BonA(skin_mask == 1)                 = skin_BonA;
medium.BonA(skull_outer_mask == 1)          = skull_BonA;
medium.BonA(skull_inner_mask == 1)          = brain_BonA;

% absorption
medium.alpha_coeff                          = water_a0 * ones(Nx, Ny);
medium.alpha_coeff(skin_mask == 1)          = skin_a0;
medium.alpha_coeff(skull_outer_mask == 1)   = skull_a0;
medium.alpha_coeff(skull_inner_mask == 1)   = brain_a0;

% --------------------
% PLOT
% --------------------

% plot the material property maps
figure;
subplot(1, 4, 1);
imagesc(kgrid.y_vec - kgrid.y_vec(1), kgrid.x_vec - kgrid.x_vec(1), medium.sound_speed);
axis image;
colorbar;
title('Sound Speed');

subplot(1, 4, 2);
imagesc(kgrid.y_vec - kgrid.y_vec(1), kgrid.x_vec - kgrid.x_vec(1), medium.density);
axis image;
colorbar;
title('Density');

subplot(1, 4, 3);
imagesc(kgrid.y_vec - kgrid.y_vec(1), kgrid.x_vec - kgrid.x_vec(1), medium.BonA);
axis image;
colorbar;
title('BonA');

subplot(1, 4, 4);
imagesc(kgrid.y_vec - kgrid.y_vec(1), kgrid.x_vec - kgrid.x_vec(1), medium.alpha_coeff);
axis image;
colorbar;
title('\alpha_0');

scaleFig(1.5, 1);