%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Version 4.0 MATLAB   
% 2-2 Finite-Difference Solution to elastic wave equation
% using a staggered grid method.
%
% with a Free Surface at x=0 and a Ricker Wavelet Source.
% Output data are pressure snapshots and particle acceleration seismograms.
% Stability is vmax dt/dx<.6 and dispersion is about 12 pts per minimum
% wavelength ..which is that of the shear wave. Make sure any near surface source is
% placed at the 3rd grid point below the top node.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; hold off; 
nx=320;nz=120; fr=50.; 
c=ones(nx,nz)*1500;for j=round(nz/2):nz;c(:,j)=c(:,j)*1.3;end;
[nzm,nxm,nt,tmin,dx,p2,p1,p0,seismo,cns,nsp,xxx,yyy,dt,alpha,cm,ca]=param(nx,nz,c,fr) ;
dens=ones(nx,nz);u1=dens*0;px1=u1;v1=u1;pz1=u1;xy1=u1;
cl=ca-2*cm;
%
%    Define Ricker Wavelet
%
[r,ricker]=calrick(dt,tmin,fr);dtx=dt/dx;
%
nt2=round(nt/5);
seismo=zeros(nx,nt2);seismoy=zeros(nx,nt2);seismotzz=zeros(nx,nt2);
for it=2:nt2
%
% Interior 2-2 FD of the Wave Equation (cns = 4*{c dt/dx}^2 )
%
it
k=[2:nz-1];
j=[2:nx-1];
u1(j,k)=u1(j,k)+ dtx./dens(j,k).*(px1(j,k)-px1(j-1,k)+xy1(j,k)-xy1(j,k-1));
%
v1(j,k)=v1(j,k)+ dtx./dens(j,k).*(pz1(j,k+1)-pz1(j,k)+xy1(j+1,k)-xy1(j,k));
%
% Compute Stress Components
%
au=u1(j+1,k)-u1(j,k);
av=v1(j,k)-v1(j,k-1);
px1(j,k)=px1(j,k)+dtx*( ca(j,k).*au+cl(j,k).*av);
pz1(j,k)=pz1(j,k)+dtx*( cl(j,k).*au+ca(j,k).*av);
xy1(j,k)=xy1(j,k)+(cm(j,k)*dtx).*(u1(j,k+1)-u1(j,k)+v1(j,k)-v1(j-1,k));
%  
% Free-surface boundary condition at the top boundary
%
k=1;
for j=1:nx;
  xy1(j,k)=0.;
  pz1(j,k)=-pz1(j,k+1);
end;
%
% Add Source Term at nsp points to vertical particle velocity
%
if it<=r;
    for ns=1:nsp;
        v1(xxx(ns),yyy(ns))=v1(xxx(ns),yyy(ns))+ricker(it);
    end;
end
%
% Save Particle Velocity and tzz Seismograms at Free Surface.
%
seismo(1:nx,it)=u1(1:nx,3);seismoy(1:nx,it)=v1(1:nx,3);
seismotzz(1:nx,it)=pz1(1:nx,3);
if it/30 == round(it/30); plotit(dx,pz1+px1,nx,nz,it,dt); end   
it
end;
%
% Plot Seismograms 3 pts below free surface
%
fac=.1;
ymx=max(abs(seismo(:)));clim=[-fac*ymx,fac*ymx];
figure(2);imagesc([1:nx]*dx,[1:nt2]*dt,seismo',clim);xlabel('X (m)');ylabel('Time (s)');
title('Common Shot Gather (horizontal component phone)');colorbar
ymx=max(abs(seismoy(:)));clim=[-fac*ymx,fac*ymx];
figure(3);imagesc([1:nx]*dx,[1:nt2]*dt,seismoy',clim);xlabel('X (m)');ylabel('Time (s)');
title('Common Shot Gather (vertical component phone)');colorbar
ymx=max(abs(seismotzz(:)));clim=[-fac*ymx,fac*ymx];
figure(4);imagesc([1:nx]*dx,[1:nt2]*dt,seismotzz',clim);xlabel('X (m)');ylabel('Time (s)');
title('Common Shot Gather (tauzz component phone)');colorbar
figure(5);imagesc([1:nx]*dx,[1:nz]*dx,c');xlabel('X (m)');ylabel('Z (m)');
title('Velocity Model (m/s)');colorbar

