
      subroutine m0dcf(mtrep,mf,mt,d,v,miso,mdc,mclvd,m0,mw,axis,
     1           plunge,azimuth,strike,dip,slip,pcdc,pcclvd,pciso)
c
c ...... decompose moment tensor into various representations
c
c     input parameters -
c
c         mtrep - moment tensor representation
c                 = 0 for spherical coordinates (r,d,p)
c                   ( Up, Away, South )
c                 = 1 for cartesian coordinates (x,y,z)
c                   ( North, East, Down )
c                 = 2 for f1,...,f6 notation
c
c         mf or mt - components of moment tensor
c                 | mrr mrd mrp |    | mxx mxy mxz |    | f1 f2 f3 |
c                 | mdr mdd mdp | or | myx myy myz | or | f2 f4 f5 |
c                 | mpr mpd mpp |    | mzx mzy mzz |    | f3 f5 f6 |
c
c              NOTE:
c                   f1 = mrr =  mzz
c                   f2 = mdd =  mxx
c                   f3 = mpp =  myy
c                   f4 = mrd =  mxz
c                   f5 = mrp = -myz
c                   f6 = mdp = -mxy
c     
c     output parameters -
c
c         d       - eigenvalues
c         v       - principal eigenvectors
c         miso    - isotropic moment tensor
c         mdc     - double-couple moment tensor
c         mclvd   - compensated linear vector dipole moment tensor
c         m0      - scalar seismic moment
c         axis    - component axis ( P, N, T )
c         plunge  - principle eigenvector plunge 
c                   (degrees below horizontal)
c         azimuth - principle eigenvector azimuth
c                   (degrees clockwise from North)
c         strike  - strike of fault plane
c                   (degrees clockwise from North)
c         dip     - dip of fault plane
c                   (degrees below horizontal)
c         slip    - slip of fault plane
c                   (degrees from horizontal)
c         pciso   - percent isotropic source
c         pcdc    - percent double-couple source
c         pcclvd  - percent CLVD source
c  
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      character*1 axis(3)
      real*8 mt(3,3),d(3),v(3,3),miso(3,3),mdc(3,3),mclvd(3,3),m0,mw
      real*8 plunge(3),azimuth(3),strike(3),dip(3),slip(3)
      real*8 pcdc,pcclvd,a(3,3),p(3),t(3),mf(6)
      n=3
      np=3
      rad2deg=180.0/3.1415
      deg2rad=1./rad2deg
c
      call mtsrce(mtrep,mf,mt,a)
c     do 666 j=1,3
c     write(*,*) (a(j,i),i=1,3)
c666  continue
      call mdcmp(a,n,np,v,d,miso,mdc,mclvd)
      call pcnts(d,pcdc,pcclvd)
      call plaz(d,v,plunge,azimuth,axis)
      call m0find(d,m0)
c
c ...... calculate percentage of DC, CLVD & Iso in source
c
      pcdc=pcdc/100.d0
      pcclvd=pcclvd/100.d0
      pciso=dabs(miso(1,1))/m0
      pcsum=pcdc+pcclvd+pciso
      pcdc=100.d0*pcdc/pcsum
      pcclvd=100.d0*pcclvd/pcsum
      pciso=100.d0*pciso/pcsum
c
c ...... compute moment magnitude
c
      mw=dlog10(m0)/1.5-10.7
c
      call ptfind(v,axis,p,t)
      call fpsol(p,t,strike,dip,slip) 
