% Dipole Point Source In A Homogeneous Propagation Medium Example
%
% This example provides a simple demonstration of using k-Wave for the
% simulation and detection of a time varying pressure dipole source within
% a two-dimensional homogeneous propagation medium. It builds on the
% Monopole Point Source In A Homogeneous Propagation Medium Example.
%
% author: Bradley Treeby
% date: 25th August 2010
% last update: 25th August 2010
%  
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009, 2010, 2011 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
% 
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details. 
% 
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>. 

clear all;

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
Nx = 128;           % number of pixels in the x (column) direction
Nz = 128;           % number of pixels in the z (row) direction
dx = 50e-3/Nx;    	% pixel width [m]
dz = dx;            % pixel height [m]
kgrid = makeGrid(Nx, dx, Nz, dz);

% define the properties of the propagation medium
medium.sound_speed = 1500;  % [m/s]
medium.density = 1000;      % [kg/m^3]
medium.alpha_power = 1.5;   % [dB/(MHz^y cm)]
medium.alpha_coeff = 0.75;  % [dB/(MHz^y cm)]

% create the time array
[kgrid.t_array dt] = makeTime(kgrid, medium.sound_speed);

% define a single source element
source.u_mask = zeros(Nz, Nx);
source.u_mask(end - Nz/4, Nx/2) = 1;

% define a time varying sinusoidal velocity source in the z-direction
source_freq = 0.25e6;
source_mag = 2/(medium.sound_speed*medium.density);
source.uz = source_mag*sin(2*pi*source_freq*kgrid.t_array);
source.ux = zeros(size(source.uz));

% smooth the source
source.uz = filterTimeSeries(kgrid, medium, source.uz);

% define a single sensor element
sensor.mask = zeros(Nz, Nx);
sensor.mask(Nz/4, Nx/2) = 1;

% run the simulation
[sensor_data, field_final] = kspaceFirstOrder2D(kgrid, medium, source, sensor, 'PlotScale', [-0.5, 0.5]);

% =========================================================================
% VISUALISATION
% =========================================================================

% plot the final wave-field
figure;
imagesc(kgrid.x(1,:)*1e3, kgrid.z(:,1)*1e3, field_final + source.u_mask + sensor.mask, [-0.5, 0.5]);
colormap(getColorMap);
ylabel('z-position [mm]');
xlabel('x-position [mm]');
axis image;

% plot the simulated sensor data
figure;
[t_sc scale prefix] = scaleSI(max(kgrid.t_array(:)));

subplot(2, 1, 1), plot(kgrid.t_array*scale, source.uz, 'k-');
xlabel(['Time [' prefix 's]']);
ylabel('Signal Amplitude');
axis tight;
title('Input Velocity Signal');

subplot(2, 1, 2), plot(kgrid.t_array*scale, sensor_data, 'r-');
xlabel(['Time [' prefix 's]']);
ylabel('Signal Amplitude');
axis tight;
title('Sensor Pressure Signal');