afd_bc_inner.m

基于matlab的反演程序,用于地球物理勘探中射线追踪及偏移成像程序.

function [snapshot]=afd_bc_inner(delx,delt,velocity,snap1,snap2,snapshot,boundary)
% AFD_BC_INNER ... apply absorbin bcs to inner boundary
%
% [snapshot]=afd_bc_inner(nx,nz,delx,delt,velocity,snap1,snap2,snapshot,boundary)
% AFD_BC_INNER applies absorbing boundary conditions to the inner layer of rows
% and columns in a situation where there are 2 layers of rows and columns on
% the boundary.  This will occur when a 9 point approximation to the Laplacian 
% operator is used.  The output matrix has absorbing boundary conditions along its
% inner edge.   In order to be stable, vmax*delt/delx must be < 0.7.  
%
% delx = the grid spacing for horizontal and vertical bins in consistent units
% delt = the time interval in seconds
% velocity = the input velocity matrix in consisnent units
%          = has a size of floor(zmax/delx)+1 by floor(xmax/delx)+1
% snap1 = the wavefield at time=t-2*delt (same size as velocity matrix)
%        = will be based on the source array desired i.e. the position
%          of the sources will be one, and the rest of the positions
%          will be zero
% snap2 = the wavefield at time = t-delt(same size as velocity matrix)
% snapshot = the wavefield at time = t (same size as velocity matrix)
% boundary = indicate whether all sides of the matrix are absorbing
%          = 0 indicates that no absorbing boundaries are desired
%          = 1 indicates all four sides are absorbing
%          = 2 choses three sides to be absorbing, and the top one not to be
%             this enables sources to be put on the surface
%
% snapshot = the wavefield at time = t with absorbing boundary conditions
%
% by Carrie Youzwishen, Febraury 1999 
%
% NOTE: It is illegal for you to use this software for a purpose other
% than non-profit education or research UNLESS you are employed by a CREWES
% Project sponsor. By using this software, you are agreeing to the terms
% detailed in this software's Matlab source file.
 
% BEGIN TERMS OF USE LICENSE
%
% This SOFTWARE is maintained by the CREWES Project at the Department
% of Geology and Geophysics of the University of Calgary, Calgary,
% Alberta, Canada.  The copyright and ownership is jointly held by 
% its author (identified above) and the CREWES Project.  The CREWES 
% project may be contacted via email at:  [e-mail]
% 
% The term 'SOFTWARE' refers to the Matlab source code, translations to
% any other computer language, or object code
%
% Terms of use of this SOFTWARE
%
% 1) Use of this SOFTWARE by any for-profit commercial organization is
%    expressly forbidden unless said organization is a CREWES Project
%    Sponsor.
%
% 2) A CREWES Project sponsor may use this SOFTWARE under the terms of the 
%    CREWES Project Sponsorship agreement.
%
% 3) A student or employee of a non-profit educational institution may 
%    use this SOFTWARE subject to the following terms and conditions:
%    - this SOFTWARE is for teaching or research purposes only.
%    - this SOFTWARE may be distributed to other students or researchers 
%      provided that these license terms are included.
%    - reselling the SOFTWARE, or including it or any portion of it, in any
%      software that will be resold is expressly forbidden.
%    - transfering the SOFTWARE in any form to a commercial firm or any 
%      other for-profit organization is expressly forbidden.
%
% END TERMS OF USE LICENSE

[nz,nx]=size(snap1);

if boundary ==1

  % apply absorbing boundary for top 

  snapshot(2,4:nx-3) = (2.*velocity(2,4:nx-3).*delx*delt.^2)./(delx+velocity(2,4:nx-3).*delt).*...
         ((snapshot(3,4:nx-3)./(2*delt*delx) - snap1(3,4:nx-3)./(2*delx*delt) + ....
         snap1(2,4:nx-3)./(2*delt*delx)) + ...
    (-1./(2.*delt.^2.*velocity(2,4:nx-3))).*...
         (-2.*snap2(2,4:nx-3) + snap1(2,4:nx-3 ) -2.*snap2(3,4:nx-3)+ ...
         snap1(3,4:nx-3) + snapshot(3,4:nx-3)) + ...
    velocity(2,4:nx-3)./(4*delx.^2).*...
         (snapshot(3,5:nx-2) + snap1(2,5:nx-2) + snapshot(3,3:nx-4) - ...
          2.* snapshot(3,4:nx-3) - 2.*snap1(2,4:nx-3) + snap1(2,3:nx-4)));
