% 2D SG PSV 
%clim=[-2e-15 2e-15];
clim=[-1e-13 1e-13];
%
%  S.M. Day, 02 Jan 2009
%
%%Problem Inputs..........................................................
%
%    computational inputs
Nx=200;Nz=200;dx=100;        % x runs from 0 to (Nx-1)*dx; likewise z
CFL=0.5;                     % Courant number for max P speed
Nt=300;                      % Nt*Nplot   = total number of time steps
%    volume inputs
Vp=6000;Vs=3464;rho=2700;    % P speed, S speed, density
%
%.........................................................................
%%Set constants...........................................................
%
%
%     General constants
%dt=CFL*dx/Vp;
dt=0.008;
%dt=0.0119;
lam=rho*(Vp^2-2*Vs^2);
lam2mu=rho*Vp^2;
mu=rho*Vs^2;
%     Fault plane constants
%     Grid coordinates
xgu=(0:Nx-1)*dx;
zgu=(0:Nz-1)*dx;
%
%%Initialize Main Arrays..................................................
%
%
%     volume arrays
u=zeros(Nx,Nz); w=zeros(Nx,Nz);
sxx=zeros(Nx,Nz); szz=zeros(Nx,Nz); 
sxz=zeros(Nx,Nz);
rxx=zeros(Nx,Nz); rzz=zeros(Nx,Nz); 
rxz=zeros(Nx,Nz);
dudx=zeros(Nx-1,Nz); dudz=zeros(Nx,Nz-1);
dwdx=zeros(Nx-1,Nz); dwdz=zeros(Nx,Nz-1);
%
%
%%Time Step (explicit central differences)................................
%
%
alpha=15;for t=1:Nt;s(t) = exp(-alpha*(t*dt-0.7)^2);end;%definition of seismic
%source
for j=1:Nt
j
    time=(j-1)*dt;
%add source
%sxx(100,100)=sxx(100,100)-s(j)*dt/(dx*dx);
%szz(100,100)=szz(100,100)-s(j)*dt/(dx*dx);
sxz(100,100)=sxz(100,100)+s(j)*dt/dx;
%
%    calculate derivatives at interior points
%
dudx(1:Nx-1,1:Nz)=(1/dx)*(u(2:Nx,:)-u(1:Nx-1,:));
dwdx(1:Nx-1,1:Nz)=(1/dx)*(w(2:Nx,:)-w(1:Nx-1,:));
dudz(1:Nx,1:Nz-1)=(1/dx)*(u(:,2:Nz)-u(:,1:Nz-1));
dwdz(1:Nx,1:Nz-1)=(1/dx)*(w(:,2:Nz)-w(:,1:Nz-1));
%
%    calculate stresses at interior points
%
rxx(1:Nx-1,2:Nz)=lam2mu*dudx(1:Nx-1,2:Nz)+lam*dwdz(1:Nx-1,1:Nz-1);
rzz(1:Nx-1,2:Nz)=lam*dudx(1:Nx-1,2:Nz)+lam2mu*dwdz(1:Nx-1,1:Nz-1);
rxz(2:Nx,1:Nz-1)=mu*dudz(2:Nx,1:Nz-1)+mu*dwdx(1:Nx-1,1:Nz-1);
sxx(1:Nx-1,2:Nz)=sxx(1:Nx-1,2:Nz)+dt*rxx(1:Nx-1,2:Nz);
szz(1:Nx-1,2:Nz)=szz(1:Nx-1,2:Nz)+dt*rzz(1:Nx-1,2:Nz);
sxz(2:Nx,1:Nz-1)=sxz(2:Nx,1:Nz-1)+dt*rxz(2:Nx,1:Nz-1);
%
%    update interior velocities
%

u(2:Nx,2:Nz)=u(2:Nx,2:Nz)+...
    (1/rho)*dt*(1/dx)*(sxx(2:Nx,2:Nz)-sxx(1:Nx-1,2:Nz)+...
    sxz(2:Nx,2:Nz)-sxz(2:Nx,1:Nz-1));
%    gamma*dt*(rxx(2:Nx,2:Nz)-rxx(1:Nx-1,2:Nz)+...
%    rxz(2:Nx,2:Nz)-rxz(2:Nx,1:Nz-1)));
w(1:Nx-1,1:Nz-1)=w(1:Nx-1,1:Nz-1)+...
    (1/rho)*dt*(1/dx)*(sxz(2:Nx,1:Nz-1)-sxz(1:Nx-1,1:Nz-1)+...
    szz(1:Nx-1,2:Nz)-szz(1:Nx-1,1:Nz-1));
%    gamma*dt*(rxz(2:Nx,1:Nz-1)-rxz(1:Nx-1,1:Nz-1)+...
%    rzz(1:Nx-1,2:Nz)-rzz(1:Nx-1,1:Nz-1)));

imagesc(xgu,-zgu,u',clim);axis xy;axis image;colorbar;
pause(.1);
end