%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Version 4.0 MATLAB   
% 2-2 Finite-Difference Solution to
%
%     2      2 2
%    d P    c d P 
%    --   -   --    =   F
%      2        2
%    dt       dx
%
% with a Free Surface at x=0 and a Ricker Wavelet Source.
% Output data are pressure snapshots and particle acceleration seismograms.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; hold off;
ttotal=1.0;fre=55.;ntrace=30; % total time (s), center frequency Ricker (Hz), # traces

model;  % Compute FD and model parameters

[zz,ricke,ricker,c,r,npt,t]=rickercal(fre,nx,dx,dt); % Generate Ricker wavelet
zz=dx*(1:1:nx);ricker=[ricker zeros(1,nt-length(ricker))];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  1-D Finite Differencing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for it=2:nt;
  p2(2:nx1)=2*p1(2:nx1) - p0(2:nx1) + ...
   alpha*c22(2:nx1).*(p1(1:nx2)-2*p1(2:nx1)+p1(3:nx) );
%
  p2(1)=p2(1)+ricker(it)-ricker(it-1);            % Free Surface Boundary Condition at top z=0;
  p2(nx)=p1(nx)+(p1(nx-1)-p1(nx))*cc(nlay)*dt/dx; % Absorbing Boundary Condition at bottom z=nx;
                                                  % Plot snapshots at every time step
  plot(zz,p2,xxx,yyy,'--');
  axis([0 2000 -.25 1.2]);
  title(['PRESSURE SNAPSHOT  T= ',num2str(dt*it),' sec']);
  xlabel('Depth (ft)'); pause(.01)

  icnt=0;
  for ii=1:20:20*ntrace;icnt=icnt+1;
    seismo(it,icnt)=(p2(ii)-p2(ii+1))/dx;
  end; % Estimate particle acceleration at z=0
    p0=p1; p1=p2;icnt=0;                            % Refresh arrays
end;

figure(2);
mx=max(seismo(:));mn=min(seismo(:));clim=[mn mx];
imagesc([0:ntrace-1]*20*dx,[0:nt-1]*dt,seismo,clim*.07);  % Plot VSP profile
title('VSP Profile');xlabel('Z (ft)');
ylabel('Time (s)');colormap(gray)

 

