

% The most simple FD algorithm for acoustic waves you ever find!


close all
clear

nx=200;         % grid points in x
nz=200;         % grid points in z
nt=1000;         % number of time steps
dx=10.;         % grid increment in x
dt=.002;        % Time step
c0=3000.;       % velocity (can be an array 
isx=nx/2;       % source index x
isz=nz/2;       % source index z
ist=10;        % shifting of source time function
f0=100.;         % dominant frequency of source (Hz)
T=1/f0;         % dominant period

% Receiver locations
irx=[60 80 100 120 140];
irz=[5 5 5 5 5];
nr=length(irx);
seis=zeros([nr nt]);

% initialize pressure at different time steps
p=zeros([nx nz]);
pold=p;
pnew=p;
pxx=p;
pzz=p;

% initialize velocity model
c=zeros([nx nz]);
c=c+c0;
% c(nx/2-5:nx/2+5,:)=.8*c(nx/2-5:nx/2+5,:); % this is an example for a
% fault zone

% source time function (1st deriavtive of Gauss)
%for it=1:nt,
%   src(it)=exp(-1/T^2*((it-ist)*dt)^2);
%end
%src=diff(src)/dt;src(nt)=0;
src = zeros(nt);
src(ist) = 1.0;

figure
plot((1:nt)*dt,src)
xlabel(' Time (s) ')
ylabel(' Source amplitude ')


figure
% Time extrapolation

for it=1:nt,   
   % calculate partial derivatives
   for i=2:nx-1,
      for j=2:nz-1,
         pxx(i,j)=(p(i-1,j)-2*p(i,j)+p(i+1,j))/dx^2;
         pzz(i,j)=(p(i,j+1)-2*p(i,j)+p(i,j-1))/dx^2;
      end
   end
       
      % time extrapolation
      
      pnew=2*p-pold+dt*dt*c.*c.*(pxx+pzz);   
      pnew(isx,isz)=pnew(isx,isz)+src(it);
      
      imagesc(pnew')
      xlabel('ix')
      ylabel('iz')
      for ir=1:nr
          text(irx(ir),irz(ir),'+')
      end
      text(isx,isz,'o')
      drawnow
     % surfc(pnew),shading interp, axis off, drawnow
      
      pold=p;
      p=pnew;
      
      % Save seismograms
      for ir=1:nr
          seis(ir,it)=p(irx(ir),irz(ir));
      end
              
end
 
figure
      
% plot seismograms
hold on
for ir=1:nr,
    plot((1:nt)*dt,seis(ir,:)+max(max(seis))*(ir-1))
    xlabel(' Time (s) ')
    ylabel(' Amplitude ')
end
      
% Here it would be good to save the data to a file readible by python

save seismogram seis
