trisubs.f

linux环境下利用fortran语言开发

         x = xr-xs
         y = yr-ys

         do20 ig=1,ng

            if(verbose .EQ. 1) write(stderr,"('    receiver #',i3)") ig

            call search(xx,yy,x,y,g1,g2,dg2,sg1,sg2,error,
     &                    ng2,maxa)

*           updates raytube angular aperture

            cdg = amin1(ddip/10.e0,dg2)

            if(error .NE. 0) then
               if(verbose .EQ. 1 .AND. error .EQ. 1) then
                  write(stderr,"(' ...receiver #',i3,': offset
     &            too large (out of map).',/)") ig
                  goto 101
               else if(verbose .EQ. 1 .AND. error .EQ. 2) then
                  write(stderr,"(' ...receiver #',i3,':
     &            no real ray to that receiver.',/)") ig
                  goto 101
               endif
            endif

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

            if(error .NE. 0 .AND. verbose .EQ. 1) then
               write(stderr,"(' ...receiver #',i3,
     &         ': undefined scattering.',/)") ig
            endif

            time = t/dt

            tempil = disp(1)*cail+disp(2)*sail
            tempcl = disp(1)*cacl+disp(2)*sacl

            imin = max0(1,time-leng/2)
            imax = min0(tsamp,time+leng/2)


            do100 i=imin,imax

               ii = i-time+(leng+1)/2

               uil(ig,i) = uil(ig,i)+real(tempil)
     &             *wave(ii)+aimag(tempil)*hilwav(ii)

               ucl(ig,i) = ucl(ig,i)+real(tempcl)
     &             *wave(ii)+aimag(tempcl)*hilwav(ii)

               uv(ig,i) = uv(ig,i)+real(disp(3))
     &             *wave(ii) +aimag(disp(3))*hilwav(ii)

100         continue

101         x = x+dxx
            y = y+dyy

20       continue

10    continue


      if(verbose .EQ. 1) 
     &  write(stderr,"(/,'    writing output to files...')")

      open(unit=temp1,file='uh1.dat',status='unknown',
     &     form='unformatted')
      rewind(temp1)
      open(unit=temp2,file='uh2.dat',status='unknown',
     &     form='unformatted')
      rewind(temp2)
      open(unit=temp3,file='uv.dat',status='unknown',
     &     form='unformatted')
      rewind(temp3)

      do30 ig=1,ng
         write(temp1) (uil(ig,j),j=1,tsamp)
         write(temp2) (ucl(ig,j),j=1,tsamp)
         write(temp3) (uv(ig,j),j=1,tsamp)
30    continue

      close(temp1)
      close(temp2)
      close(temp3)

      if(verbose .EQ. 1) write(stderr,"(27x,'...done',/)")
      if(verbose .EQ. 1) 
     &  write(stderr,"('             *** END OF JOB ***',/)")

      return
      end



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

      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 j0, x(3), x1(3), x2(3), x3(3), x4(3), t, t1, t2, t3, t4
      real w0(3), w10(3), w20(3), w30(3), w40(3), j
      complex amp10, amp20, amp1, amp2
      integer mode, start, end

*
*   This routine traces a ray from one interface to the next.
*
*   INPUT VARIABLES:
*
*       p0(3)   ray slowness vector
*       p10(3)  first tube-ray slowness
*       p20(3)  second   "    "    "
*       p30(3)  third    "    "    "
*       p40(3)  fourth   "    "    "
*       x0(3)   ray starting point
*       x10(3)  first tube-ray starting point
*       x20(3)  second    "    "    "   "
*       x30(3)  third     "    "    "   "
*       x40(3)  fourth    "    "    "   "
*       t0      ray initial traveltime
*       t10     first tube-ray initial traveltime
*       t20     second    "    "    "   "
*       t30     third     "    "    "   "
*       t40     fourth    "    "    "   "
*       j0      initial ray jacobian
*       amp10   initial amplitude
*       amp20   initial QS amplitude if medium is isotropic
*               and mode is S (mode=1).
*       start   index of starting interface
*       end     index of destination interface
*       mode    type of wave: 1 for SP, 2 for QP, 3 for QS
*
*   OUTPUT VARIABLES:
*
*       x(3)    ray destination point
*       x1(3)   first tube ray destination point
*       x2(3)   second    "    "    "    " 
*       x3(3)   third     "    "    "    "
*       x4(3)   fourth    "    "    "    "
*       t       ray final traveltime
*       t1      first tube ray final traveltime
*       t2      second    "    "    "    "
*       t3      third     "    "    "    "
*       t4      fourth    "    "    "    "
*       j       final ray jacobian
*       amp1    final amplitude
*       amp2    final QS amplitude if medium is isotropic
*               and mode is S (mode=1).
*   
*   GLOBAL VARIABLES:
*
*       block "model"   geometric and elastic model parameters
*
*   EXTERNAL ROUTINES:
*
*       groupvel        group velocity calculator
*

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

      real axis(3), elast(6), spread
      integer i,ii,layer

      layer = max0(start,end)
      ii = (layer-1)*9

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

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

*   compute the five ray directions

      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 the arrival point for the central ray

      x(3) = geom(end)
      t = t0 + (x(3)-x0(3))/w0(3)
      x(1) = x0(1)+w0(1)*(t-t0)
      x(2) = x0(2)+w0(2)*(t-t0)

*   compute the jacobian and amplitude at the arrival point

*   trace the four extra rays to traveltime t

      do30 i=1,3
        x1(i) = x10(i)+w10(i)*(t-t10)
        x2(i) = x20(i)+w20(i)*(t-t20)
        x3(i) = x30(i)+w30(i)*(t-t30)
        x4(i) = x40(i)+w40(i)*(t-t40)
30      continue

*   compute the jacobian itself from four tube-rays
      
      j = w0(1)*((x2(2)-x1(2))*(x4(3)-x3(3))
     &    -(x2(3)-x1(3))*(x4(2)-x3(2)))
      j = j-w0(2)*((x2(1)-x1(1))*(x4(3)-x3(3))
     &    -(x2(3)-x1(3))*(x4(1)-x3(1)))
      j = j+w0(3)*((x2(1)-x1(1))*(x4(2)-x3(2))
     &    -(x2(2)-x1(2))*(x4(1)-x3(1)))
      j = abs(j)

*   compute the amplitude

      if(j .EQ. 0.) then
        return
      endif

      spread = sqrt(abs(j0/j))
      amp1 = amp10*spread
      amp2 = amp20*spread


*      continue the four extra rays all the way to the interface

      x1(3) = geom(end)
      x2(3) = geom(end)
      x3(3) = geom(end)
      x4(3) = geom(end)

      t1 = t10+(x1(3)-x10(3))/w10(3)
      t2 = t20+(x2(3)-x20(3))/w20(3)
      t3 = t30+(x3(3)-x30(3))/w30(3)
      t4 = t40+(x4(3)-x40(3))/w40(3)

      x1(1) = x10(1)+w10(1)*(t1-t10)
      x1(2) = x10(2)+w10(2)*(t1-t10)
      x2(1) = x20(1)+w20(1)*(t2-t20)
      x2(2) = x20(2)+w20(2)*(t2-t20)
      x3(1) = x30(1)+w30(1)*(t3-t30)
      x3(2) = x30(2)+w30(2)*(t3-t30)
      x4(1) = x40(1)+w40(1)*(t4-t40)
      x4(2) = x40(2)+w40(2)*(t4-t40)

      return
      end



      subroutine ctraction(p,normal,u,axis,elast,traction)
      real normal(3), axis(3), elast(6)
      complex p(3), u(3), traction(3)

*
*   computes the traction traction(3) on the surface element 
*   specified by the normal unit vector normal(3).
*
*   INPUT VARIABLES:
*   
*       p(3)            slowness vector (complex)
*       normal(3)       unit vector normal to traction surface
*       u(3)            displacement vector (complex)
*       elast(6)        elastic parameters in this order:
*                       A, C, N, L, F, rho
*       axis(3)         anisotropy axis unit vector
*
*   OUTPUT VARIABLES:
*
*       traction(3)     traction vector (complex)
*

      real A, C, N, L, F
      real r1(3), r2(3), nprinc(3), dot
      complex epsilon(3,3), sigma(3,3), uprinc(3), pprinc(3)
      complex fprinc(3)
      integer i,j

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

*   Computation of principal axes of anisotropy (i.e. 2 unit vectors 
*    r1(3) and r2(3) to complement axis(3) into a cartesian basis)

      r1(1) = 0.e0
      r1(2) = -axis(3)
      r1(3) = axis(2)
      dot = 0.e0
      
      do 10 i=1,3
        dot = dot+r1(i)*r1(i)
10      continue

      if(dot .NE. 0.e0) then
        dot = sqrt(dot)

        do 20 i=1,3
                r1(i) = r1(i)/dot
20         continue

      else
        r1(1) = -axis(3)
        r1(2) = 0.e0
        r1(3) = axis(1)
        dot = 0.e0
        
        do 30 i=1,3
                dot = dot+r1(i)*r1(i)
30         continue

        if(dot .NE. 0.e0) then
                dot = sqrt(dot)
        
                do 40 i=1,3
                        r1(i) = r1(i)/dot
40                 continue

        else
                r1(1) = axis(2)
                r1(2) = -axis(1)
                r1(3) = 0.e0
                dot = 0.e0

                do 50 i=1,3
                        dot = dot+r1(i)*r1(i)
50                 continue
                
                dot = sqrt(dot)

                do 60 i=1,3
                        r1(i) = r1(i)/dot
60                 continue

        endif
      endif

      r2(1) = axis(2)*r1(3) - axis(3)*r1(2)
      r2(2) = axis(3)*r1(1) - axis(1)*r1(3)
      r2(3) = axis(1)*r1(2) - axis(2)*r1(1)
      dot = 0.e0

      do 70 i=1,3
        dot = dot+r2(i)*r2(i)
70      continue

      dot = sqrt(dot)

      do 80 i=1,3
        r2(i) = r2(i)/dot
80      continue


*   project displacement, unit normal, and slowness onto 
*   principal directions        

      uprinc(1) = u(1)*axis(1) + u(2)*axis(2) + u(3)*axis(3)
      uprinc(2) = u(1)*r1(1) + u(2)*r1(2) + u(3)*r1(3)
      uprinc(3) = u(1)*r2(1) + u(2)*r2(2) + u(3)*r2(3)

      pprinc(1) = p(1)*axis(1) + p(2)*axis(2) + p(3)*axis(3)
      pprinc(2) = p(1)*r1(1) + p(2)*r1(2) + p(3)*r1(3)
      pprinc(3) = p(1)*r2(1) + p(2)*r2(2) + p(3)*r2(3)

      nprinc(1) = normal(1)*axis(1) + normal(2)*axis(2) 
     &            + normal(3)*axis(3)
      nprinc(2) = normal(1)*r1(1) + normal(2)*r1(2) + normal(3)*r1(3)
      nprinc(3) = normal(1)*r2(1) + normal(2)*r2(2) + normal(3)*r2(3)

*   compute deformation tensor epsilon in principal axes 
*   (scaled down by i*omega)

      do 90 i=1,3
        do 100 j=1,3
        epsilon(i,j) = (uprinc(i)*pprinc(j)+uprinc(j)*pprinc(i))/2.
100        continue
90      continue

*   computes stress tensor in principal directions

      sigma(1,1) = C*epsilon(1,1) + F*epsilon(2,2) + F*epsilon(3,3)     
      sigma(2,2) = F*epsilon(1,1) + A*epsilon(2,2) 
     &             + (A-2.e0*N)*epsilon(3,3)    
      sigma(3,3) = F*epsilon(1,1) + (A-2.e0*N)*epsilon(2,2) 
     &             + A*epsilon(3,3)    
      sigma(1,2) = 2.e0*L*epsilon(1,2)
      sigma(1,3) = 2.e0*L*epsilon(1,3)
      sigma(2,3) = 2.e0*N*epsilon(2,3)
      sigma(2,1) = sigma(1,2)
      sigma(3,1) = sigma(1,3)
      sigma(3,2) = sigma(2,3)

*   computes traction in principal directions

      do 110 i=1,3
        fprinc(i) = 0.e0
        do 120 j=1,3
                fprinc(i) = fprinc(i) + sigma(i,j)*nprinc(j)
120        continue
110      continue

*   expresses traction back in reference coordinates

      do 130 i=1,3
        traction(i) = fprinc(1)*axis(i) + fprinc(2)*r1(i) 
     &                + fprinc(3)*r2(i)
130      continue

      return
      end



      subroutine czoeppritz(pi,modei,normal,axisi,elasti,axist,
     &   elastt,coefs,error)
      real pi(3), normal(3), axisi(3), elasti(6), axist(3), elastt(6)
      complex coefs(6)
      integer error, modei

*
*   computes complex scattering coefficients.
*
*   INPUT VARIABLES:
*
*       pi(3)           incident slowness vector (real)
*       modei           incident mode
*       normal(3)       unit normal to interface
*       axisi(3)        anisotropy axis in incident medium
*       elasti(6)       A, C, N, L, F, rho in incident medium
*       axist(3)        anisotropy axis in transmission medium  
*       elastt(6)       A, C, N, L, F, rho in transmission medium
*
*   OUTPUT VARIABLES:
*
*       coefs(6)        complex scattering coefficient in order:
*                       R-SP, R-QP, R-QS, T-SP, T-QP, T-QS
*       error           0 if solutions were found
*                       1 if algorithm did not converge
*
*   EXTERNAL ROUTINES:
*
*       cscatslow       computes complex scattering slownesses
*       cpolar          computes complex polarizations
*       ctraction       computes complex tractions
*       cgausselim      complex linear system solver
*

      integer i
      complex prSP(3), prQP(3), prQS(3), ptSP(3), ptQP(3), ptQS(3)
      complex urSP(3), urQP(3), urQS(3), utSP(3), utQP(3), utQS(3)
      complex frSP(3), frQP(3), frQS(3), ftSP(3), ftQP(3), ftQS(3)
      complex mat(6,6), vec(6), ui(3), fi(3), cpi(3)

      error = 0

*   computation of scattered slownesses from incident slowness

      call cscatslow(pi,0,normal,axisi,elasti,prSP,prQP,prQS,error)

      if(error .NE. 0) then
        return
      endif

      call cscatslow(pi,1,normal,axist,elastt,ptSP,ptQP,ptQS,error)

      if(error .NE. 0) then
        return
      endif

*   compute displacement unit vectors for each wave

      do40 i=1,3
        cpi(i) = cmplx(pi(i),0.e0)
40      continue
      call cpolar(cpi,axisi,modei,elasti,ui)
      call cpolar(prSP,axisi,1,elasti,urSP)
      call cpolar(prQP,axisi,2,elasti,urQP)
      call cpolar(prQS,axisi,3,elasti,urQS)
      call cpolar(ptSP,axist,1,elastt,utSP)
      call cpolar(ptQP,axist,2,elastt,utQP)
      call cpolar(p