3-stolt偏移程序-p48.txt

地球物理大师Claebout的几个有用的程序

 #Test case for Stolt migration.
 integer it,nt,ix, real vdtodx; complex cp(256,64)
 open(4,file='plotfile',status='new',access='direct',form='unformatted',recl=1)
  nx=64;nt=256; vdtodx=1./4.      #vdtodx=v dt/dx
  do it=1,nx
     do ix=1,nx
        cp(it,ix)=0.
  cp(32,9)=1.; cp(64,17)=1.; cp(128,33)=1.
  call stolt(nt,nx,cp,vdtodx)
  write(4,rec=1)((real(cp(it,ix)),it=1,nt),ix=1,nx)
  stop;  end

 #Stolt migration subroutine without cosine weight.
 subroutine stolt(nt,nx,cp,vdtodx)
 integer ikx,nx,nt,nth,iktau,iom
 real om, vkx, wl,wh,aktau,pi,pionth,vdtodx
 complex cp(nt,nx),cbf(1025)
 pi=3.14159265; nth=nt/2;pionth=pi/nth;
 call ft2d(nt,nx,cp,1.,-1.,cbf)
 do ikx=1,nx {
     vkx=(ikx-1)*2*pi*vdtodx/nx
     if(ikx>nx/2)vkx=2.*pi*vdtodx-vkx   #negative k
     cbf(1)=0.;     cbf(nt+1)=0         #cbf=working buffer
     do iom=1.nt
         cbf(iom)=cp(iom,ikx)           # Omit weighting
     cp(1,ikx)=0.                       #Ignore Zero freq
     do iktau=2,nth+1  {                #Stretch
         aktau=(iktau-1.01)*pionth
         om=sqrt(aktau*aktau+vkx*vkx); iom=1+om/pionth
         if(iom<nth) {
              wl=iom-om/pionth; wh=1. -wl
              cp(iktau,ikx)=wi*cbf(iom)+wh*cbf(iom+1)
              cp(nt-iktau+2,ikx)=wi*cbf(nt-iom+2)+wh*cbf(nt-iom+1)
              }
         else
              cp(iktau,ikx)=0.
	 }
 }
 call ft2d(nt,nx,cp,-1.,1.,cbf)
 
 relurn; end