clear all
%% load Data
Polar = 1 % displacement poralization
load('Deltaevent_030000.mat')
NbSources = 30000  ; %
freqech = 1/dt % 'pulse repetition frequency'
dx_vec = [0 cumsum(dx)] ; dx = dx_vec(2)
dz_vec = [0 cumsum(dz)] ;  dz = dz_vec(2)

Dpl_rs = permute(vz_vx,[2 1 3 4]) ; 
% extract one component of displacement
DPL = squeeze(((Dpl_rs(:,:,Polar,:)))) ; %
% retrieve apparent particle velocities
dtDpl =  diff(DPL,[],3)./dt ; 

% extract the area containting sources and is reflection free
dtDpl_Layer1  = dtDpl(size(dtDpl,1)./2+1:3*size(dtDpl,1)./4-1,...
    size(dtDpl,2)./4:size(dtDpl,2)./4 *3,:) ;  
   
%% 3D cross correlation in space

CrossCorr = zeros(size(dtDpl_Layer1,1)*2-1,size(dtDpl_Layer1,2)*2-1,size(dtDpl_Layer1,3)) ;

% best at around 12, you can see part of wavefront
Ex = 20 ;
zvecEx = round(size(CrossCorr,1)/2)-Ex*8:round(size(CrossCorr,1)/2)+Ex*8 ;
xvecEx =  round(size(CrossCorr,2)/2)-Ex:round(size(CrossCorr,2)/2)+Ex ;



tRange = 5
CrossCorrShow =zeros([size(CrossCorr(zvecEx, xvecEx,:),1),size(CrossCorr(zvecEx, xvecEx,:),2),  tRange*2]) ;

for ptt = 12 
    for trun = 1:size(dtDpl_Layer1,3)
        [CrossCorr(:,:,trun)] = xcorr2(squeeze(dtDpl_Layer1(:,:,ptt)),squeeze(dtDpl_Layer1(:,:,trun)));
    end
    
    for i =  -tRange:tRange
        % ImShow = vit2D(round(size(vit2D,1)/2)-20:round(size(vit2D,1)/2)+20,...
        %     round(size(vit2D,2)/2)-20:round(size(vit2D,2)/2)+20,ptt+i).^2;
        CrossCorrShow(:,:,i+tRange+1)=CrossCorr(zvecEx, xvecEx,ptt+i);
    end
    % figure1(ImShowVec,tRange,dx_vec, xvecEx,dz_vec, zvecEx, ptt,PrC)
    % figure1Plot(ImShowVec,tRange,dx_vec, xvecEx,dz_vec, zvecEx, ptt,PrC)
    
end

%% Visualize

PrC1 = prctile(abs(dtDpl_Layer1(:)),99.5) ;
PrC2 = prctile(abs(CrossCorrShow(:)),99.5) ;
figure,
nexttile, imagesc(dtDpl_Layer1(:,:,2),[-PrC1 PrC1]) %snapshot particle velocity
nexttile, imagesc(dtDpl_Layer1(:,:,3),[-PrC1 PrC1]) %snapshot particle velocity
nexttile, imagesc(CrossCorrShow(:,:,6),[-PrC2 PrC2]) %autocorrelation
nexttile, imagesc(CrossCorrShow(:,:,3),[-PrC2 PrC2]) %snapshot cross-correlation