trisubs.f

linux环境下利用fortran语言开发

*       t30     initial traveltime for third tube ray
*       t40     initial traveltime for fourth tube ray
*       amp10   initial ray amplitude for mode
*       amp20   initial ray amplitude for QS mode if mode=1 and 
*               top layer is isotropic. 
*       j0      initial ray jacobian
*
*   GLOBAL VARIABLES:
*
*   common block "model"           contains elastic parameters 
*                               of upper layer.
*

      integer nlayer, triso(100)
      real param(900), geom(0:99)
      common /model/ param, geom, nlayer, triso

      real deg2rad, w0(3), w10(3), w20(3), w30(3), w40(3)
      real azi, dip, g1, g2, dg, sg1, cg1
      real elast(6), axis(3), force(3)
      complex u0(3), cu0(3), cp0(3)
      integer i

*   initializes useful variables

      deg2rad = 1.7453293e-02
      azi = razi*deg2rad
      dip = rdip*deg2rad
      g1 = rg1*deg2rad
      g2 = rg2*deg2rad
      dg = rdg*deg2rad
      cg1 = cos(g1)
      sg1 = sin(g1)

      do10 i=1,6
         elast(i) = param(i)
10    continue

      do20 i=1,3
         axis(i) = param(i+6)
20    continue

*     assemble (unit) force vector 

      force(1) = cos(azi)*cos(dip)
      force(2) = sin(azi)*cos(dip)
      force(3) = sin(dip)

*     computes initial traveltime

      if(mode .EQ. 2) then
         t0 = r*sqrt(param(6)/param(1))
      else
         t0 = r*sqrt(param(6)/param(3))
      endif

      t10 = t0
      t20 = t0
      t30 = t0
      t40 = t0

*   compute slowness direction of central ray

      p0(1) = cg1*cos(g2)
      p0(2) = sg1*cos(g2)
      p0(3) = sin(g2)

*   compute direction of four extra slownesses for ray tube
*   so that the SLOWNESS tube is square, parallel to the 
*   central slowness, and has solid angle aperture dg*dg.

      p10(1) = cg1*cos(g2+dg/2.)
      p10(2) = sg1*cos(g2+dg/2.)
      p10(3) = sin(g2+dg/2.)

      p20(1) = cg1*cos(g2-dg/2.)
      p20(2) = sg1*cos(g2-dg/2.)
      p20(3) = sin(g2-dg/2.)

      p30(1) = p0(1)+(p10(2)-p0(2))*p0(3)-(p10(3)-p0(3))*p0(2)
      p30(2) = p0(2)+(p10(3)-p0(3))*p0(1)-(p10(1)-p0(1))*p0(3)
      p30(3) = p0(3)+(p10(1)-p0(1))*p0(2)-(p10(2)-p0(2))*p0(1)

      p40(1) = p0(1)+(p20(2)-p0(2))*p0(3)-(p20(3)-p0(3))*p0(2)
      p40(2) = p0(2)+(p20(3)-p0(3))*p0(1)-(p20(1)-p0(1))*p0(3)
      p40(3) = p0(3)+(p20(1)-p0(1))*p0(2)-(p20(2)-p0(2))*p0(1)

*   compute group velocity of all rays

      call groupvel(p0,axis,mode,elast,w0)
      call groupvel(p10,axis,mode,elast,w10)
      call groupvel(p20,axis,mode,elast,w20)
      call groupvel(p30,axis,mode,elast,w30)
      call groupvel(p40,axis,mode,elast,w40)

*   compute radiated amplitude

      call freerad(force,p0,axis,elast,mode,t0,u0,x0)

      do60 i=1,3
         cp0(i) = cmplx(p0(i),0.e0)
60    continue

      call cpolar(cp0,axis,mode,elast,cu0)

      amp10 = u0(1)*cu0(1)+u0(2)*cu0(2)+u0(3)*cu0(3)

