%% demo_WaveEq
% WaveEq : finite difference acoustic wave simulation toolbox demo
%
%     Copyright (C) Srdan Kitic, Inria Rennes
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU Affero General Public License as
%     published by the Free Software Foundation, either version 3 of the
%     License, or (at your option) any later version.
% 
%     This program 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 Affero General Public License for more details.
% 
%     You should have received a copy of the GNU Affero General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% 
% Contact: nancy.bertin@irisa.fr or remi.gribonval@inria.fr
% PANAMA Research Group,  IRISA  -  (Inria & CNRS)
% Campus de Beaulieu, 35042 Rennes cedex, France
%


clc
clear all

% Setup path
setupwave

%% Initialization

verbose = true; % display output

sizeX = 15; % 15m
sizeY = 15; % 15m
T = 1;      % 1s
sStep = 1;  % spatial discretization stepsize (1m)
tStep = []; % temporal stepsize left to be automatically determined (CFL)

% "Drawing" the shape of domain (a rectangle with Robin boundary cond.)
shape.type = zeros(round([sizeX,sizeY]./sStep));
shape.type([1,end],:) = 4; % "4" indicates Robin BC (check the WaveField.m header)
shape.type(:,[1,end]) = 4;
% Set the boundary parameter for each voxel on the boundary
% (for the Robin BC the parameter "a" is a weight in  "p + a*dp/dn = 0")
shape.param = shape.type;   
shape.param(shape.param ~= 0) = 10; % set it high, to approximate Neumann BC

sources = false; % no sources for now, just create a wavefield object

c = 343; % speed of sound (m/s)

wf = WaveField([sizeX sizeY],T,sStep,tStep,shape,sources,c,verbose);

%% Data

k = 2;                  % number of sources
src_type = wf.wideband; % generate wideband (low-passed white noise sources)
src_cutoff = 0.1*wf.Ft; % source cutoff frequency (wf.Ft is temporal sampling frequency)
start_t = 0.01;         % sources start at 0.01s
active_t = 0.1;         % sources active for 0.1s
magnitude = 1;          % source magnitude

% Generate sources and the wavefield
wf.randomizeSources(k,src_type,start_t,active_t,magnitude,src_cutoff);

% Take some measurements
m = 5;                  % number of microphones
[y,M] = wf.quicksensing(m);

%% Some output

% Wavefield tensor
p = wf.data;
% Source signal tensor
z = wf.sourceData;

% The forward (analysis or D'Alembertian + BC) operator
A = wf.operators{1};

% Verify that Ap = z;
tic
disp(['||Ap-z||=' num2str(norm(A*p(:) - z(:)))])
toc
% Verify that p = inv(A)*z is of O(n) cost
tic
disp(['||p-inv(A)z||=' num2str(norm(p(:) - A\z(:)))])
toc

% Verify that the measurements are correct
disp(['||y-Mp||=' num2str(norm(y-M*p(:)))])

% Animation
disp('Animation (Ctrl+C to break)...')
wf.animate

close(gcf)

%% Restoring old path
setupwave(0)