%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2-2 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;
c=4.0;FRE=20;nx=300;nz=nx; % Define model size, src. freq., velocity
dx=c/FRE/20;  % Set spatial sampling dx

dt=.5*dx/c;   % Set dt for stability condition
xs=round(nx/2.3); zs=round(nx/2);nt=600; % Define pt src coord.
nxm=nx-1;nzm=nz-1; 

t=[0:1:nt-1]*dt-0.95/FRE;RICKER=zeros(length(t));
RICKER= (1-t .*t * FRE^2 *pi^2  ) .*exp(- t.^2 * pi^2 * FRE^2 ) ;
plot([0:nt-1]*dt,RICKER);title('Ricker Wavelet');xlabel('Time (s)')
pause(3)

BVEL=ones(nx,nz)*c;
BVEL(nx-round(nx/2):nx,:)= BVEL(nx-round(nx/2):nx,:)*1.4; % 2-Layer velocity model
REALDATA=zeros(nx,nt); % Define shot gather matrix
p0=zeros(nx,nz);p1=p0;p2=p0;
cns=4*(dt/dx*BVEL).^2; 
recdepth=7;            % Receiver line 7 grid pts below source
% LOOP OVER TIME STEPS
for it=1:1:nt
   p2 = 2*p1 - p0 +  cns.*del2(p1);    % FD of Acoustic WE
   p2(xs,zs) = p2(xs,zs) + RICKER(it); % Add bodypoint src term
   a = p2(xs+recdepth,:); REALDATA(:,it) = a'; % Save seismograms
   alpha=dt/dx;
   
%  ABC BOUNDARY CONDITIONS
   p2(1,:)=p1(1,:)   + alpha*BVEL(1,:).*(p1(2,:)-p1(1,:));     %Top ABC
   p2(nx,:)=p1(nx,:) - alpha*BVEL(nx,:).*(p1(nx,:)-p1(nxm,:)); % Bottom ABC
   p2(:,nz)=p1(:,nz) + alpha*BVEL(:,nz).*(p1(:,nzm)-p1(:,nz)); % RHS ABC

   p0=p1;p1=p2; % UPDATE PRESSURE FIELDS
   if round(it/20)*20==it;plt;pause(1);end % PLOT OUT RESULTS
end;