*   if top layer is isotropic and mode is S, then compute 
*   QS amplitude. 

      if(triso(1) .EQ. 0 .AND. mode .EQ. 1) then

         call freerad(force,p0,axis,elast,3,t0,u0,x0)

         do70 i=1,3
            cp0(i) = cmplx(p0(i),0.e0)
70       continue

         call cpolar(cp0,axis,3,elast,cu0)

         amp20 = u0(1)*cu0(1)+u0(2)*cu0(2)+u0(3)*cu0(3)

      else

         amp20 = cmplx(0.e0,0.e0)

      endif

*   compute slowness vector for all rays

      call slowness(p0,axis,mode,elast)
      call slowness(p10,axis,mode,elast)
      call slowness(p20,axis,mode,elast)
      call slowness(p30,axis,mode,elast)
      call slowness(p40,axis,mode,elast)

*   compute ray starting points

      do30 i=1,3
        x0(i) = w0(i)*t0
        x10(i) = w10(i)*t0
        x20(i) = w20(i)*t0
        x30(i) = w30(i)*t0
        x40(i) = w40(i)*t0
30      continue

*   compute ray jacobian

      j0 = w0(1)*( (x20(2)-x10(2))*(x40(3)-x30(3)) - 
     &   (x20(3)-x10(3))*(x40(2)-x30(2)) )
      j0 = j0 + w0(2)*( (x20(3)-x10(3))*(x40(1)-x30(1)) - 
     &   (x20(1)-x10(1))*(x40(3)-x30(3)) )
      j0 = j0 + w0(3)*( (x20(1)-x10(1))*(x40(2)-x30(2)) - 
     &   (x20(2)-x10(2))*(x40(1)-x30(1)) )
      j0 = abs(j0)

      return
      end



      subroutine cpolar(p,axis,mode,elast,u)
      real axis(3), elast(6)
      complex p(3), u(3)
      integer mode

*
*   computes displacement polarization of given mode
*
*   INPUT VARIABLES:
*
*       p(3)            slowness (complex)
*       elast(6)        elastic coefficients in order:
*                       A, C, N, L, F, rho
*       axis(3)         anisotropy axis direction
*       mode            1 for SP, 2 for QP, 3 for QS
*
*   OUTPUT VARIABLES:
*
*       u(3)            displacement polarization (complex)
*

      complex dot, pr2, pa2, pa(3), pr(3), alpha, beta, epsp
      complex epss, norm, u0(3)
      real A, C, N, L, F
      integer i


      if(mode .EQ. 1) then

*       determining the S-parallel polarization
        
        u(1) = axis(2)*p(3)-axis(3)*p(2)
        u(2) = axis(3)*p(1)-axis(1)*p(3)
        u(3) = axis(1)*p(2)-axis(2)*p(1)

        dot = cmplx(0.e0,0.e0)
      
        do 40 i=1,3
                dot = dot + u(i)*u(i)
40         continue
      
        
        if(cabs(dot) .EQ. 0.e0) then
      
*       ...undetermination case: propagation along axis
      
                u(1) = 0.e0
                u(2) = axis(3)
                u(3) = -axis(2)
                dot = sqrt(axis(2)*axis(2)+axis(3)*axis(3))
                if(cabs(dot) .NE. 0.e0) then
                        u(2) = u(2)/dot
                        u(3) = u(3)/dot
                else
                        u(1) = -axis(3)
                        u(2) = 0.e0
                        u(3) = axis(1)
                        dot = sqrt(axis(1)*axis(1)+axis(3)*axis(3))
                        if(cabs(dot) .NE. 0.e0) then
                                u(1) = u(1)/dot
                                u(3) = u(3)/dot
                        else
                                dot = sqrt(axis(2)*axis(2)+
     &                                  axis(1)*axis(1))
                                u(1) = axis(2)/dot
                                u(2) = -axis(1)/dot
                                u(3) = 0.e0
                        endif
                endif
        else
                dot = csqrt(dot)
      
                do 30 i=1,3
                        u(i) = u(i)/dot
