trisubs.f

linux环境下利用fortran语言开发

* Copyright (c) Colorado School of Mines, 2007.
* All rights reserved.

*****************************************************************
*								*
*	Trisubs: Routines for TRISO				*	
*								*
*	Author: Sebastien Geoltrain				*
*								*
*	Copyrights: Center for Wave Phenomena,			*
*		    Mathematics Department,			*
*		    Colorado School of Mines,			*
*		    Golden, CO 80401.				*
*								*
*	All Rights Reserved.					*
*								*
*****************************************************************

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

      complex amp10, amp20, amp1, amp2
      real p0(3), p10(3), p20(3), p30(3), p40(3), j0
      real p(3), p1(3), p2(3), p3(3), p4(3), j
      integer error, iseg, ipath

*
*   This routine continues the ray parameters from layer number
*   "inc" to layer number "scat".
*
*   INPUT VARIABLES:
*
*       p0(3)   incident ray slowness
*       p10(3)  first tube-ray incident slowness
*       p20(3)  second tube-ray incident slowness
*       p30(3)  third tube-ray incident slowness
*       p40(3)  fourth tube-ray incident slowness
*       j0      incident ray jacobian
*       amp10   incident amplitude
*       amp20   incident amplitude of QS mode if incident
*               medium is isotropic and mode is S.
*       iseg    index of incident ray segment.
*       ipath   index of raypath
*
*   OUTPUT VARIABLES:
*
*       p(3)    scattered ray slowness
*       p1(3)   first tube-ray scattered slowness
*       p2(3)   second tube-ray scattered slowness
*       p3(3)   third tube-ray scattered slowness
*       p4(3)   fourth tube-ray scattered slowness
*       j       scattered ray jacobian
*       amp1    scattered amplitude
*       amp2    scattered amplitude of QS mode if scattering
*               medium is isotropic and scattering mode is S.
*       error   0 for no error.
*               1 if postcritical scattering occurs.
*               2 if jacobian is zero (caustic)
*               3 if scattering coefficients are not computed
*
*   GLOBAL VARIABLES:
*       
*       block "model"   elastic and geometric parameters
*       block "path"    raypath parameters
*
*   EXTERNAL ROUTINES:
*
*       scatslow        real scattering slowness finder
*       groupvel        real group velocity evaluator
*       czoeppritz      complex zoeppritz system solver
*

      integer nlayer, triso(100), nsegment(200)
      integer mode(200,100), interface(200,0:100), npath
      real param(900), geom(0:99)
      common /model/ param,geom,nlayer,triso
      common /path/ mode, interface, nsegment, npath

      integer inc, scat, trans, imode, smode
      real elasti(6), elasts(6), elastt(6), axisi(3), axiss(3)
      real dot, dot0, w(3), w0(3), axist(3), normal(3)
      complex coefs(6), coefs2(6)
      integer i, iinc, iscat, isol, isol1, isol2, isol3, isol4
      integer which(2), itrans

*   initializing some variables

      error = 0
      
*   determining incident and scattering modes, incident, 
*   scattering and transmission layer indexes.
      
      imode = mode(ipath,iseg)
      smode = mode(ipath,min0(iseg+1,nsegment(ipath)))

      inc = max0(interface(ipath,iseg-1),interface(ipath,iseg))
      scat = max0(interface(ipath,iseg),interface(ipath,min0(iseg+1,
     &          nsegment(ipath))))

      if(inc .NE. scat) then
*   ...transmission
        trans = scat
      else if(interface(ipath,iseg-1) .LT. interface(ipath,iseg)) then
*   ..."v" type reflection
        trans = interface(ipath,iseg)+1
      else
*   ..."^" type reflection
        trans = interface(ipath,iseg)
      endif

*   loading in proper parameters

      iinc = (inc-1)*9
      iscat = (scat-1)*9
      itrans = (trans-1)*9


      do10 i=1,6
        elasti(i) = param(iinc+i)
        elasts(i) = param(iscat+i)
        elastt(i) = param(itrans+i)