c
c      write(*,*) ' '
c      write(*,'(a,1p3e12.2)') 'eigenvalues = ',(d(i),i=1,n)
c      write(*,*) ' '
c      write(*,'(a,1p3e12.2)') 'eigenvector = ',(v(1,i),i=1,n)
c      write(*,'(a,1p3e12.2)') 'eigenvector = ',(v(2,i),i=1,n)
c      write(*,'(a,1p3e12.2)') 'eigenvector = ',(v(3,i),i=1,n)
c      write(*,*) ' '
c      write(*,'(a,1p3e12.2)') 'miso        = ',(miso(1,i),i=1,n)
c      write(*,'(a,1p3e12.2)') 'miso        = ',(miso(2,i),i=1,n)
c      write(*,'(a,1p3e12.2)') 'miso        = ',(miso(3,i),i=1,n)
c      write(*,*) ' '
c      write(*,'(a,1p3e12.2)') 'mdc         = ',(mdc(1,i),i=1,n)
c      write(*,'(a,1p3e12.2)') 'mdc         = ',(mdc(2,i),i=1,n)
c      write(*,'(a,1p3e12.2)') 'mdc         = ',(mdc(3,i),i=1,n)
c      write(*,*) ' '
c      write(*,'(a,1p3e12.2)') 'mclvd       = ',(mclvd(1,i),i=1,n)
c      write(*,'(a,1p3e12.2)') 'mclvd       = ',(mclvd(2,i),i=1,n)
c      write(*,'(a,1p3e12.2)') 'mclvd       = ',(mclvd(3,i),i=1,n)
c      write(*,*) ' '
c      write(*,*) ' '
c      write(*,'(a,1pe10.2)') 'Scalar Moment    = ',m0
c      write(*,*) ' '
c      write(*,'(a,  f10.2)') 'Moment Magnitude = ',mw
c      write(*,*) ' '
c      write(*,'(a,7x,a1,2(9x,a1))') 'Axis             = ',
c     1         (axis(i),i=1,n)
c      write(*,*) ' '
c      write(*,'(a,3f10.2)') 'Plunge           = ',(plunge(i),i=1,n)
c      write(*,*) ' '
c      write(*,'(a,3f10.2)') 'Azimuth          = ',(azimuth(i),i=1,n)
c      write(*,*) ' '
c      write(*,*) ' '
c      write(*,'(2a)')       'Plane            = ',
c     1           '      P1        P2'
c      write(*,*) ' '
c      write(*,'(a,2f10.2)') 'Strike           = ',(strike(i),i=1,n-1)
c      write(*,*) ' '
c      write(*,'(a,2f10.2)') 'Dip              = ',(dip(i),i=1,n-1)
c      write(*,*) ' '
c      write(*,'(a,2f10.2)') 'Slip             = ',(slip(i),i=1,n-1)
c      write(*,*) ' '
c      write(*,*) ' '
c      write(*,*)            'Source Composition     Percent'
c      write(*,*) ' '
c      write(*,'(a, f10.2)') 'Isotropic        = ',pciso
c      write(*,*) ' '
c      write(*,'(a, f10.2)') 'Double-Couple    = ',pcdc
c      write(*,*) ' '
c      write(*,'(a, f10.2)') 'CLVD             = ',pcclvd
c      write(*,*) ' '
c
      return
      end

      subroutine pcnts(d,pcdc,pcclvd)
c
c ...... calculate percentages of double-couple, CLVD & iso
c        in seismic source
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      real*8 d(3),dd(3)
      do 1 i=1,3
      dd(i)=dabs(d(i))
    1 continue
      do 3 i=1,3
      do 2 j=i,3
      if(dd(j).lt.dd(i)) then
        temp=dd(i)
        dd(i)=dd(j)
        dd(j)=temp
      endif
    2 continue
    3 continue
      eps=dd(1)/dd(3)
      pcdc=100.d0*(1.d0-2.d0*eps)
      pcclvd=200.d0*eps
      return
      end

      subroutine m0find(d,m0)
c
c ...... find scalar moment M0 
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      real*8 d(3),m0
      emin=d(1)
      emax=d(1)
      do 1 i=2,3
      if(d(i).lt.emin) emin=d(i)
      if(d(i).gt.emax) emax=d(i)
    1 continue
      m0=(dabs(emin)+dabs(emax))/2.d0
      return
      end

      subroutine ptfind(v,axis,p,t)
c
c ...... find P & T vectors
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      real*8 v(3,3),p(3),t(3)
      character*1 axis(3),pch,tch
      pch='P'
      tch='T'
      do 3 i=1,3
      if(axis(i).eq.pch) then
        do 1 j=1,3
        p(j)=v(j,i)
    1   continue
      elseif(axis(i).eq.tch) then
        do 2 j=1,3
        t(j)=v(j,i)
    2   continue
      endif
    3 continue
      return
      end

      subroutine mtsrce(mtrep,mf,mt,a)
