%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% chap3b.m performs forward modeling (via 2-2 staggered grid FD 1st-order 
% acoustic eqns) to generate line-src response of velocity model for a  line
% src at A and receivers along rectangle boundary. Then it
% backpropagates these data (via 2-2 FD algorithm for 2nd-order acoustic
% wave eqn) back to receivers at A and B using the reciprocity thm of 
% correlation type.

% Theory:  The reciprocity equation of correlation type for positive t>0 values is given by
%
% G(B,t|A,0) = sum_x [dG(A,t|x,0)/dn @ G(x,-t|B,0) - G(A,t|x,0) @ dG(x,-t|B,0)/dn]dx
%
% where @ denotes temporal convolution.  The code is structured so that 
%  1). dG(A,t|x,0)/dn and G(A,t|x,0) are computed by a FD code
%      for a source at A and receivers at x along a square boundary. 
%  2). The dG(A,t|x,0)/dn and G(A,t|x,0) values are time reversed to give dG(A,-t|x,0)/dn
%      and G(A,-t|x,0). These time reversed values will be
%      used by another FD code as the temporal history of sources on the square boundary, to 
%      generate the pressure field at B. This new pressure field at B according to the above 
%      formula should give us G(B,t|A,0).
% 
% (nt,he,wi,c,fr)  - input- (# time steps, height model, width model, velocity, peak frequency)
% (dx,dt,nx,nz)    - input- (gridpoint interval, time sample interval, # model grid pts in x and z directions)
% (nxsrc,nzsrc)    - input- (x,z) coordinates of source at A
% (nxsrc1,nzsrc1)  - input- (x,z) coordinates of source at B
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
%
% Input modeling parameters
%
[nt,he,wi,c,fr,vmax,vmin,invrho,kappa,dx,dt,nx,nz,u,w,px,pz,p,nxsrc,nzsrc,seismo,alpha]=indata1;
nxsrc1=nxsrc;nzsrc1=nzsrc+100; % Locations of A and B in the model interior
indata2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     FORWARD MODELING STEP TO GET G and dG/dn along square boundary
%
% 2-2 FD of the 1st-order acoustic equations of motion to generate data at
% the recatngular boundary for sources at A=(nxsrc,nzsrc) and at
% B=(nxsrc1,nzsrc1) (staggered grid code by Mark Fan of China). Data
% created are G(A|x), dG(A|x)/dn, G(B|x), dG(B|x)/dn
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for is=1:2; % Loop over sources at A and B
%for is=1:1; % Loop over sources at A and B
    xs=nxsrc;zs=nzsrc;xr=nxsrc1;zr=nzsrc1;
    if is==2;
        xs=nxsrc1;zs=nzsrc1;xr=nxsrc;zr=nzsrc;u=u.*0;w=u;px=u;pz=u;p=u;
    end;

    for it=2:nt % Loop over time steps for forward modeling
        u(1:nx-1,:) = (u(1:nx-1,:).* dampx(1:nx-1,:) - kappa1.*diff(p,1,1))./dampxa(1:nx-1,:);
        w(:,1:nz-1) = (w(:,1:nz-1).* dampz(:,1:nz-1) - kappa2.*diff(p,1,2))./dampza(:,1:nz-1);
        px(2:nx,:)  = (px(2:nx,:).* dampx(2:nx,:) - invrho1.*diff(u(:,:),1,1) )./dampxa(2:nx,:);
        pz(:,2:nz)  = (pz(:,2:nz).* dampz(:,2:nz) - invrho2.*diff(w(:,:),1,2) )./dampza(:,2:nz);
        p = px + pz;
        %Adding source
        p(xs,zs)=p(xs,zs)+rick(it);
        store1;plot6;   % Store field values dG(A|x)/dn and G(A|x) along square boundary at x
        pAB(it,is)=p(xr,zr);
    end
    plot5;
end

p0=zeros(nx,nz);p1=p0;p2=p0;bc=zeros(nx,1);
cns=4*(c*dt/dx).^2;
UU=diff(UU);WW=diff(WW);seismo1=seismo(nt,2);
%---------------------------------------------------------------------------
 
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       BACKWARD MODELING STEP TO BACKPROPAGATE FIELD FROM BOUNDARY TO A and B
%
% 2-2 FD of the 2nd-order acoustic wave equation  to backpropagate data
% from the rectangular boundary to receivers at  A=(nxsrc,nzsrc) and at
% B=(nxsrc1,nzsrc1). Reciprocity Eqn is
% G(B|A)=int_B [dG(A,t|x)/dn@G(x,-t|B) - G(x,t|A)@dG(B,-t|x)/dn]dx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for is=1:2
    x=nxsrc1;z=nzsrc1;
    if(is==2)
        x=nxsrc;z=nzsrc;p1=p1.*0.0;p2=p1.*0.0;p0=p1.*0.0;
    end
    for it=nt-2:-1:2 % Loop over reverse time for backpropagation backward into medium

        for ns=1:nsp;p1(xsr(ns),zsr(ns))=UU(it,ns,is)*nxx(ns)+WW(it,ns,is)*nzz(ns);end; % dG/dn at t
        for ns=1:nsp;p0(xsr(ns),zsr(ns))=UU(it+1,ns,is)*nxx(ns)+WW(it+1,ns,is)*nzz(ns);end;% dG/dn at t-1
        p2 = 2*p1 - p0 +  cns.*del2(p1);
        if it==round(it/30)*30;imagesc([1:nx]*dx,[1:nz]*dx,p2');title(num2str(it));pause(0.1);end;
        p0=p1;p1=p2;  %xs=nxsrc;zs=nzsrc;xr=nxsrc1;zr=nzsrc1;


        seismo1(it,is) = p2(x,z);
        %    seismo1(it,2) = p2(x1,z1);
    end;
end
%---------------------------------------------------------------------------



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             PLOT RESULTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nt1]=length(pAB(:,1));
subplot(111);a=diff(pAB(:,1));amx=max(abs(a));plot([0:nt1-1]*dt,[0; a(1:nt1-1)/amx],'-');
%hold on;nt1=length(seismo1(:,2));aa2=seismo1(:,2);mx=max(abs(aa2(:)));
 hold on;nt1=length(seismo1(:,2));aa1=seismo1(:,1);aa2=seismo1(:,2);aa12=aa1+aa2;mx=max(abs(aa12(:)));
plot([0:nt1-1]*dt,aa12/mx,'--');hold off
xlabel('Time (s)');ylabel('Normalized Amplitude');title('G(B,t|A,0)');
text(.05,.6,'Solid  = Analytic');text(.05,.7,'Dashed = Interferometric');
%---------------------------------------------------------------------------