% Example to show recording and re-transmitting a signal using a Dirichlet
% boundary condition
%
% author: Bradley Treeby
% date: 19th July 2015
% last update: 9 September 2016

clear all;

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

% perfectly matched layer
pml_size = 20;
pml_alpha = 2;

% spatial grid
Nx = 256 - 2*pml_size;
dx = 0.1e-3;

% temporal grid
CFL = 0.25;
t_end = 40e-6;

% medium properties
c0 = 1500;

% source
source_freq    = 1.245e6; % [Hz]
source_amp     = 1e6;     % [Pa]
source_cycles  = 15;
source_padding = 100;     % [time points]

% positions
source_pos = 1;           % [grid points]
split_pos = Nx/2;         % [grid points]
record_pos = Nx;          % [grid points]

% width of the recording region
overlap = 1;              % [grid points]  

% options
plot_sim = false;

% =========================================================================
% RUN FULL SIMULATION
% =========================================================================

% create grid
kgrid = makeGrid(Nx, dx);

% define the medium properties
medium.sound_speed = c0;

% define the time array
t_array = makeTime(kgrid, c0, CFL, t_end);
kgrid.t_array = t_array;

% define the source mask
source.p_mask = zeros(Nx, 1);
source.p_mask(source_pos) = 1;

% define the source
source.p = [toneBurst(1/kgrid.dt, source_freq, source_cycles, 'Envelope', 'Rectangular'), zeros(1, source_padding)];
source.p = filterTimeSeries(kgrid, medium, source.p);
source.p = source_amp * source.p ./ max(source.p(:));

% define the sensor mask
sensor.mask = zeros(Nx, 1);
sensor.mask(record_pos, 1) = 1;

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

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

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

% define sensor mask
sensor.mask = zeros(Nx, 1);
sensor.mask(split_pos:split_pos + overlap - 1, :) = 1;

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

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

% take the previously recorded pressure and enforce as a Dirichlet boudnary
% condition at the same position as it was recorded
source.p_mask = zeros(Nx, 1);
source.p_mask(split_pos:split_pos + overlap - 1, :) = 1;
source.p = split_sim_part_1;
source.p_mode = 'dirichlet';

% define the sensor mask
sensor.mask = zeros(Nx, 1);
sensor.mask(record_pos, 1) = 1;

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

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

% plot the signals
figure;
subplot(2, 1, 1);
plot(full_sim * 1e-6, 'b-')
hold on;
plot(split_sim * 1e-6, 'r--')
legend('Full Simulation', 'Split Simulation', 'Location', 'NorthWest');
title('Time Series');
xlabel('time index');
ylabel('pressure [MPa]');

% plot the difference
subplot(2, 1, 2);
plot(100 * abs(full_sim - split_sim) / max(full_sim(:)))
title('Difference');
xlabel('time index');
ylabel('difference [% of max]');