
clear
clf

% Simple finite difference solver
%
% Acoustic wave equation  p_tt   = c^2 (p_xx +p_zz) + src
% 2-D regular grid


% Basic parameters

nt=500;    % number of time steps
c0=1.;       % acoustic velocity
eps=.6;     % stability 
isnap=10;    % snapshot sampling


nx=400;     % number of grid points in x 
nz=nx;


p=zeros([nx nz]); pnew=p; pold=p; d2px=p; d2pz=p; 
c=p; c=c+c0;% initialization of pressure and velocity fields
dx=1./(nx-1);  % calculate space increment
x=(0:(nx-1))*dx;
z=(0:(nz-1))*dx;  % initialize space coordinates
dt=eps*dx/max(max(c)); % calculate time step from stability criterion
 
% source time function
f0=1./(10.*dt);                             % dominant frequency
t=(1:nt)*dt; t0=4/f0;                       % initialize time axis
src=exp(-f0^2*(t-t0).*(t-t0));              % source time function
src=diff(diff(src))/dx^2;					% time derivative to obtain Gaussian
                                            % normalized with the cell
                                            % volume
% Source location
xs=nx/2;
zs=nz/2;

% Time stepping

for it=1:nt-2,
      
   % FD
   disp(sprintf(' Time step : %i',it));
   
   for j=2:nx-1,
   
      for k=2:nz-1,
         d2px(j,k)=(p(j+1,k)-2*p(j,k)+p(j-1,k))/dx^2; % space derivative 
         d2pz(j,k)=(p(j,k+1)-2*p(j,k)+p(j,k-1))/dx^2; % space derivative 
      end
   end
   pnew=2*p-pold+c.*c.*(d2px+d2pz)*dt^2;                % time extrapolation
   pnew(xs,zs)=pnew(xs,zs)+src(it)*dt^2;     % add source term
   pold=p;											% time levels
   p=pnew;
   
   if rem(it,isnap)== 0,
   % Display 
   imagesc(x,z,p), axis equal
   colormap gray
   xlabel('x'),ylabel('z')
   title(sprintf(' Time step: %i - Max ampl: %g ',it,max(max(p))))
   drawnow
   end
end

   

