10-总程序p198-200.txt

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

#representations of exp(-z*sqrt(((0.,-1.)*w/v)**2+k**2))

integer output,outfd,fetch
integer iw,nw,ik,nk,omhat,kxhat,kzhat,degree,tfilt,xfilt
real v,dt,dx,dz,xf,x0,tf,tau0,rho,bi,r0,eps,pi,omega,k,vk2
complex cz,cs,cikz,cexp,cmplx,csqrt,cp(1024)
outfd=output()
              call putch("esize","i",8)  #complex numbers
nw=256;       call putch(  "n1","i",nw)  #inner index is 
nk=64;        call putch(  "n2","i",nk)  #outer index is k
              call putch(  "n3","i",1)   #one frame movie

if(fetch(  "v","f", v)==0)  v  =3.754    #rock vel
if(fetch(  "dt,"f", dt)==0) dt =.004     #delta t,sec
if(fetch(  "dx,"f", dx)==0) dx =.025     #delta x,km
if(fetch(  "dz,"f", dz)==0) dz =.004     #delta z,sec
if(fetch(  "xf","f", xf)==0) xf =.25;    x0=xf*nk*dx
if(fetch(  "tf","f", tf)==0) tf =.5;     tau0=tf*nw*dt

if(fetch("omhat","i", omhat)==0) omhat=0 
if(fetch("kxhat","i", kxhat)==0) kxhat=0
if(fetch("kzhat","i", kzhat)==0) kzhat=0
if(fetch("degree","i", degree)==0) degree=90
if(fetch("tfilt","i", tfilt)==0) tfilt=1
if(fetch("xfilt","i", xfilt)==0) xfilt=1

if(fetch(  "rho","i", rho)==0) rho=1-4./nw
if(fetch(  "bi","i", bi)==0) bi=6.726	# b**-1
if(fetch(  "ro","f", ro)==0) r0=0.7071
if(fetch(  "eps","f", eps)==0) eps=0.

call putch("d1","f",dt);   call putch("label1","s","sec")
call putch("d2","f",dx);   call putch("label2","s","kilometers")
call hclose()                    #close date description file 
pi=3.14159265
do ik=1, nk {                    #loop over all 
       k=2*pi*(ik-1.)/nk
       if(k>pi)     k=k-2*pi
       k=k/dx
       if(kxhat==0)
              vk2=(v/2)**2*k*k
       else
              vk2=(v/2)**2*(2/dx)**2*sin(k*dx/2)**2/(1-
                              (4./bi)*sin(k*dx/2)**2)

       do iw=1,nw{                      #loop over all 
             omega=2*pi*(iw-1.)/nw
             if(omega>pi)  omega=omega-2*pi
             omega=omega/dt
             cz=cexp(cmplx(0., omega*dt))
             if(omhat==0)
                cs=cmplx(1.e-5/dt, -omega)
             else
                 cs=(2./dt)*(1.-rho*cz)/(1.+rho*cz)
             if(degree==90)
                cikz=vk2/(csqrt(cs*cs+vk2)+cs)
             if(degree==15 | degree ==45)
                cikz=vk2/(eps+(r0+1.)*cs)
             if(degree==45)
                cikz=vk2/(2.*cs+cikz)
             if(real(cikz)<0.) call erexit("cikz not positive real")

             if(kzhat==0)
                 cp(iw)=cexp(-tau0*cikz)
             else
                 cp(iw)=((1.-cikz*dz/2)/(1.+cikz*dz/2))**(tau0/dz)
             cp(iw)=cp(iw)*cexp(cmplx(0., omega*tau0))   #unretard

             if(tfilt>=1) cp(iw)=cp(iw)*(1+cz)/(1-.8*cz)
             if(tfilt>=2) cp(iw)=cp(iw)*(1-cz)/(1-.8*cz)
             if(xfilt==1) cp(iw)=cp(iw)*(1+cos(k*dx))/(1+.85*
                                 cos(k*dx))
             cp(iw)=cp(iw)*cexp(cmplx(0.,k*x0))
             }
       call rite(outfd,cp,8*nw)  #write
       }
stop;      end

#Finally,you must 2-D Foutier Transform(Section 1.7), take real part,
 and plot.