% Example to show recording and re-transmitting a signal using a Dirichlet
% boundary condition in 3D.
%
% author: Bradley Treeby
% date: 30th November 2017
% last update: 6th June 2020

clearvars;

% =========================================================================
% DEFINE LITERALS
% =========================================================================

% spatial grid (the PML is outside the grid)
Nx = 50;
Ny = 50;
Nz = Ny;
dx = 0.1e-3;
dy = dx;
dz = dx;

% temporal grid
CFL = 0.15;
t_end = 20e-6;

% medium properties
c0 = 1500;
rho0 = 1000;

% source
source_width    = 20;   % [grid points]
source_freq     = 1e6;  % [Hz]
source_amp      = 1e6;  % [Pa]
source_phase    = 0;    % [rads]

% number of periods to record
record_periods = 3;

% =========================================================================
% RUN FULL SIMULATION (FOR REFERENCE ONLY)
% =========================================================================

% create grid
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);

% define the medium properties
medium.sound_speed = c0;
medium.density = rho0;

% set the time step based on the CFL
dt = CFL * dx / c0;

% recalculate dt to force time sampling to be an integer number of points
% per period
period = 1 / source_freq;
points_per_period = round(period / dt);
dt = period / points_per_period;

% calculate number of time steps
Nt = round(t_end / dt);

% assign time array
kgrid.setTime(Nt, dt);

% calculate corresponding CFL, PPW, and PPP for reference
disp(['CFL = ' num2str(c0 * dt / dx)]);
disp(['PPP = ' num2str(1 / (kgrid.dt * source_freq))]);
disp(['PPW = ' num2str(c0 / (dx * source_freq))]);

% define the source mask
source.p_mask = zeros(Nx, Ny, Nz);
source.p_mask = makeBowl([Nx, Ny, Nz], [1, Ny/2, Nz/2], 50, 41, [Nx, Ny/2, Nz/2]);

% define the source as a sinusoid
source.p = createCWSignals(kgrid.t_array, source_freq, source_amp, source_phase);

% set the sensor mask
sensor.mask = zeros(Nx, Ny, Nz);
sensor.mask(1) = 1;
sensor.record = {'p_max_all'};

% only record the last few periods in steady state
sensor.record_start_index = kgrid.Nt - record_periods * points_per_period + 1;

% set the input arguments
input_args = {'PMLInside', false, 'PMLSize', 'auto',...
    'PlotPML', false, 'PlotScale', [-1, 1] * source_amp};

% run the simulation
sensor_data_all = kspaceFirstOrder3DG(kgrid, medium, source, sensor, input_args{:});

% =========================================================================
% RUN SPLIT SIMULATION PART 1
% =========================================================================

% create grid
kgrid = kWaveGrid(Nx/2, dx, Ny, dy, Nz, dz);

% assign time array, running the simulation for half the time
kgrid.setTime(Nt/2, dt);

% redefine the source mask
source.p_mask = source.p_mask(1:Nx/2, :, :);

% define sensor mask to record across final plane in the simulation
sensor.mask = zeros(Nx/2, Ny, Nz);
sensor.mask(end, :, :) = 1;
sensor.record = {'p', 'p_max_all'};

% only record the last few periods in steady state
sensor.record_start_index = kgrid.Nt - record_periods * points_per_period + 1;

% run the simulation
sensor_data_split_1 = kspaceFirstOrder3DG(kgrid, medium, source, sensor, input_args{:});

% extract amplitude from the sensor data
[amp, phase] = extractAmpPhase(sensor_data_split_1.p, 1/kgrid.dt, source_freq, ...
    'Dim', 2, 'FFTPadding', 1, 'Window', 'Rectangular');

% =========================================================================
% RUN SPLIT SIMULATION PART 2
% =========================================================================

% take the previously recorded pressure and velocity and enforce as a
% Dirichlet boudnary condition
source.p_mask = zeros(Nx/2, Ny, Nz);
source.p_mask(1, :, :) = 1;
source.p = createCWSignals(kgrid.t_array, source_freq, amp, phase);
source.p_mode = 'dirichlet';

% set the sensor mask
sensor.mask = zeros(Nx/2, Ny, Nz);
sensor.mask(1) = 1;
sensor.record = {'p_max_all'};

% run the simulation
sensor_data_split_2 = kspaceFirstOrder3DG(kgrid, medium, source, sensor, input_args{:});

% =========================================================================
% DISPLAY
% =========================================================================

% full simulation
full_sim = sensor_data_all.p_max_all(1:Nx-1, :, Nz/2);

% split simulation
split_sim = [sensor_data_split_1.p_max_all(:, :, Nz/2); sensor_data_split_2.p_max_all(2:end, :, Nz/2)];

% plot
figure;
subplot(1, 3, 1);
imagesc(full_sim * 1e-6)
axis image
title('Full Simulation');
h = colorbar;
title(h, '[MPa]');

subplot(1, 3, 2);
imagesc(split_sim * 1e-6)
axis image
title('Split Simulation');
h = colorbar;
title(h, '[MPa]');

subplot(1, 3, 3);
imagesc(100 * abs(full_sim - split_sim) / max(full_sim(:)))
axis image
title('Difference');
h = colorbar;
title(h, '[%]');