end

% apply absorbing boundary for bottom

snapshot(nz-1,4:nx-3) = -2.*delx.*delt.^2.*velocity(nz-1,4:nx-3)./(delx+velocity(nz-1,4:nx-3).*delt).*...
       ((-snapshot(nz-2,4:nx-3)./(2*delt*delx) - snap1(nz-1,4:nx-3)./(2*delt*delx) +...
       snap1(nz-2,4:nx-3)./(2*delt*delx)) +...
  1./(2*delt.^2.*velocity(nz-1,4:nx-3)).*...
       (-2.*snap2(nz-1,4:nx-3) + snap1(nz-1,4:nx-3) + snapshot(nz-2,4:nx-3) - ...
       2.*snap2(nz-2,4:nx-3) + snap1(nz-2,4:nx-3)) + ...
  (-velocity(nz-1,4:nx-3)./(4*delx.^2)).* ...
       (snapshot(nz-2,5:nx-2) - 2.*snapshot(nz-2,4:nx-3) + snapshot(nz-2,3:nx-4) + ...
       snap1(nz-1,5:nx-2) - 2.*snap1(nz-1,4:nx-3) + snap1(nz-1,3:nx-4)));

% apply absorbing boundary for right hand side

snapshot(4:nz-3,nx-1) =  -2.*delx.*delt.^2.*velocity(4:nz-3,nx-1)./(delx+velocity(4:nz-3,nx-1).*delt).*...
       ((-snapshot(4:nz-3,nx-2)./(2*delt*delx) - snap1(4:nz-3,nx-1)./(2*delt*delx) + ...
       snap1(4:nz-3,nx-2)./(2*delt*delx)) + ...
  1./(2*delt.^2.*velocity(4:nz-3,nx-1)).*...
       (-2.*snap2(4:nz-3,nx-1) + snap1(4:nz-3,nx-1) + snapshot(4:nz-3,nx-2) - ...
       2.*snap2(4:nz-3,nx-2) + snap1(4:nz-3,nx-2)) + ...
  (-velocity(4:nz-3,nx-1)./(4*delx.^2)).* ...
       (snapshot(5:nz-2,nx-2) - 2.*snapshot(4:nz-3,nx-2) + snapshot(3:nz-4,nx-2) + ...
       snap1(5:nz-2,nx-1) - 2.*snap1(4:nz-3,nx-1) + snap1(3:nz-4,nx-1)));

% apply absorbing boundary for left hand side

snapshot(4:nz-3,2) =(2.*velocity(4:nz-3,2).*delx*delt.^2)./(delx+velocity(4:nz-3,2).*delt).*...
       ((snapshot(4:nz-3,3)./(2*delt*delx) - snap1(4:nz-3,3)./(2*delt*delx) + ...
       snap1(4:nz-3,2)./(2*delt*delx)) + ...
  (-1./(2*delt.^2.*velocity(4:nz-3,2))).* ...
       (-2.*snap2(4:nz-3,2) + snap1(4:nz-3,2) + snapshot(4:nz-3,3) - ...
       2.*snap2(4:nz-3,3) + snap1(4:nz-3,3)) + ...
  (velocity(4:nz-3,2)./(4*delx.^2)).* ...
       (snapshot(5:nz-2,3) - 2.*snapshot(4:nz-3,3) + snapshot(3:nz-4,3) +...
       snap1(5:nz-2,2) - 2.*snap1(4:nz-3,2) + snap1(3:nz-4,2)));

%%%%%% Absorbing boundary conditions for corners

% for lower right hand corner

snapshot(nz-2,nx-1) = velocity(nz-2,nx-1).*delt.*delx./(2.*velocity(nz-2,nx-1).*delt + (2^1/2)*delx).*...
       (snapshot(nz-3,nx-1)./delx + snapshot(nz-2,nx-2)./delx + ...
       (2^1/2)/(velocity(nz-2,nx-1).*delt)*snap2(nz-2,nx-1));

