%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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;
nx  = 300; 
nz  = nx/2; 
nt  = 800;  %Define pt src coord.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Doing RTM
shots = 24;
rtm_image = zeros([nz nx]);
if (matlabpool('size')>0)
    matlabpool close;
end
% matlabpool local;
tic;
% parfor is = 1:shots 
for is = 12:12
    xs = round(nx/25)*(is-1) + (nx/25);
    fw = zeros([nz nx]);
    bw = zeros([nz nx]);        
    for it = 3:nt-3
        fdwave = strcat('./fd_data/shot',num2str(xs),'iter', num2str(nt-it),'.dat');%nt-it
        fw = load(fdwave,'-ascii');                      
        bkwave = strcat('./bk_data/bkseismo_shot',num2str(xs),...      
                    'iter',num2str(it),'.dat');   % it
        bw = load(bkwave,'-ascii');
        rtm_image = rtm_image + fw.*bw;   %zero-lag cross correlation     
    end
end
toc;
if (matlabpool('size')>0)
    matlabpool close;
end
imagesc(del2(rtm_image), [-0.006 0.006]); %laplacian filter
colormap(gray);
colorbar;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