c
c ...... set up appropriate moment tensor components from input
c        moment tensor representation
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      real*8 mf(6),mt(3,3),a(3,3)
      if(mtrep.eq.0) then
        a(1,1)= mt(2,2)
        a(1,2)=-mt(2,3)
        a(1,3)= mt(2,1)
        a(2,1)=-mt(2,3)
        a(2,2)= mt(3,3)
        a(2,3)=-mt(3,1)
        a(3,1)= mt(2,1)
        a(3,2)=-mt(3,1)
        a(3,3)= mt(1,1)
      elseif(mtrep.eq.1) then
        a(1,1)=mt(1,1)
        a(1,2)=mt(1,2)
        a(1,3)=mt(1,3)
        a(2,1)=mt(1,2)
        a(2,2)=mt(2,2)
        a(2,3)=mt(2,3)
        a(3,1)=mt(1,3)
        a(3,2)=mt(2,3)
        a(3,3)=mt(3,3)
      else
        a(1,1)= mf(2)
        a(1,2)=-mf(6)
        a(1,3)= mf(4)
        a(2,1)=-mf(6)
        a(2,2)= mf(3)
        a(2,3)=-mf(5)
        a(3,1)= mf(4)
        a(3,2)=-mf(5)
        a(3,3)= mf(1)
      endif
      return
      end
 
      subroutine fpsol(p,t,strike,dip,slip)
c
c ...... calculate strike, dip & slip of fault planes
c        from the P & T vectors
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      real*8 p(3),t(3),strike(2),dip(2),slip(2)
      real*8 u(2,3),nu(2,3),lambda
      real*8 rad2deg, deg2rad
      integer*4 idarg(2)

      rad2deg=180.0/3.1415
      deg2rad=1./rad2deg

      con=1.d0/dsqrt(2.d0)
      do 1 i=1,3
      u(1,i)=con*(t(i)+p(i))
      nu(1,i)=con*(t(i)-p(i))
      u(2,i)=con*(t(i)-p(i))
      nu(2,i)=con*(t(i)+p(i))
    1 continue
c      write(*,11) (p(i),i=1,3)
c   11 format(1h ,'P-vector = ',3f10.6)
c      write(*,12) (t(i),i=1,3)
c   12 format(1h ,'T-vector = ',3f10.6)
      idarg(1)=0
      idarg(2)=0      
      do 2 i=1,2
      dip(i)=rad2deg*dacos(-nu(i,3))
      if((nu(i,1).eq.0.d0).and.(nu(i,2).eq.0.d0)) then
        strike(i)=0.d0
      else
        strike(i)=rad2deg*datan2(-nu(i,1),nu(i,2))
      endif
    2 continue
      do 3 i=1,2
      sstr=dsin(deg2rad*strike(i))
      cstr=dcos(deg2rad*strike(i))
      sdip=dsin(deg2rad*dip(i))
      cdip=dcos(deg2rad*dip(i))
      if(dabs(sdip).gt.0.d0) then
        lambda=rad2deg*dasin(-u(i,3)/dsin(deg2rad*dip(i)))
      else
        arg1=1.d0
        arg2=u(i,3)
        arg=dsign(arg1,arg2)
        if(arg.lt.0.d0) then
          lambda=180.d0
        else
          lambda=0.d0
        endif
      endif
      slambda=dsin(deg2rad*lambda)
      cdsl=cdip*slambda
      if(dabs(sstr).gt.dabs(cstr)) then
        clambda=(u(i,2)+cdsl*cstr)/sstr
      else
        clambda=(u(i,1)-cdsl*sstr)/cstr
      endif
      if((slambda.eq.0.d0).and.(clambda.eq.0.d0)) then
        slip(i)=0.d0
      else
        slip(i)=rad2deg*datan2(slambda,clambda)
      endif
      if(dip(i).gt.90.d0) then
        dip(i)=180.d0-dip(i)
        strike(i)=strike(i)+180.d0
        slip(i)=360.d0-slip(i)
      endif
      if(strike(i).lt.0.d0) strike(i)=strike(i)+360.d0
      if(slip(i).ge.180.d0) slip(i)=slip(i)-360.d0
    3 continue
      return
      end

      subroutine plaz(d,v,plunge,azimuth,axis)
