trisubs.f

linux环境下利用fortran语言开发

                return
        endif

        t0 = t
        t10 = t1
        t20 = t2
        t30 = t3
        t40 = t4

        do20 i=1,3
                x0(i) = x(i)
                x10(i) = x1(i)
                x20(i) = x2(i)
                x30(i) = x3(i)
                x40(i) = x4(i)
                p0(i) = p(i)
                p10(i) = p1(i)
                p20(i) = p2(i)
                p30(i) = p3(i)
                p40(i) = p4(i)
20         continue
      
10      continue

      call ctracr(p0,p10,p20,p30,p40,x0,x10,x20,x30,x40,t0,t10,t20,
     &   t30,t40,j0,amp10,amp20,interface(ipath,nsegment(ipath)-1),
     &   interface(ipath,nsegment(ipath)),mode(ipath,nsegment(ipath)),
     &   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 cendray(amp1,amp2,p0,mode(ipath,nsegment(ipath)),disp,error)

      return
      end



      subroutine croskin(p0,inc,scat,smode,p,error)

      real p0(3), p(3)
      integer inc, scat, smode, error

*
*   This routine continues the ray parameters through a boundary
*
*   INPUT VARIABLES:
*
*       p0(3)   incident ray slowness
*       inc     incident layer index
*       scat    scattering layer index
*       smode   scattered wave type: 1 for SP, 2 for QP, 3 for QS
*
*   OUTPUT VARIABLES:
*
*       p(3)    scattered ray slowness
*       error   0 for no error. 1 if postcritical scattering occurs.
*
*   GLOBAL VARIABLES:
*
*       block "model"   geometric and elastic model parameters.
*
*   EXTERNAL ROUTINES:
*
*       scatslow        real scattering slowness finder
*

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

      real elasti(6), elasts(6), axisi(3), axiss(3), normal(3)
      integer i,ii, is, which(2)

*   initializes some useful variables

      error = 0
      ii = (inc-1)*9
      is = (scat-1)*9

*   load in proper parameters

      do10 i=1,6
        elasti(i) = param(ii+i)
        elasts(i) = param(is+i)
10      continue

      do20 i=1,3
        axisi(i) = param(ii+i+6)
        axiss(i) = param(is+i+6)
20      continue

      which(1) = smode
      normal(1) = 0.
      normal(2) = 0.

*   normal unit vector must be pointing towards incident medium

      if(inc .LT. scat) then
        normal(3) = -1.
        which(2) = 2
      else if(inc .GT. scat) then
        normal(3) = 1.
        which(2) = 2
      else if(p0(3) .GE. 0.) then
        normal(3) = -1.
        which(2) = 1
      else
        normal(3) = 1.
        which(2) = 1
      endif

*   compute the scattered slowness

      call scatslow(p0,normal,axiss,which,elasts,p,error)

      return
      end



      subroutine cscatslow(pinc,iscat,normal,axis,elast,
     &                     pscat1,pscat2,pscat3,error)
      real pinc(3), normal(3), axis(3), elast(6)
      complex pscat1(3), pscat2(3), pscat3(3)
      integer error

*
*   computes the complex scattering slowness of the given mode.
*
*   INPUT VARIABLES:
*   
*       pinc(3)         incident slowness (real)
*       iscat           0 for reflection, 1 for transmission
*       normal(3)       normal unit vector to interface
*                       pointing towards incident medium.
*       axis(3)         anisotropy direction in scattering medium
*       elast(6)        array of elastic coefficients in scattering
*                       medium: A, C, N, L, F, rho
*
*   OUTPUT VARIABLES:
*
*       pscat1(3)       scattered slowness SP mode (complex)
*       pscat2(3)       scattered slowness QP mode (complex)
*       pscat3(3)       scattered slowness QS mode (complex)
*       error           0 if algorithm converged
*                       1 if root finding algorithm did not converge
*                       2 if QP and QS roots could not be separated
*                       3 if triplication occured
*
*   EXTERNAL ROUTINES:
*
*       roots           root finder
*       slowness        real slowness evaluator
*

      real t(3), ptan, ptan2, vtan(3)
      complex del, del1, del2, z1, z2, z3, z4, poly(5), pnorm
      integer iscat, sign, isol, ip, is, i, ireal
      complex pa2, pr2, alpha, beta, zp1, zs1, zp2, zs2, z(4), zz
      real AL, CL, LAF, ca, dot1, dot2, p1(3), p2(3), w1(3), w2(3)
      real A, C, N, L, F, rho, dot, pinc1(3), pinc2(3)
      real pinc3(3), na, ta, na1, ta1, na2, ta2, c2, c1, c0

      error = 0


*   transfers data to more explicit variables

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

*   computes tangential slowness and tangent direction

      dot = 0.e0

      do10 i=1,3
        dot = dot + pinc(i)*normal(i)
10      continue

      ptan2 = 0.e0

      do20 i=1,3
        vtan(i) = pinc(i) - dot*normal(i)
        ptan2 = ptan2 + vtan(i)*vtan(i)
20      continue

      dot = pinc(1)*pinc(1)+pinc(2)*pinc(2)+pinc(3)*pinc(3)     
      ptan = sqrt(ptan2)

*       gets rid of trivial case (normal incidence)

      if(ptan2/dot .EQ. 0.e0) then

        if(iscat .EQ. 0) then
                sign = 1
        else
                sign = -1
        endif

        do40 i=1,3
                pinc1(i) = sign*normal(i)
                pinc2(i) = sign*normal(i)
                pinc3(i) = sign*normal(i)
40         continue

        call slowness(pinc1,axis,1,elast)
        call slowness(pinc2,axis,2,elast)
        call slowness(pinc3,axis,3,elast)

        do50 i=1,3
                pscat1(i) = cmplx(pinc1(i),0.e0)
                pscat2(i) = cmplx(pinc2(i),0.e0)
                pscat3(i) = cmplx(pinc3(i),0.e0)
50         continue

        return

      endif

      do30 i=1,3
        t(i) = vtan(i)/ptan
30      continue

*   computes some usefull quantities for later use

      na = 0.e0

      do60 i=1,3
        na = na + normal(i)*axis(i)
60      continue

      ta = 0.e0

      do70 i=1,3
        ta = ta + t(i)*axis(i)
70      continue

      na2 = na*na
      ta2 = ta*ta

      na1 = 1.e0 - na2
      ta1 = 1.e0 - ta2

*   SP case (solved analytically)

        c2 = L*na2 + N*na1
        c1 = 2.e0*ptan*ta*na*(L-N)
        c0 = ptan2*(L*ta2+N*ta1)-rho

        del = csqrt(cmplx(c1*c1-4.e0*c2*c0,0.e0))
        z1 = (-c1+del)/c2/2.e0
        z2 = (-c1-del)/c2/2.e0

*       determines appropriate solution

        ca = cabs(z1)

        if(ca .EQ. 0.0) then

                pnorm = cmplx(0.,0.)

        else if(abs(aimag(z1))/ca .LE. 1.e-6) then

*       scattering is real: test for scattering direction

                do888 i=1,3
                        p1(i) = vtan(i)+real(z1)*normal(i)
                        p2(i) = vtan(i)+real(z2)*normal(i)
888                continue
                call groupvel(p1,axis,1,elast,w1)
                call groupvel(p2,axis,1,elast,w2)
                dot1 = 0.e0
                dot2 = 0.e0
                do999 i=1,3
                        dot1 = dot1 + w1(i)*normal(i)
                        dot2 = dot2 + w2(i)*normal(i)
999                continue

                if(iscat .EQ. 0) then
                        if(dot1 .GT. 0.e0) then
                                pnorm = z1
                        else
                                pnorm = z2
                        endif
                else
                        if(dot1 .GT. 0.e0) then
                                pnorm = z2
                        else
                                pnorm = z1
                        endif
                endif

        else

*       scattering is complex: test for decay

                if(iscat .EQ. 0) then
                        if(aimag(z1) .GT. 0.e0) then
                                pnorm = z1
                        else
                                pnorm = z2
                        endif
                else
                        if(aimag(z1) .GT. 0.e0) then
                                pnorm = z2
                        else
                                pnorm = z1
                        endif
                endif

        endif

      do100 i=1,3
        pscat1(i) = vtan(i)+pnorm*normal(i)
100      continue


*   QP and QS cases (numerical solution)

*   computes four scattered slownesses (roots of eikonal determinant)

      AL = A*L
      CL = C*L
      LAF = L*L+A*C-(F+L)*(F+L)

      poly(5) = AL*na1*na1 + CL*na2*na2 + LAF*na2*na1
      poly(4) = 2.e0*ptan*na*ta*(LAF*(1.e0-2.e0*na2)
     &          +2.e0*(CL*na2-AL*na1))
      poly(3) = 6.e0*CL*na2*ta2+2.e0*AL*(ta1*na1+2.e0*ta2*na2)
      poly(3) = poly(3)+LAF*(na2+ta2-6.e0*na2*ta2)
      poly(3) = poly(3)*ptan2-rho*((C+L)*na2+(A+L)*na1)
      poly(2) = 2.e0*(CL*ta2-AL*ta1)+LAF*(1.e0-2.e0*ta2)
      poly(2) = poly(2)*ptan2+rho*(A-C)
      poly(2) = poly(2)*2.e0*ptan*na*ta
      poly(1) = CL*ta2*ta2+AL*ta1*ta1+LAF*ta2*ta1
      poly(1) = poly(1)*ptan2-rho*((L+C)*ta2+(L+A)*ta1)
      poly(1) = poly(1)*ptan2+rho*rho

      call roots(poly,z1,z2,z3,z4,isol)

      z(1) = z1
      z(2) = z2
      z(3) = z3
      z(4) = z4

      if(isol .EQ. 0) then
        error = 1
        return
      endif

*   determines number of real roots and reshuffles them 
*   so that if there are any real roots, they are at the 
*   first ranks in the array.

      ireal = 0

      do111 i=1,4
        ca = cabs(z(i))
        if(ca .EQ. 0.e0) then
                ca = 1.e0
        endif
        if(abs(aimag(z(i)))/ca .LE. 1.e-4) then
                zz = z(1)
                z(1) = z(i)
                z(i) = zz
                ireal = ireal + 1
                goto 222
        endif
111      continue

222   do333 i=2,4
        ca = cabs(z(i))
        if(ca .EQ. 0.e0) then
                ca = 1.e0
        endif
        if(abs(aimag(z(i)))/ca .LE. 1.e-4) then
                zz = z(2)
                z(2) = z(i)
                z(i) = zz
                ireal = ireal + 1
        endif
333      continue          

      if(ireal .EQ. 1 .OR. ireal .EQ. 3) then

*   An odd number of real roots, something is wrong.
        error = 2
        return
      endif

      if(ireal .EQ. 4) then

*   There is real QP and QS scattering

      ip = 0
      is = 0

      pa2 = 2.e0*ptan*na*ta*real(z(1))
      pr2 = pa2
      pa2 = pa2+ptan2*ta2+real(z(1))*real(z(1))*na2
      pr2 = ptan2*ta1-pr2+real(z(1))*real(z(1))*na1

      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)

      del1 = cmplx(1.e0,0.e0) - 2.e0*rho/(alpha+beta)
      del2 = cmplx(1.e0,0.e0) - 2.e0*rho/(alpha-beta)

      if(cabs(del1) .LE. cabs(del2)) then
        zp1 = real(z(1))
        ip = ip+1               
      else
        zs1 = real(z(1))
        is = is+1
      endif

      pa2 = 2.e0*ptan*na*ta*real(z(2))
      pr2 = pa2
      pa2 = pa2+ptan2*ta2+real(z(2))*real(z(2))*na2
      pr2 = ptan2*ta1-pr2+real(z(2))*real(z(2))*na1

      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)

      del1 = cmplx(1.e0,0.e0) - 2.e0*rho/(alpha+beta)
      del2 = cmplx(1.e0,0.e0) - 2.e0*rho/(alpha-beta)

      if(cabs(del1) .LE. cabs(del2)) then
        if(ip .EQ. 1) then
                zp2 = real(z(2))
        else
                zp1 = real(z(2))
        endif
        ip = ip+1
      else
        if(is .EQ. 1) then
                zs2 = real(z(2))
        else
                zs1 = real(z(2))
        endif
        is = is+1
      endif

      if(ip .EQ. 2) then
        zs1 = real(z(3))
        zs2 = real(z(4))
      else if(is .EQ. 2) then
        zp1 = real(z(3))
        zp2 = real(z(4))
      else if(is .EQ. 1 .AND. ip .EQ. 1) then
        pa2 = 2.e0*ptan*na*ta*real(z(3))
        pr2 = pa2
        pa2 = pa2+ptan2*ta2+real(z(3))*real(z(3))*na2
        pr2 = ptan2*ta1-pr2+real(z(3))*real(z(3))*na1

        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)

        del1 = cmplx(1.e0,0.e0) - 2.e0*rho/(alpha+beta)
        del2 = cmplx(1.e0,0.e0) - 2.e0*rho/(alpha-beta)

        if(cabs(del1) .LE. cabs(del2)) then
                zp2 = real(z(3))
                ip = ip+1
        else
                zs2 = real(z(3))
                is = is+1
        endif

        if(ip .EQ. 2) then