%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2-4 DF Soln to ACOUSTIC WAVE EQUATION (cnst density)
% (nx,nz,nt)      - input- (Horizontal,Vertical) gridpt dimens. of vel 
%                            model & # Time Steps
% FRE             - input- Peak frequency of Ricker wavelet
% BVEL            - input- NXxNZ matrix of background velocity model
% (dx,dt)         - input- (space, time) sample intervals
% (xs,zs)         - input- (x,z) coordinates of line source
% RICKER(nt)      - input- nt vector of source time histories
% (p2,p1,p0)      -calcul- (future,present,past) NXxNZ matrices of 
%                           modeled pressure field
% (p0,p1)         -output- Old and present pressure panels at time nt.
% REALDATA(nx,nt) -output- CSG seismograms at z=2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
c   = 4.0; %km/s
FRE = 20; 
nx  = 300; 
nz  = nx/2; % Define model size, src. freq., velocity
dx  = c/FRE/20;  % Set spatial sampling dx
dz  = dx;
dt  = 0.6*(sqrt(dx*dx+dz*dz)/2.0)/c; %MUST BE THE SAME AS IN THE FORWARD MODELING

nt   =  800;
nxm  =  nx-1;
nzm  =  nz-1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%show the Ricker wavelet
t       =   [0:1:nt-1]*dt - 0.95/FRE;
RICKER  =   zeros(length(t));
%RICKER(2) = 1;
RICKER  =   (1 - t.*t * FRE^2 *pi^2  ) .*exp(- t.^2 * pi^2 * FRE^2 );

plot([0:nt-1]*dt, RICKER);
title('Ricker Wavelet');
xlabel('Time (s)');
pause(1.5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Testing the reciprocal
input  = [0:1:nt-1]*0; %receiver becomes source
output = zeros([nt 1]); %source becomes receiver
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
BVEL  =  ones(nz,nx)*c;
BVEL(nz-round(nz/2):nz, :) = BVEL(nz-round(nz/2):nz, :)*1.4; % 2-Layer velocity model
% velocity in the 2nd layer is 1.4 time of the first layer above it
REALDATA  = zeros(nx,nt); % Define shot gather matrix
% p0 = zeros(nz,nx);
% p1 = p0;
% p2 = p0;

cns = (dt/dx*BVEL).^2;
alpha = dt/dx;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Doing bk_propagate
zs = 2;
zr = 2;


tic;
%parfor is = 1:24     
for is  = 12:12
    xs = (nx/25)*(is-1) + (nx/25);

    fd_seis_path = strcat('./fd_data/seismo_shot',num2str(xs),'.dat');
    fd_seismo    = load(fd_seis_path,'-ascii');       
    
    p0 = zeros(nz,nx);
    p1 = p0;
    p2 = p0;
    p0(zs,:) = fd_seismo(:,nt)';
    p1(zs,:) = fd_seismo(:,nt-1)';
    laplace_p1 = zeros([nz nx]);
    
    for it = nt-2:-1:1
        laplace_p1 = Laplacian(p1,nz,nx,p0,alpha,BVEL);
        p2 = 2*p1 - p0 + cns.*laplace_p1;                    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %POINT SOURCEs at receivers, back propagate
        p2(zs,:) = p2(zs,:) + (fd_seismo(:,it))';
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % ABC BOUNDARY CONDITIONS
        % if p2(1, :) == 0, then it is a free boundary surface, row = 1,         
        %p2(1, :)  = p1(1,:) - alpha*BVEL(1,:).*(p1(1,:) - p1(2,:));        % Top ABC, just perpendicular to boundary       
        p2(1, :)   = 0;                                                     % free surface
        p2(nz, :)  =  p1(nz, :) - alpha*BVEL(nz, :).*(p1(nz,:) - p1(nzm,:));% Bottom ABC
        %p2(nz, :) =  p1(nz, :);                                           % Bottom reflecting
        p2(:, nx)  = p1(:, nx) - alpha*BVEL(:, nx).*(p1(:,nx) - p1(:,nxm)); % RHS ABC                
        p2(:, 1)   = p1(:, 1) + alpha*BVEL(:, 1).*(p1(:, 2) - p1(:, 1));    % LHS ABC
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
        REALDATA(:, nt-it) = p2(zr,:)';
        p0 = p1;
        p1 = p2; %update pressure fields
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        bkwavefield = strcat('./bk_data/bkseismo_shot',num2str(xs),...
                'iter', num2str(it), '.dat');
        save_survey(bkwavefield,p2);
        if (it/10 == round(it/10))            
            plt;pause(0.01);           
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end    
end
toc;