c
c ...... calculate plunge & azimuth of eigenvectors
c        in degrees
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      real*8 d(3),v(3,3),plunge(3),azimuth(3)
      character*1 axis(3)
      do 1 i=1,3
      if(v(3,i).lt.0.) then
	v(1,i)=-v(1,i)
	v(2,i)=-v(2,i)
	v(3,i)=-v(3,i)
      endif
    1 continue
      do 2 i=1,3
      if((v(2,i).eq.0.d0).and.(v(1,i).eq.0.d0)) then
        azimuth(i)=0.d0
      else
        azimuth(i)=datan2(v(2,i),v(1,i))
      endif
      if(azimuth(i).lt.0.) azimuth(i)=azimuth(i)+360.0
      r=dsqrt(v(1,i)*v(1,i)+v(2,i)*v(2,i))
      if((v(3,i).eq.0.d0).and.(r.eq.0.d0)) then
        plunge(i)=0.d0
      else
        plunge(i)=datan2(v(3,i),r)
      endif
    2 continue
      axis(1)='N'
      if(d(2).gt.d(3)) then
        axis(2)='T'
        axis(3)='P'
      else
        axis(2)='P'
        axis(3)='T'
      endif
      return
      end

      subroutine mdcmp(mij,n,np,v,d,miso,mdc,mclvd)
c
c ...... calculate components of an arbitrary moment tensor
c
      implicit real*8 (a-h,m,o-z)
      implicit integer*4 (i-l,n)
      real*8 mij(np,np),miso(np,np),mdc(np,np),mclvd(np,np)
      real*8 a(3,3),v(3,3),d(3)
      real*8 a1a1(3,3),a2a2(3,3),a3a3(3,3)
      do 2 i=1,n
      do 1 j=1,n
      a(i,j)=mij(i,j)
    1 continue
    2 continue
      trace=a(1,1)+a(2,2)+a(3,3)
      tr3=trace/3.0d0
      do 3 i=1,n
      a(i,i)=a(i,i)-tr3
    3 continue
      do 5 i=1,n
      do 4 j=1,n
      if(i.eq.j) then
        miso(i,j)=tr3
      else
        miso(i,j)=0.d0
      endif
    4 continue
    5 continue
      call jacobi(a,n,np,d,v,nrot)
c      write(*,*) 'jacobi - nrot = ',nrot
      call eigsrt(d,v,n,np)
      f=-d(1)/d(3)
      c=d(3)*(1.d0-2.d0*f)
      call dyadic(v,3,3,c,n,np,a3a3)
      c=-c
      call dyadic(v,2,2,c,n,np,a2a2)
      call matsum(a3a3,a2a2,n,np,mdc)
      c=2.d0*d(3)*f
      call dyadic(v,3,3,c,n,np,a3a3)
      c=-d(3)*f
      call dyadic(v,2,2,c,n,np,a2a2)
      call dyadic(v,1,1,c,n,np,a1a1)
      call matsum(a3a3,a2a2,n,np,mclvd)
      call matsum(mclvd,a1a1,n,np,mclvd)
      return
      end

      subroutine dyadic(v,n1,n2,c,n,np,d)
c
c ...... calculate dyadic of eigenvectors v(i,n1)*v(j,n2)
c        scaled by c
c
      implicit real*8 (a-h,o-z) 
      implicit integer*4 (i-n)
      real*8 v(np,np),d(np,np) 
      do 2 i=1,n
      do 1 j=1,n
      d(i,j)=v(i,n1)*v(j,n2)*c
    1 continue
    2 continue
      return
      end

      subroutine matsum(a,b,n,np,c)
c
c ...... calculate matrix sum c = a + b
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      real*8 a(np,np),b(np,np),c(np,np)
      do 2 i=1,n
      do 1 j=1,n
      c(i,j)=a(i,j)+b(i,j)
    1 continue 
    2 continue
      return
      end

      subroutine rotate(v,n,np,a,b)
c
c ...... rotate matrix
c                     T
c        b = v * a * v
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      real*8 v(np,np),a(np,np),b(np,np),c(3,3)
      do 3 i=1,n
      do 2 j=1,n
      sum=0.d0
      do 1 k=1,n
      sum=sum+a(i,k)*v(j,k)
    1 continue
      c(i,j)=sum
    2 continue
    3 continue
      do 6 i=1,n
      do 5 j=1,n
      sum=0.d0
      do 4 k=1,n
      sum=sum+v(i,k)*c(k,j)
    4 continue
      b(i,j)=sum
    5 continue
    6 continue
      return
      end

      subroutine jacobi(a,n,np,d,v,nrot)