30                 continue        
      
        endif
      
        return
      
      endif
        
*   computes axial and radial slownesses and related quantities

      A = elast(1)
      C = elast(2)
      N = elast(3)
      L = elast(4)
      F = elast(5)

      dot = cmplx(0.e0,0.e0)

      do50 i=1,3
        dot = dot+p(i)*axis(i)
50      continue
      
      do10 i=1,3
        pa(i) = dot*axis(i)
        pr(i) = p(i)-pa(i)
10      continue

      pr2 = cmplx(0.e0,0.e0)
      pa2 = cmplx(0.e0,0.e0)

      do20 i=1,3
        pr2 = pr2+pr(i)*pr(i)
        pa2 = pa2+pa(i)*pa(i)
20      continue

      alpha = (A+L)*pr2 + (C+L)*pa2
      beta = (A-L)*pr2 - (C-L)*pa2
      beta = beta*beta + 4.e0*pr2*pa2*(L+F)*(L+F)
      beta = csqrt(beta)

      if(mode .EQ. 2) then

*       computes Quasi-P polarization
      
        if(cabs(pr2) .LT. 1.e-6*cabs(pa2)) then
*       propagation is along anisotropy axis: polarization is
*       longitudinal.
                u(1) = p(1)
                u(2) = p(2)
                u(3) = p(3)
        else
                epsp = (F+L)*pr2/((alpha+beta)/2.e0-L*pr2-C*pa2)
      
                do80 i=1,3
                        u(i) = pr(i) + epsp*pa(i)
80                 continue
      
        endif
      
        norm = 0.e0
      
        do90 i=1,3
                norm = norm + u(i)*u(i)
90         continue
      
        norm = csqrt(norm)
      
        do 70 i=1,3
                u(i) = u(i)/norm
70         continue
      
        return

      endif

      if(mode .EQ. 3) then

*       computes Quasi-S polarization
      
        if(cabs(pa2) .LT. 1.e-6*cabs(pr2)) then
*       propagation within fracture planes: polarization is
*       along the axis (epss -> - infinity).
        
                u(1) = axis(1)
                u(2) = axis(2)
                u(3) = axis(3)

        else if(cabs(pr2) .LT. 1.e-6*cabs(pa2)) then
*       propagation is along axis: polarization has one degree of
*       freedom, but must be orthogonal to S-parallel so there
*       is no ambiguity after all.

                u0(1) = 0.e0
                u0(2) = axis(3)
                u0(3) = -axis(2)
                dot = sqrt(axis(2)*axis(2)+axis(3)*axis(3))
                if(dot .NE. 0.e0) then
                        u0(2) = u0(2)/dot
                        u0(3) = u0(3)/dot
                else
                        u0(1) = -axis(3)
                        u0(2) = 0.e0
                        u0(3) = axis(1)
                        dot = sqrt(axis(1)*axis(1)+
     &                          axis(3)*axis(3))
                        if(dot .NE. 0.e0) then
                                u0(1) = u0(1)/dot
                                u0(3) = u0(3)/dot
                        else
                                dot = sqrt(axis(2)*axis(2)+
     &                                  axis(1)*axis(1))
                                u0(1) = axis(2)/dot
                                u0(2) = -axis(1)/dot
                                u0(3) = 0.e0
                        endif
                endif
                
                u(1) = -p(2)*u0(3)+p(3)*u0(2)
                u(2) = -p(3)*u0(1)+p(1)*u0(3)
                u(3) = -p(1)*u0(2)+p(2)*u0(1)
                
        else
                epss = (F+L)*pr2/((alpha-beta)/2.e0-L*pr2-C*pa2)

                do110 i=1,3
                        u(i) = pr(i) + epss*pa(i)
110                continue
      
        endif   

        norm = cmplx(0.e0,0.e0)
      
        do130 i=1,3
                norm = norm + u(i)*u(i)