10      continue

      do20 i=1,3
        axisi(i) = param(iinc+i+6)
        axiss(i) = param(iscat+i+6)
        axist(i) = param(itrans+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
*   ...downward transmission
        normal(3) = -1.
        which(2) = 2
      else if(inc .GT. scat) then
*   ...upward transmission
        normal(3) = 1.
        which(2) = 2
      else if(interface(ipath,iseg-1) .LT. interface(ipath,iseg)) then
*   ..."v" type reflection
        normal(3) = -1.
        which(2) = 1
      else
*   ..."^" type reflection
        normal(3) = 1.
        which(2) = 1
      endif

*   compute the five scattered slownesses


      call scatslow(p0,normal,axiss,which,elasts,p,isol)
      call scatslow(p10,normal,axiss,which,elasts,p1,isol1)
      call scatslow(p20,normal,axiss,which,elasts,p2,isol2)
      call scatslow(p30,normal,axiss,which,elasts,p3,isol3)
      call scatslow(p40,normal,axiss,which,elasts,p4,isol4)

*   test for post-critical scattering and interrupts tracing if needed

      if((isol+isol1+isol2+isol3+isol4) .NE. 0) then
        error=1
        return
      endif

*   compute the scattered ray jacobian from its incident value


      call groupvel(p0,axisi,imode,elasti,w0)
      call groupvel(p,axiss,smode,elasts,w)

      dot0 = 0.
      dot = 0.

      do30 i=1,3
        dot0 = dot0+w0(i)*normal(i)
        dot = dot+w(i)*normal(i)
30      continue

      if(dot0.EQ.0.) then
*   ...caustic: should not happen, hence a fatal error
        error = 2
        return
      endif

      j = j0*dot/dot0
      j = abs(j)


*   compute the scattering coefficients

      call czoeppritz(p0,imode,normal,axisi,elasti,
     &                axist,elastt,coefs,error)


*   exits if any error or unsuccessful computation occured

      if(error .NE. 0) then
        error = 3
        return
      endif

      if((imode.EQ.1) .AND. (triso(inc).EQ.0)) then

*   if incident medium is isotropic and incident mode is S, then 
*   combines the SP and QS contributions together. SP scattering 
*   coefficients have already been computed. Now compute QS scattering 
*   coefficients (array coefs2(6)) and combine.

        i = 3
        call czoeppritz(p0,i,normal,axisi,elasti,axist,elastt,
     &          coefs2,error)

        if(error .NE. 0) then
                error = 3
                return
        endif

        if(inc .EQ. scat) then
                amp1 = coefs(smode)*amp10 + coefs2(smode)*amp20
                amp2 = coefs(3)*amp10 + coefs2(3)*amp20
        else
                amp1 = coefs(smode+3)*amp10 + coefs2(smode+3)*amp20
                amp2 = coefs(6)*amp10 + coefs2(6)*amp20
        endif
                
      else if(inc .EQ. scat) then

*   computes scattering amplitude for the appropriate mode, 
*   and saves QS scattering amplitude in case scattering layer 
*   is isotropic.

        amp1 = coefs(smode)*amp10
        amp2 = coefs(3)*amp10
      else
        amp1 = coefs(smode+3)*amp10
        amp2 = coefs(6)*amp10
      endif

      return
      end



      subroutine cendray(amp1,amp2,p0,mode,disp,error)

      complex amp1, amp2, disp(3)
      real p0(3)
      integer mode, error

*
*   computes the total displacement at the free elastic surface
*   of a transversely isotropic medium with arbitrary axis of symmetry. 
*
*   INPUT VARIABLES
*
*       p0(3)           incident slowness vector
*       amp1            amplitude of incident mode
*       amp2            amplitude of incident QS wave
*       mode            incident wave type (1 for S, 2 for P)
*
*   OUTPUT VARIABLES
*
*       disp(3)         total displacement vector at surface
*       error           0 if no error, 1 if wrong mode,
*       
*   GLOBAL VARIABLES
*
*       block "model"   elastic parameters
*

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

      real axis(3), elast(6), normal(3)
      complex ui(3), urSP(3), urQP(3), urQS(3)
      complex coefs(3)
      integer i

      do5 i=1,6
         elast(i) = param(i)
5     continue

      do10 i=1,3
         axis(i) = param(6+i)
10    continue

      normal(1) = 0.e0
      normal(2) = 0.e0
      normal(3) = 1.e0

*     computes free surface scattering parameters for incident mode

      call freezoep(p0,mode,normal,axis,elast,coefs,error
     &              ,ui,urSP,urQP,urQS)

      if(error .NE. 0) then
         error = 1
         return
      endif         

*     assembles total free-surface displacement as sum of 
*     incident and reflected displacements

      do20 i=1,3
         disp(i) = amp1*(ui(i)+coefs(1)*urSP(i)
     &                    +coefs(2)*urQP(i)+coefs(3)*urQS(i))
20    continue

*     if top layer is isotropic and incident mode is S, then must 
*     add QS contribution (SP just computed above)

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

         call freezoep(p0,3,normal,axis,elast,coefs,error
     &                 ,ui,urSP,urQP,urQS)
         if(error .NE. 0) then
            error=1
            return
         endif

         do30 i=1,3
            disp(i) = disp(i)+amp2*(ui(i)+coefs(1)*urSP(i)
     &              + coefs(2)*urQP(i)+coefs(3)*urQS(i))
30       continue

      endif

      return
      end



      subroutine cgausselim(mat,vec,n)
      integer n
      complex mat(n,n), vec(n)

*
*   complex gaussian elimination routine with partial pivoting
*   and line index permutation for systems of degree less than 100.
*   
*   INPUT VARIABLES:
*
*       mat(n,n)        matrix of linear system
*       vec(n)          forcing vector
*       n               dimension of system
*
*   OUTPUT VARIABLES:
*       mat(n,n)        reduced matrix (useless)
*       vec(n)          solution of system
*
*   NOTE... 
*   ... this routine will modify the contents of both 
*   mat and vec. Moreover, mat should be declared statically
*   in the calling routine EXACTLY with the same dimension 
*   that is used when this routine is called.
*

      complex pivot, try, res(100), coef
      integer index(100), i, j, k, l, ind, jpiv

      if(n .GT. 100) then
        return
      endif

*   INITIALIZATION STEP

      do 10 i=1,n
        index(i) = i
10      continue

*   ELIMINATION LOOP

      do 20 i=1,n
        ind = index(i)
        jpiv = i
        pivot = mat(ind,i)

*       FINDING PIVOT

        do 30 j=i+1,n
                k = index(j)
                try = mat(k,i)
                if(cabs(try) .GT. cabs(pivot)) then
                        pivot = try
                        ind = k
                        jpiv = j
                endif
30         continue

*       TESTING FOR SINGULARITY AND PERMUTING INDEXES

        if(cabs(pivot) .EQ. 0.) then
                return
        else
                index(jpiv) = index(i)
                index(i) = ind
        endif

*       REACTUALIZING mat AND vec  (ELIMINATION) 

        do 40 j=i+1,n
                k = index(j)
                vec(k) = vec(k) - mat(k,i)*vec(ind)/pivot
                coef = mat(k,i)/pivot
                do 50 l=i+1,n
                        mat(k,l) = mat(k,l) - mat(ind,l)*coef
50                 continue
*               mat(k,i) = 0.
40         continue
20      continue                   

*   BACKSUBSTITUTION

      do 60 i=n,1,-1
        k = index(i)
        res(i) = 0.e0
        do 70 j=i+1,n
                res(i) = res(i) + res(j)*mat(k,j)
70         continue
        res(i) = (vec(k)-res(i))/mat(k,i)
60      continue

*   PUTS RESULT BACK IN vec

      do 80 i=1,n
        vec(i) = res(i)
80      continue

      return
      end



      subroutine cinitray(razi,rdip,mode,rg1,rg2,r,rdg,x0,x10,x20,
     &   x30,x40,p0,p10,p20,p30,p40,t0,t10,t20,t30,t40,j0,amp10,amp20)

      real razi, rdip, rg1, rg2, r, rdg, x0(3), x10(3), x20(3), x30(3)
      real x40(3), p0(3), p10(3), p20(3), p30(3), p40(3), t0, t10
      real t20, t30, t40, j0
      complex amp10, amp20
      integer mode

*
*   computes the initial parameters for a ray emanating from 
*   a point source at the free surface of a homogeneous isotropic
*   medium. The elastic parameters are in param.
*
*   INPUT VARIABLES:
*
*       razi    azimuth of point force in degrees (angle from 
*               geographic North in a horizontal plane)
*       rdip    dip of force in degrees (inclination of force
*               with respect to horizontal)
*       mode    type of excitation: 1 for S, 2 for P
*       rg1     take-off azimuth angle for the ray in degrees
*       rg2     take-off dip angle for the ray in degrees
*       r       radius of take-off sphere
*       rdg     spherical angle defining the angular aperture
*               of the ray tube in degrees
*
*   OUTPUT VARIABLES:
*
*       x0(3)   vector location of ray starting point with respect
*               to source location              
*       x10(3)  vector location of first tube ray starting point 
*               with respect to source location         
*       x20(3)  vector location of second tube ray starting point 
*               with respect to source location         
*       x30(3)  vector location of third second tube ray starting point 
*               with respect to source location         
*       x40(3)  vector location of fourth tube ray starting point 
*               with respect to source location         
*       p0(3)   initial ray slowness
*       p10(3)  initial slowness for first tube ray
*       p20(3)  initial slowness for second tube ray
*       p30(3)  initial slowness for third first tube ray
*       p40(3)  initial slowness for fourth tube ray
*       t0      initial ray traveltime
*       t10     initial traveltime for first tube ray
*       t20     initial traveltime for second tube ray