snapshot(nz-1,nx-2) = velocity(nz-1,nx-2).*delt.*delx./(2.*velocity(nz-1,nx-2).*delt + (2^1/2)*delx).*...
       (snapshot(nz-2,nx-2)./delx + snapshot(nz-1,nx-3)./delx + ...
       (2^1/2)/(velocity(nz-1,nx-2).*delt)*snap2(nz-1,nx-2));

snapshot(nz-1,nx-1) = velocity(nz-1,nx-1).*delt.*delx./(2.*velocity(nz-1,nx-1).*delt + (2^1/2)*delx).*...
       (snapshot(nz-2,nx-1)./delx + snapshot(nz-1,nx-2)./delx + ...
       (2^1/2)/(velocity(nz-1,nx-1).*delt)*snap2(nz-1,nx-1));

% for lower left hand corner

snapshot(nz-2,2) = velocity(nz-2,2).*delt.*delx./(2.*velocity(nz-2,2).*delt + (2^1/2)*delx).*...
       (snapshot(nz-3,2)/delx + snapshot(nz-2,3)/delx +...
       (2^1/2)/(velocity(nz-2,2)*delt)*snap2(nz-2,2));

snapshot(nz-1,3) = velocity(nz-1,3).*delt.*delx./(2.*velocity(nz-1,3).*delt + (2^1/2)*delx).*...
       (snapshot(nz-2,3)/delx + snapshot(nz-1,4)/delx +...
       (2^1/2)/(velocity(nz-1,3)*delt)*snap2(nz-1,3));

snapshot(nz-1,2) = velocity(nz-1,2).*delt.*delx./(2.*velocity(nz-1,2).*delt + (2^1/2)*delx).*...
       (snapshot(nz-2,2)/delx + snapshot(nz-1,3)/delx +...
       (2^1/2)/(velocity(nz-1,2)*delt)*snap2(nz-1,2));

%%%%%%%%%%%%%%
if boundary ==1 

  % for upper right hand corner

  snapshot(3,nx-1) = velocity(3,nx-1).*delt.*delx./(2.*velocity(3,nx-1).*delt + (2^1/2)*delx).*...
         (snapshot(4,nx-1)/delx + snapshot(3,nx-2)/delx +...
         (2^1/2)/(velocity(3,nx-1)*delt)*snap2(3,nx-1)); 

  snapshot(2,nx-2) = velocity(2,nx-2).*delt.*delx./(2.*velocity(2,nx-2).*delt + (2^1/2)*delx).*...
         (snapshot(3,nx-2)/delx + snapshot(2,nx-3)/delx +...
         (2^1/2)/(velocity(2,nx-2)*delt)*snap2(2,nx-2));

  snapshot(2,nx-1) = velocity(2,nx-1).*delt.*delx./(2.*velocity(2,nx-1).*delt + (2^1/2)*delx).*...
         (snapshot(3,nx-1)/delx + snapshot(2,nx-2)/delx +...
         (2^1/2)/(velocity(2,nx-1)*delt)*snap2(2,nx-1));

  % for upper left hand corner
  
  snapshot(3,2) = velocity(3,2).*delt.*delx./(2.*velocity(3,2).*delt + (2^1/2)*delx).*...
         (snapshot(4,2)/delx + snapshot(3,3)/delx + ...
         (2^1/2)/(velocity(3,2)*delt)*snap2(3,2));
  
  snapshot(2,3) = velocity(2,3).*delt.*delx./(2.*velocity(2,3).*delt + (2^1/2)*delx).*...
         (snapshot(3,3)/delx + snapshot(2,4)/delx + ...
         (2^1/2)/(velocity(2,3)*delt)*snap2(2,3)); 

  snapshot(2,2) = velocity(2,2).*delt.*delx./(2.*velocity(2,2).*delt + (2^1/2)*delx).*...
         (snapshot(3,2)/delx + snapshot(2,3)/delx + ...
         (2^1/2)/(velocity(2,2)*delt)*snap2(2,2));
end