c
c ...... Computes all eigenvalues and eigenvectors of a real 
c        symmetric matrix A, which is of size N by N, stored
c        in a physical np by np array.  On output, elements
c        of A above the diagonal are destroyed.  D returns
c        the eigenvalues of A in its first N elements.  V is
c        a matrix with the same logical and physical dimensions
c        as A whose columns contain, on output, the normalized
c        eigenvectors of A.  NROT returns the number of Jacobi
c        rotations which were required.
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      parameter (nmax=100)
      real*8 a(np,np),d(np),v(np,np),b(nmax),z(nmax)
      do 12 ip=1,n
        do 11 iq=1,n
          v(ip,iq)=0.d0
   11   continue
        v(ip,ip)=1.d0
   12 continue
      do 13 ip=1,n
        b(ip)=a(ip,ip)
        d(ip)=b(ip)
        z(ip)=0.d0
   13 continue
      nrot=0
      do 24 i=1,500
        sm=0.d0
        do 15 ip=1,n-1
          do 14 iq=ip+1,n
            sm=sm+dabs(a(ip,iq))
   14     continue
   15   continue
        if(sm.eq.0.d0) return
        if(i.lt.1) then
          thresh=0.2d0*sm/dble(n**2)
        else
          thresh=0.d0
        endif
        do 22 ip=1,n-1
          do 21 iq=ip+1,n
            g=100.d0*dabs(a(ip,iq))
            if((i.gt.4).and.(dabs(d(ip))+g.eq.dabs(d(ip)))
     1        .and.(dabs(d(iq))+g.eq.dabs(d(iq)))) then
              a(ip,iq)=0.d0
            elseif(dabs(a(ip,iq)).gt.thresh) then
              h=d(iq)-d(ip)
              if(dabs(h)+g.eq.dabs(h)) then
                t=a(ip,iq)/h
              else
                theta=0.5d0*h/a(ip,iq)
                t=1.d0/(dabs(theta)+dsqrt(1.d0+theta**2))
                if(theta.lt.0.d0) t=-t
              endif
              c=1.d0/dsqrt(1.d0+t**2)
              s=t*c
              tau=s/(1.d0+c)
              h=t*a(ip,iq)
              z(ip)=z(ip)-h
              z(iq)=z(iq)+h
              d(ip)=d(ip)-h
              d(iq)=d(iq)+h
              a(ip,iq)=0.d0
              do 16 j=1,ip-1
                g=a(j,ip)
                h=a(j,iq)
                a(j,ip)=g-s*(h+g*tau)
                a(j,iq)=h+s*(g-h*tau)
   16         continue
              do 17 j=ip+1,iq-1
                g=a(ip,j)
                h=a(j,iq)
                a(ip,j)=g-s*(h+g*tau)
                a(j,iq)=h+s*(g-h*tau)
   17         continue
              do 18 j=iq+1,n
                g=a(ip,j)
                h=a(iq,j)
                a(ip,j)=g-s*(h+g*tau)
                a(iq,j)=h+s*(g-h*tau)
   18         continue
              do 19 j=1,n
                g=v(j,ip)
                h=v(j,iq)
                v(j,ip)=g-s*(h+g*tau)
                v(j,iq)=h+s*(g-h*tau)
   19         continue
            nrot=nrot+1
            endif
   21     continue
   22   continue
        do 23 ip=1,n
          b(ip)=b(ip)+z(ip)
          d(ip)=b(ip)
          z(ip)=0.d0
   23   continue
   24 continue
      return
      end
      
      subroutine eigsrt(d,v,n,np)
c
c ...... Given the eigenvalues D and eigenvectors V as output from
c        JACOBI, this routine sorts the eigenvalues into increasing
c        absolute order and rearranges the columns of V 
c        correspondingly.  The mothod is straight insertion.
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      real*8 d(np),v(np,np)
      do 13 i=1,n-1
        k=i
        p=d(i)
        do 11 j=i+1,n
          if(dabs(d(j)).lt.dabs(p)) then
            k=j
            p=d(j)
          endif
   11   continue
        if(k.ne.i) then
          d(k)=d(i)
          d(i)=p
          do 12 j=1,n
            p=v(j,i)
            v(j,i)=v(j,k)
            v(j,k)=p
   12     continue
        endif
   13 continue
      return
      end
