junk

seismic software,very useful

1c1
< c******************************************************************************
---
> *******************************************************************************
3,4c3,5
<   	subroutine makeray(nx,ny,nz,fx,fy,fz,dxv,dyv,dzv,vel,vxx,vxy,vxz,
<      1    vyy,vyz,vzz,x0,y0,z0,a0,azh,amin,amax,nr,nt,dt,
---
>   	subroutine makeray(nx,ny,nz,fx,fy,fz,dxv,dyv,dzv,vel,vxx,vxy,
>      1    vxz,vyy,vyz,vzz,x0,y0,z0,a0,azh,amin,amax,nr,nt,dt,
>      5	  vg,vo,del,eps,axa,aya,nag,ndel,neps,dcg,ddel,deps,fdel,feps,
7,29c8
<      4	  rq211,rq212,rq221,rq222,rv,rdvdx,rdvdy,rdvdz,nrs,
<      5    map,v,dvdx,dvdy,dvdz,uxx,uxy,uxz,uyy,uyz,uzz,tzt,
<      6	  x,y,z,px,py,pz,e1x,e1y,e1z,e2x,e2y,e2z,
<      7	  p111,p112,p121,p122,q111,q112,
<      8	  q121,q122,p211,p212,p221,p222,
<      9    q211,q212,q221,q222,dx,dy,dz,
<      a    dpx,dpy,dpz,de1x,de1y,de1z,de2x,
<      b	  de2y,de2z,dp111,dp112,dp121,
<      c    dp122,dq111,dq112,dq121,dq122,
<      d	  dp211,dp212,dp221,dp222,dq211,dq212,
<      e    dq221,dq222,xt,yt,zt,pxt,pyt,pzt,
<      f    e1xt,e1yt,e1zt,e2xt,e2yt,e2zt,
<      g	  p111t,p112t,p121t,p122t,
<      h    q111t,q112t,q121t,q122t,p211t,
<      i	  p212t,p221t,p222t,q211t,
<      j	  q212t,q221t,q222t,dxt,dyt,
<      k	  dzt,dpxt,dpyt,dpzt,de1xt,
<      l	  de1yt,de1zt,de2xt,de2yt,de2zt,
<      m	  dp111t,dp112t,dp121t,dp122t,
<      n	  dq111t,dq112t,dq121t,dq122t,
<      o    dp211t,dp212t,dp221t,dp222t,
<      p	  dq211t,dq212t,dq221t,dq222t,n1)
< 
---
>      4	  rq211,rq212,rq221,rq222,rv,rdvdx,rdvdy,rdvdz,rdt,nrs)
43,45c22,24
< c******************************************************************************
< c*  Trace rays for uniformly sampled vel(nz,ny,nx).
< c******************************************************************************
---
> ******************************************************************************
> c  Trace rays for uniformly sampled v(x,y,z).
> ******************************************************************************
49d27
< c  nz		number of z samples of velocity
52a31
> c  nz		number of z samples of velocity
55a35
> c  vel 		sampled velocity
65,67c45,46
< c               measured from y axis 
< c  amin		lower limit of emergence polar angle 
< c  amax		upper limit of emergence polar angle  
---
> c  amin		lower limit of emergence angle measured from z-axis
> c  amax		upper limit of emergence angle measured from z-axis
71,72c50,66
< c 
< c****************************************************************************
---
>  
> c additional input variables for TI medium 
> c vg 		table of normalized group velocity  --- array(nag,neps,ndel)
> c nag		number of group angle samples in the table 
> c neps		number of epsilon samples in the table 
> c ndel		number of delta samples in the table 
> c dcg		sampling interval of cosin of group angle in the table 
> c deps		sampling interval of epsilon in the table 
> c ddel		sampling interval of delta in the table 
> c feps		first sample of epsilon in the table 
> c fdel		first sample of delta in the table 
> c vo		sampled array(nz,ny,nx) of vertical velocity
> c del		sampled array(nz,ny,nx) of delta
> c eps		sampled array(nz,ny,nx) of epsilon
> c axa		sampled array(nz,ny,nx) of x-directivity of symmetric axis
> c aya		sampled array(nz,ny,nx) of y-directivity of symmetric axis
> ******************************************************************************
85c79,80
< cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
---
> c  rdt				residual traveltime due to anisotropy
> ******************************************************************************
90c85
< c****************************************************************************
---
> *****************************************************************************/
92,101c87,93
< 	integer nrs(nr) 
< 	real vel(nz*ny*nx),vxx(nz*ny*nx),vxy(nz*ny*nx),vxz(nz*ny*nx),
<      1    vyy(nz*ny*nx),vyz(nz*ny*nx),vzz(nz*ny*nx),
<      2	  a0(nr),azh(nr),rx(nt*nr),ry(nt*nr),rz(nt*nr),rpx(nt*nr),
<      3	  rpy(nt*nr),rpz(nt*nr),re1x(nt*nr),re1y(nt*nr),
<      4	  re1z(nt*nr),re2x(nt*nr),re2y(nt*nr),re2z(nt*nr),
<      5	  rq111(nt*nr),rq112(nt*nr),rq121(nt*nr),rq122(nt*nr),
<      6    rp211(nt*nr),rp212(nt*nr),rp221(nt*nr),rp222(nt*nr),
<      7	  rq211(nt*nr),rq212(nt*nr),rq221(nt*nr),rq222(nt*nr),
<      8	  rv(nt*nr),rdvdx(nt*nr),rdvdy(nt*nr),rdvdz(nt*nr)
---
> 	integer nrs(*) 
> 	real vel(*),vxx(*),vxy(*),vxz(*),vyy(*),vyz(*),vzz(*),
>      1	  a0(*),azh(*),rx(*),ry(*),rz(*),rpx(*),rpy(*),rpz(*),
>      2	  re1x(*),re1y(*),re1z(*),re2x(*),re2y(*),re2z(*),
>      3	  rq111(*),rq112(*),rq121(*),rq122(*),rp211(*),rp212(*),
>      4	  rp221(*),rp222(*),rq211(*),rq212(*),rq221(*),rq222(*),
>      5	  rv(*),rdvdx(*),rdvdy(*),rdvdz(*)
109,110c101
< cccc	parameter (n1=128)
< 	integer n1
---
> 	parameter (n1=128)
137,141c128,132
< 	real fxr,lxt,fyr,lyr
< 
<    
< cc	call msgsci(" makeray n1= ",n1)
< cc	call msgsci(" makeray nr= ",nr)
---
> c    additional variabels for anisotropy and perturbation
> 	real vg(*),vo(*),del(*),eps(*),axa(*),aya(*),rdt(*) 
> 	real axa0(n1),aya0(n1),vg0(n1),
>      1	     vo0(n1),del0(n1),eps0(n1),cg0(n1)
> 									   
168d158
< ccc	call msgsci(" loop 5 nr= ",nr)
204d193
< ccc	call msgsci(" loop 10 nr= ",nr)
213d201
< ccc	call msgsci(" loop 20 nr= ",nr)
240d227
< ccc	call msgsci(" loop 50 nr= ",nr)
270a258
> 		rdt(ir) = 0.
272,273d259
< 
< ccc	call msgsci(" after loop 50 nr= ",nr)
290,294d275
< c       boundary tolerance
< 	fxr = fx - 0.001*dxv
< 	lxr = lx + 0.001*dxv
< 	fyr = fy - 0.001*dyv
< 	lyr = ly + 0.001*dyv
296,304d276
< c	call msgscf("tzmin=",tzmin)
< c	call msgscf("tzmax=",tzmax)
< c	call msgscf("fx=",fx)
< c	call msgscf("lx=",lx)
< c	call msgscf("fy=",fy)
< c	call msgscf("ly=",ly)
< c	call msgscf("fz=",fz)
< c	call msgscf("lz=",lz)
< 
311,317d282
< 
< cccc round up to the boudaries if x,y,z are very close 
< 		if(x(ir).lt.fx .and. x(ir).gt.fxr) x(ir) = fx
< 		if(x(ir).gt.lx .and. x(ir).lt.lxr) x(ir) = lx
< 		if(y(ir).lt.fy .and. y(ir).gt.fxr) y(ir) = fy
< 		if(y(ir).gt.ly .and. y(ir).lt.lyr) y(ir) = ly
< 
322,331d286
< 
< c		if(x(ir).lt.fx) call msgscf("x=",x(ir))
< c		if(x(ir).gt.lx) call msgscf("x=",x(ir))
< c		if(y(ir).lt.fy) call msgscf("y=",y(ir))
< c		if(y(ir).gt.ly) call msgscf("y=",y(ir))
< c		if(z(ir).lt.fz) call msgscf("z=",z(ir))
< c		if(z(ir).gt.lz) call msgscf("z=",z(ir))
< c		if(tzt(ir).lt.tzmin) call msgscf("tzt=",tzt(ir))
< c		if(tzt(ir).gt.tzmax) call msgscf("tzt=",tzt(ir))
< 
517a473,476
> 
> c	interpolate anisotropic parameters
>  	    call anip3Interp(nx,ny,nz,fx,fy,fz,dxv,dyv,dzv,
>      1    	vo,del,eps,axa,aya,x,y,z,nr,vo0,del0,eps0,axa0,aya0)
518a478,489
> c	compute group angle
> 
> 	    do 490 ir=1,nr
>    		cg0(ir) = v(ir)*(px(ir)*axa0(ir)+py(ir)*aya0(ir)+pz(ir)
>      1			 *sqrt(1-axa0(ir)*axa0(ir)-aya0(ir)*aya0(ir)))
> 		cg0(ir) = abs(cg0(ir))
> 490	    end do
> 
> c	look up group velocity table
>  	    call vg3Interp(ndel,neps,nag,fdel,feps,1.0,ddel,deps,-dcg,
>      1    	vg,del0,eps0,cg0,nr,vg0)
>     
549a521,523
> c	      calculate traveltime perturbation due to anisotropy
>  		rdt(it*nrsave+jr) = (v(ir)/(vg0(ir)*vo0(ir))-1.0)*dt
>