130        continue
      
        norm = csqrt(norm)
      
        do150 i=1,3
                u(i) = u(i)/norm
150        continue

        return

      endif

      return
      end



      subroutine cray(g1,g2,ipath,disp,t,error)

      complex disp(3)
      real g1, g2, t
      integer ipath,error

*
*   traces one raypath from source to receiver
*
*   INPUT VARIABLES:
*
*       g1              ray azimuthal take-off angle (degrees)
*       g2              ray dip take-off angle (degrees)
*       ipath           raypath index
*
*   OUTPUT VARIABLES:
*
*       disp(3)         total complex displacement at receiver
*       t               arrival time at receiver
*       error           0 if no error,
*                       1 if receiver is not illuminated etc...
*
*   GLOBAL VARIABLES:
*
*       block "model"   elastic and geometric model parameters
*       block "array"   source-receiver configuration
*       block "path"    raypath parameters
*       block "output"  graphics/verbose parameters
*       block "io"      input/output logical unit numbers
*
*   EXTERNAL ROUTINES:
*
*       shootmap        creates kinematic interpolation map
*       cinitray        initializes ray parameters at the source
*       ctracr          traces a ray between two interfaces
*       ccrossbound        continues ray parameters through interfaces
*       cendray         incorporates free surface effects to the 
*                       field at the receiver.
*


*      ray color
      parameter (icolor=2)


      real param(900), geom(0:99), dt, scale, xmin, ymin
      real xs, ys, cazi, cdip, cg1, dx, cr, cdg, ail, acl, xr, yr
      integer nlayer, nsegment(200), mode(200,100), tsamp, trace
      integer interface(200,0:100), triso(100), ng, npath, verbose
      integer stdin, stderr, stdout, temp1, temp2, temp3
      common /model/ param, geom, nlayer, triso
      common /path/ mode, interface, nsegment, npath
      common /array/ xs,ys,cazi,cdip,cr,cdg,xr,yr,cg1,
     &               ail,acl,dx,ng,tsamp,dt
      common /output/ verbose, trace, scale, xmin, ymin
      common /io/ stdin, stdout, stderr, temp1, temp2, temp3
      
      real p0(3), p10(3), p20(3), p30(3), p40(3), x0(3), x10(3)
      real x20(3), x30(3), x40(3), t0, t10, t20, t30, t40
      real p(3), p1(3), p2(3), p3(3), p4(3), x(3), x1(3)
      real x2(3), x3(3), x4(3), t1, t2, t3, t4
      real j0, j, xp, yp, zp
      complex amp10, amp20, amp1, amp2
      integer i, iseg


      error = 0

      if(trace .EQ. 1) then
        write(stdout,*) nsegment(ipath)+1, icolor
      endif

      call cinitray(cazi,cdip,mode(ipath,1),g1,g2,cr,cdg,x0,x10,x20,
     &   x30,x40,p0,p10,p20,p30,p40,t0,t10,t20,t30,t40,j0,amp10,amp20)

      if(trace .EQ. 1) then
        xp = xs-xmin
	yp = ys-ymin
	zp = 0. 
        write(stdout,*) xp,yp,zp
      endif

      do10 iseg=1,nsegment(ipath)-1

        call ctracr(p0,p10,p20,p30,p40,x0,x10,x20,x30,x40,t0,t10,
     &          t20,t30,t40,j0,amp10,amp20,interface(ipath,iseg-1),
     &          interface(ipath,iseg),mode(ipath,iseg),x,x1,x2,x3,x4,t,
     &          t1,t2,t3,t4,j,amp1,amp2)

      if(trace .EQ. 1) then
        xp = x(1)+xs-xmin
	yp = x(2)+ys-ymin
        zp = x(3)
        write(stdout,*) xp,yp,zp
      endif

        call ccrossbound(p0,p10,p20,p30,p40,j,amp1,amp2,iseg,ipath,
     &          p,p1,p2,p3,p4,j0,amp10,amp20,error)

        if(error .NE. 0) then