c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c+                                                                +
c+              4-node plate bending element program              +
c+                                           by T. Nishihara      +
c+                                                                +
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      program kani
      implicit real*8 (a-h,o-z)
c
      parameter( nele = 5, node = 12, ndof = 6 * node )
      parameter( nbcc = 2, nout = 1 )
c
      dimension  bc(ndof), forc(ndof),dd(ndof+1), d(ndof,24)
      dimension  s(24,24), sk(24,24), sinc(ndof+1,ndof+1)
      dimension  ddl(ndof), xb(ndof+1), d24(24), dd24(24)
      dimension  ba24(24), ba(ndof+1), alp3(3), cor(ndof,2)
      dimension  alp(ndof), snr2(ndof,ndof)
      dimension  jcon(4), nodn(ndof,4), mds(nout), mmds(nout)
      dimension  nbc(nbcc)
c
      character*20 file1
c
      open(7,file='input.dat')
      open(8,file='outkani.dat')
c
c********************************************************************
c     データの読み込み
c********************************************************************
c
      call clrv (ndof,alp)
c
      read(7,*) th,e,um
c
      do 11 i=1,node
         read(7,*) cor(i,1), cor(i,2)
   11 continue
c
      do 13 i=1,nele
         read(7,*) nodn(i,1), nodn(i,2), nodn(i,3), nodn(i,4)
   13 continue
c
      do 50 i=1,ndof
         bc(i) = 1.d0
   50 continue
c
      read(7,*) nxx, nyy, nzz, nrx, nry
      read(7,*) (nbc(i), i=1, nxx)
      do 130 i=1,nxx
         bc ( nbc(i) * 6 - 5 ) = 0.d0
  130 continue
      read(7,*) (nbc(i), i=1, nyy)
      do 131 i=1,nyy
         bc ( nbc(i) * 6 - 4 ) = 0.d0
  131 continue
      read(7,*) (nbc(i), i=1, nzz)
      do 132 i=1,nzz
         bc ( nbc(i) * 6 - 3 ) = 0.d0
  132 continue
      read(7,*) (nbc(i), i=1, nrx)
      do 133 i=1,nrx
         bc ( nbc(i) * 6 - 2 ) = 0.d0
  133 continue
      read(7,*) (nbc(i), i=1, nry)
      do 134 i=1,nry
         bc ( nbc(i) * 6 - 1 ) = 0.d0
  134 continue
      do 135 i=1,node
         bc ( i * 6 ) = 0.d0
  135 continue
c
      call clrv (ndof,forc)
      read(7,*) nf
      do 60 i=1,nf
         read(7,*) nodd, xfor, yfor, zfor, xrot, yrot, zrot
         forc( nodd * 6 - 5 ) = xfor
         forc( nodd * 6 - 4 ) = yfor
         forc( nodd * 6 - 3 ) = zfor
         forc( nodd * 6 - 2 ) = xrot
         forc( nodd * 6 - 1 ) = yrot
         forc( nodd * 6 ) = zrot
   60 continue
c
      read(7,*) np
      do 560 i=1,np
         read(7,561) mds(i), mmds(i), file1
         open(20+i,file=file1)
  560 continue
c 後藤追記（04/1/19）
c Linux 上の g77 では大丈夫なようだけど、Cygwin の g77 とかでは、
c 以下のフォーマットを入れないとちゃんと読めないようだ。
c だから入力データもここのフォーマットと合わせて書く必要がある
561   format(2i3, a6)
c
       read(7,*) ma
       if ( ma.eq.0 ) go to 695
      do 570 ie=1,nele
           read(7,*) alp(ie*3-2),alp(ie*3-1),alp(ie*3)
  570 continue
c
  695 continue
c
      write(8,*)'===================================================='
      write(8,*)'                output of program'
      write(8,*)'===================================================='
      write(8,*)
      write(8,*)
      write(8,*)'----------------------------------------------------'
      write(8,*)'                     parameter'
      write(8,*)'----------------------------------------------------'
      write(8,2005) nele, node
      write(8,2001) e, um
      write(8,2002) th
c
 2001 format('Young modulas = ',1pe15.7,5x,'Poisson ration = ',
     &       1pe15.7)
 2002 format('Thickness = ',1pe15.7)
 2005 format('Number of element:',i5,5x,'Number of node:',i5)
c
      write(8,*)
      write(8,*)'----------------------------------------------------'
      write(8,*)'               co-ordinate of nodes'
      write(8,*)'----------------------------------------------------'
      write(8,2014)
      do 12 i=1,node
         write(8,2015) i,cor(i,1), cor(i,2)
   12 continue
 2014 format(2x,'node number',12x,'x',20x,'y')
 2015 format('node = ',i5,5x,1pe15.7,5x,1pe15.7)
c
      write(8,*)
      write(8,*)'----------------------------------------------------'
      write(8,*)'                       element'
      write(8,*)'----------------------------------------------------'
      do 40 i=1,nele
         write(8,2007) i,nodn(i,1),nodn(i,2),nodn(i,3),nodn(i,4)
   40 continue
 2007 format('element = ',i5,5x,i5,i5,i5,i5)
c
      write(8,*)
      write(8,*)'----------------------------------------------------'
      write(8,*)'                 boundary condition'
      write(8,*)'----------------------------------------------------'
      write(8,2009)
      do 2008 i=1,node
         write(8,2010) i,bc(6*i-5),bc(6*i-4),bc(6*i-3),
     &                bc(6*i-2),bc(6*i-1),bc(6*i)
 2008 continue
 2009 format(3x,'node number',8x,'x',9x,'y',9x,'z',7x,
     &       'rot-x',5x,'rot-y',5x,'rot-z')
 2010 format('node = ',i5,1x,
     &       1pe10.1,1pe10.1,1pe10.1,1pe10.1,1pe10.1,1pe10.1)
c
      write(8,*)
      write(8,*)'----------------------------------------------------'
      write(8,*)'                    nodal force'
      write(8,*)'----------------------------------------------------'
      write(8,2009)
      do 2011 i=1,node
         write(8,2010) i,forc(6*i-5),forc(6*i-4),forc(6*i-3),
     &                forc(6*i-2),forc(6*i-1),forc(6*i)
 2011 continue
c
       if ( ma.eq.0 ) go to 696
      write(8,*)
      write(8,*)'----------------------------------------------------'
      write(8,*)'            initial angle of element'
      write(8,*)'----------------------------------------------------'
      write(8,2017)
      do 571 ie=1,nele
           write(8,2016) ie,alp(ie*3-2),alp(ie*3-1),alp(ie*3)
  571 continue
 2017 format(2x,'element number',10x,'alpha',15x,'beta',15x,'gamma')
 2016 format('element = ',i5,5x,1pe15.7,5x,1pe15.7,5x,1pe15.7)
  696 continue
c
      write(*,*)'start loding step'
c
c********************************************************************
c     第1増分ステップの第1収斂ステップ
c********************************************************************
      read(7,*) r
      write(8,*)
      write(8,*) 
      write(8,2012) r
      write(8,*)
      read(7,*) rn
      write(8,2013) rn
 2012 format('first arc = ',1pe15.7)
 2013 format('next arc  = ',1pe15.7)
c
      do 225 i=1,ndof
         do 225 j=1,24
            d(i,j) = 0.d0
  225 continue
      call clrv (ndof+1,dd)
      call clrv (ndof,ddl)
      call clrv (ndof+1,ba)
      dp = 0.d0
      p = 0.d0
      dpl = 0.d0
      r0 = 0.d0
      write(3) dd,d,ddl,dp,dpl,r0,p
      rewind 3
      r0 = r
      dpl = r
      nacl = 0
      nbcl = 0
      nccl = 0
      konv = 0
c
c********************************************************************
c     弧長増分過程
c********************************************************************
c
 1000 continue
      nofcs = 0
      ninc = ninc + 1
c      if(ninc.gt.1000) then
c         write(12,*)
c         write(12,*)'収束しないよ！'
c         close(12)
c         close(13)
c         stop
c      endif
c
      call clrmx (ndof+1,sinc)
c
      do 2000 ie=1,nele
      do 440 i = 1,4
         jcon(i) = nodn(ie, i)
  440 continue
         ib1 = (jcon(1) - 1) * 6
         ib2 = (jcon(2) - 1) * 6
         ib3 = (jcon(3) - 1) * 6
         ib4 = (jcon(4) - 1) * 6
c
      a = cor( jcon(2), 1 ) - cor( jcon(1), 1 )
      b = cor( jcon(4), 2 ) - cor( jcon(1), 2 )
      call smx(s,th,e,um,a,b)
c
      do 800 i=1,3
         alp3(i) = alp((ie-1)*3+i)
  800 continue
c
      do 150 i=1,6
         dd24(i) = dd(ib1+i)
         dd24(i+6) = dd(ib2+i)
         dd24(i+12) = dd(ib3+i)
         dd24(i+18) = dd(ib4+i)
         d24(i) = d(ie,i)
         d24(i+6) = d(ie,i+6)
         d24(i+12) = d(ie,i+12)
         d24(i+18) = d(ie,i+18)
  150 continue
c
      call stmx (s,sk,a,b,d24,dd24,ba24,alp3)
c
      do 160 i=1,6
         d(ie,i) = d24(i)
         d(ie,i+6) = d24(i+6)
         d(ie,i+12) = d24(i+12)
         d(ie,i+18) = d24(i+18)
  160 continue
c
      do 170 i=1,6
      do 170 j=1,6
         sinc(ib1+i,ib1+j) = sinc(ib1+i,ib1+j) + sk(i,j)
         sinc(ib1+i,ib2+j) = sinc(ib1+i,ib2+j) + sk(i,j+6)
         sinc(ib1+i,ib3+j) = sinc(ib1+i,ib3+j) + sk(i,j+12)
         sinc(ib1+i,ib4+j) = sinc(ib1+i,ib4+j) + sk(i,j+18)
         sinc(ib2+i,ib1+j) = sinc(ib2+i,ib1+j) + sk(i+6,j)
         sinc(ib2+i,ib2+j) = sinc(ib2+i,ib2+j) + sk(i+6,j+6)
         sinc(ib2+i,ib3+j) = sinc(ib2+i,ib3+j) + sk(i+6,j+12)
         sinc(ib2+i,ib4+j) = sinc(ib2+i,ib4+j) + sk(i+6,j+18)
         sinc(ib3+i,ib1+j) = sinc(ib3+i,ib1+j) + sk(i+12,j)
         sinc(ib3+i,ib2+j) = sinc(ib3+i,ib2+j) + sk(i+12,j+6)
         sinc(ib3+i,ib3+j) = sinc(ib3+i,ib3+j) + sk(i+12,j+12)
         sinc(ib3+i,ib4+j) = sinc(ib3+i,ib4+j) + sk(i+12,j+18)
         sinc(ib4+i,ib1+j) = sinc(ib4+i,ib1+j) + sk(i+18,j)
         sinc(ib4+i,ib2+j) = sinc(ib4+i,ib2+j) + sk(i+18,j+6)
         sinc(ib4+i,ib3+j) = sinc(ib4+i,ib3+j) + sk(i+18,j+12)
         sinc(ib4+i,ib4+j) = sinc(ib4+i,ib4+j) + sk(i+18,j+18)
  170 continue
 2000 continue
c
      do 180 i=1,ndof
         sinc(ndof+1,i) = ddl(i) / r0
c         ddl(i) = 0.d0
         sinc(i,ndof+1) = -forc(i)
  180 continue
      sinc(ndof+1,ndof+1) = dpl / r0
c
      call clrv (ndof+1,xb)
      xb(ndof+1) = r
      r0 = r
      dpl = 0.d0
c
      do 190 i=1,ndof
      do 190 j=1,ndof
         sinc(i,j) = bc(i) * sinc(i,j)
         sinc(j,i) = bc(i) * sinc(j,i)
  190 continue
      do 200 i=1,ndof
         if(sinc(i,i) .ne. 0.d0) go to 200
            sinc(i,i) = 1.d0
  200 continue
c
      call clrv (ndof+1,dd)
      call solve (sinc,dd,xb,dp,det,idt,ndof,ndof+1,jdt)
c
      dp = dd(ndof+1)
c
c********************************************************************
c     Newton-Raphson収斂過程
c********************************************************************
 1100 nofcs = nofcs + 1
c
      do 210 i=1,ndof
      ddl(i) = ddl(i) + dd(i)
      ddi2 = ddi2 + dd(i)*dd(i)
      ddl2 = ddl2 + ddl(i)*ddl(i)
  210 continue
c
      dpl = dpl + dp
      p = p + dp
c
      if(ddl2.eq.0) go to 1200
      arc = ddi2 / ddl2
      ddi2 = 0.d0
      ddl2 = 0.d0
      if ( nofcs.gt.6 ) then
c         write(*,*)'6ステップで収斂しないので、弧長を半分にする。'
          nlet = nlet + 1
          konv = konv + 1
          if( nlet.gt.10 ) then
             p1 = p
             go to 4020
           end if
          go to 4000
      end if
      write(*,3020) arc
 3020 format(5x,'arc = ',1pe15.7)
c
      if(arc.lt.1.d-10) go to 1300
c
 1200 call clrmx (ndof+1,sinc)
      call clrv (ndof+1,ba)
c
      do 220 i=1,ndof
         ba(i) = p * forc(i)
  220 continue
c
      do 2100 ie=1,nele
      do 240 i = 1,4
         jcon(i) = nodn(ie, i)
  240 continue
         ib1 = (jcon(1) - 1) * 6
         ib2 = (jcon(2) - 1) * 6
         ib3 = (jcon(3) - 1) * 6
         ib4 = (jcon(4) - 1) * 6
cc
      a = cor( jcon(2), 1 ) - cor( jcon(1), 1 )
      b = cor( jcon(4), 2 ) - cor( jcon(1), 2 )
      call smx(s,th,e,um,a,b)
c
      do 810 i=1,3
         alp3(i) = alp((ie-1)*3+i)
  810 continue
c
      do 250 i=1,6
         dd24(i) = dd(ib1+i)
         dd24(i+6) = dd(ib2+i)
         dd24(i+12) = dd(ib3+i)
         dd24(i+18) = dd(ib4+i)
         d24(i) = d(ie,i)
         d24(i+6) = d(ie,i+6)
         d24(i+12) = d(ie,i+12)
         d24(i+18) = d(ie,i+18)
         ba24(i) = ba(ib1+i)
         ba24(i+6) = ba(ib2+i)
         ba24(i+12) = ba(ib3+i)
         ba24(i+18) = ba(ib4+i)
  250 continue
c
      call stmx (s,sk,a,b,d24,dd24,ba24,alp3)
c
      do 260 i=1,6
         d(ie,i) = d24(i)
         d(ie,i+6) = d24(i+6)
         d(ie,i+12) = d24(i+12)
         d(ie,i+18) = d24(i+18)
         ba(ib1+i) = ba24(i)
         ba(ib2+i) = ba24(i+6)
         ba(ib3+i) = ba24(i+12)
         ba(ib4+i) = ba24(i+18)
  260 continue
c
      do 270 i=1,6
      do 270 j=1,6
         sinc(ib1+i,ib1+j) = sinc(ib1+i,ib1+j) + sk(i,j)
         sinc(ib1+i,ib2+j) = sinc(ib1+i,ib2+j) + sk(i,j+6)
         sinc(ib1+i,ib3+j) = sinc(ib1+i,ib3+j) + sk(i,j+12)
         sinc(ib1+i,ib4+j) = sinc(ib1+i,ib4+j) + sk(i,j+18)
         sinc(ib2+i,ib1+j) = sinc(ib2+i,ib1+j) + sk(i+6,j)
         sinc(ib2+i,ib2+j) = sinc(ib2+i,ib2+j) + sk(i+6,j+6)
         sinc(ib2+i,ib3+j) = sinc(ib2+i,ib3+j) + sk(i+6,j+12)
         sinc(ib2+i,ib4+j) = sinc(ib2+i,ib4+j) + sk(i+6,j+18)
         sinc(ib3+i,ib1+j) = sinc(ib3+i,ib1+j) + sk(i+12,j)
         sinc(ib3+i,ib2+j) = sinc(ib3+i,ib2+j) + sk(i+12,j+6)
         sinc(ib3+i,ib3+j) = sinc(ib3+i,ib3+j) + sk(i+12,j+12)
         sinc(ib3+i,ib4+j) = sinc(ib3+i,ib4+j) + sk(i+12,j+18)
         sinc(ib4+i,ib1+j) = sinc(ib4+i,ib1+j) + sk(i+18,j)
         sinc(ib4+i,ib2+j) = sinc(ib4+i,ib2+j) + sk(i+18,j+6)
         sinc(ib4+i,ib3+j) = sinc(ib4+i,ib3+j) + sk(i+18,j+12)
         sinc(ib4+i,ib4+j) = sinc(ib4+i,ib4+j) + sk(i+18,j+18)
  270 continue
 2100 continue
c
      do 280 i=1,ndof
         sinc(ndof+1,i) = ddl(i)
         sinc(i,ndof+1) = -forc(i)
  280 continue
      sinc(ndof+1,ndof+1) = p
      ba(ndof+1) = 0.d0
      do 290 i=1,ndof
      do 290 j=1,ndof
         sinc(i,j) = bc(i) * sinc(i,j)
         sinc(j,i) = bc(i) * sinc(j,i)
         ba(i) = bc(i) * ba(i)
  290 continue
      do 300 i=1,ndof
         if(sinc(i,i) .ne. 0.d0) go to 300
            sinc(i,i) = 1.d0
  300 continue
      do 700 i=1,ndof
      do 700 j=1,ndof
         snr2(i,j) = sinc(i,j)
  700 continue
c
      call solve (sinc,dd,ba,dp,det,idt,ndof,ndof+1,jdt)
      dp = dd(ndof+1)
      go to 1100
c
c********************************************************************
c     計算結果表示
c********************************************************************
 1300 continue
      if( abs((p1-p)/p).lt.1.d-7 ) then
         write(8,*)
         write(8,5100)
         write(8,5111) p
 5100 format('---------------------- bucled ------------------------')
 5111 format('converged p = ',1pe15.7)
         stop
      end if
c
      arc = 0.d0
      nlet = 0
      if( konv.gt.1 ) rn = rn * 1.5d0
      if( jdt.eq.0 ) nccl = 1
      write(*,*)'----------------------------------------------------'
      write(*,3010) jdt, p
 3010 format(3x,'jdt = ',i2,5x,'p = ',1pe15.7)
c
c
      do 950 i=1,np
         write(20+i,5001) d(mds(i), mmds(i)), p
  950 continue
c
 5000 format(i3,3x,1pd15.7,3x,1pd15.7,3x,1pd15.7)
 5001 format(1pe15.7,5x,1pe15.7)
c
      if( jdt.ge.1 ) then
         write(*,*)'buckled.....'
         nbcl = 1
         nacl = 0
      end if
c
 4020 continue
      p1 = p
      if( nacl.eq.1 ) rn = rn / 2.d0
      if( nbcl.eq.1 ) then
         nbcl = 0
         nacl = 1
         go to 4000
      end if
      go to 4010
c
 4000 continue
      read(3) dd,d,ddl,dp,dpl,r0,p
      rewind 3
      rn = rn / 2.d0
c
      if( nccl.eq.0 ) then
         r = r * 0.9d0
         r0 = r
         dpl = r0
      end if
c
 4010 continue
      if( nccl.eq.1 ) then
      write(3) dd,d,ddl,dp,dpl,r0,p
      rewind 3
         r = rn
      end if
c
c      if ( p.gt.8000) then
c         stop
c      end if
c
      goto 1000
      end
c
c---------------------- main program end ----------------------------
c
c====================================================================
      subroutine stmx (s,sk,a,b,d,dd,ba,alp)
      implicit real*8(a-h,o-z)
c
      dimension d(24), dl(24), ta(24,24), tb(24,24), tg(24,24),
     &          dt0(24,24,24), sk(24,24), cf1(6,6),te(24,24),
     &          cf2(6,6), cf3(6,6), cf4(6,6), cf5(6,6), cf6(6,6),
     &          tngnt1(24,24), tngnt2(24,24), tngnt3(24,24), 
     &          tl(24,24), s(24,24), alp(3),
     &          dmy1(24), dmy2(24), dmy3(24,24), dmy4(24,24),
     &          dmy5(24,24), dmy6(24,24), cf(24,24), sk1(24,24),
     &          ba(24), dd(24), sk2(24,24), ti(24,24), sk3(24,24),
     &          dmy7(24), dmy8(24)
c
c********************************************************************
c
      call clrv (24,dmy1)
      call clrmx (24,tl)
      call clrmx (24,te)
      call clrmx (24,ti)
      call trns (ti,alp(1),alp(2),alp(3))
      call euler (te,d(4),d(5),d(10),d(11),d(16),d(17),d(22),d(23))
c
      call mxtv (24,ti,dd,dmy1)
      call mxv (24,te,dmy1,dd)
c
      do 130 i=1,24
         d(i) = d(i) + dd(i)
 130  continue
c
c********************************************************************
c
      u1 = d(1)
      v1 = d(2)
      w1 = d(3)
      a1 = d(4)
      b1 = d(5)
      g1 = d(6)
      u2 = d(7)
      v2 = d(8)
      w2 = d(9)
      a2 = d(10)
      b2 = d(11)
      g2 = d(12)
      u3 = d(13)
      v3 = d(14)
      w3 = d(15)
      a3 = d(16)
      b3 = d(17)
      g3 = d(18)
      u4 = d(19)
      v4 = d(20)
      w4 = d(21)
      a4 = d(22)
      b4 = d(23)
      g4 = d(24)
c
      ca1 = cos(a1)
      cb1 = cos(b1)
      cg1 = cos(g1)
      ca2 = cos(a2)
      cb2 = cos(b2)
      cg2 = cos(g2)
      ca3 = cos(a3)
      cb3 = cos(b3)
      cg3 = cos(g3)
      ca4 = cos(a4)
      cb4 = cos(b4)
      cg4 = cos(g4)
c
      sa1 = sin(a1)
      sb1 = sin(b1)
      sg1 = sin(g1)
      sa2 = sin(a2)
      sb2 = sin(b2)
      sg2 = sin(g2)
      sa3 = sin(a3)
      sb3 = sin(b3)
      sg3 = sin(g3)
      sa4 = sin(a4)
      sb4 = sin(b4)
      sg4 = sin(g4)
c
      a12 = (a1+a2) / 2.d0
      a13 = (a1+a3) / 2.d0
      a14 = (a1+a4) / 2.d0
      b12 = (b1+b2) / 2.d0
      b13 = (b1+b3) / 2.d0
      b14 = (b1+b4) / 2.d0
      g12 = (g1+g2) / 2.d0
      g13 = (g1+g3) / 2.d0
      g14 = (g1+g4) / 2.d0
c
      caa2 = cos(a12)
      caa3 = cos(a13)
      caa4 = cos(a14)
      cbb2 = cos(b12)
      cbb3 = cos(b13)
      cbb4 = cos(b14)
      cgg2 = cos(g12)
      cgg3 = cos(g13)
      cgg4 = cos(g14)
      saa2 = sin(a12)
      saa3 = sin(a13)
      saa4 = sin(a14)
      sbb2 = sin(b12)
      sbb3 = sin(b13)
      sbb4 = sin(b14)
      sgg2 = sin(g12)
      sgg3 = sin(g13)
      sgg4 = sin(g14)
c
c********************************************************************
c
      call clrv (24,dl)
c
      dl(7)  = u2-u1+a-a*cg1*cb1
      dl(8)  = v2-v1-a*sg1*ca1-a*cg1*sb1*sa1
      dl(9)  = w2-w1-a*sg1*sa1+a*cg1*sb1*ca1
      dl(10) = (a2-a1)+(g2-g1)*sbb2
      dl(11) = -(g2-g1)*cbb2*saa2+(b2-b1)*caa2
      dl(12) = (g2-g1)*cbb2*caa2+(b2-b1)*saa2
      dl(13) = u3-u1+a-a*cg1*cb1+b*sg1*cb1
      dl(14) = v3-v1-a*sg1*ca1-a*cg1*sb1*sa1+b-b*cg1*ca1+b*sg1*sb1*sa1
      dl(15) = w3-w1-a*sg1*sa1+a*cg1*sb1*ca1-b*cg1*sa1-b*sg1*sb1*ca1
      dl(16) = (a3-a1)+(g3-g1)*sbb3
      dl(17) = -(g3-g1)*cbb3*saa3+(b3-b1)*caa3
      dl(18) = (g3-g1)*cbb3*caa3+(b3-b1)*saa3
      dl(19) = u4-u1+b*sg1*cb1
      dl(20) = v4-v1+b-b*cg1*ca1+b*sg1*sb1*sa1
      dl(21) = w4-w1-b*cg1*sa1-b*sg1*sb1*ca1
      dl(22) = (a4-a1)+(g4-g1)*sbb4
      dl(23) = -(g4-g1)*cbb4*saa4+(b4-b1)*caa4
      dl(24) = (g4-g1)*cbb4*caa4+(b4-b1)*saa4
c
c********************************************************************
c
      call clrmx (24,ta)
      call clrmx (24,tb)
      call clrmx (24,tg)
c
      ta(1,1) = 0.d0
      ta(1,2) = 0.d0
      ta(1,3) = 0.d0
      ta(2,1) = -sg1*sa1+cg1*sb1*ca1
      ta(2,2) = -cg1*sa1-ca1*sg1*sb1
      ta(2,3) = -cb1*ca1
      ta(3,1) = sg1*ca1+cg1*sb1*sa1
      ta(3,2) = ca1*cg1-sg1*sb1*sa1
      ta(3,3) = -cb1*sa1
      tb(1,1) = -cg1*sb1
      tb(1,2) = sb1*sg1
      tb(1,3) = cb1
      tb(2,1) = cg1*cb1*sa1
      tb(2,2) = -sa1*sg1*cb1
      tb(2,3) = sb1*sa1
      tb(3,1) = -cg1*cb1*ca1
      tb(3,2) = sg1*cb1*ca1
      tb(3,3) = -sb1*ca1
      tg(1,1) = -sg1*cb1
      tg(1,2) = -cg1*cb1
      tg(1,3) = 0.d0
      tg(2,1) = cg1*ca1-sg1*sb1*sa1
      tg(2,2) = -sg1*ca1-cg1*sb1*sa1
      tg(2,3) = 0.d0
      tg(3,1) = cg1*sa1+sg1*sb1*ca1
      tg(3,2) = -sg1*sa1+cg1*sb1*ca1
      tg(3,3) = 0.d0
c
      do 10 i=1,3
      do 10 j=1,3
         ta(i+3,j+3) = ta(i,j)
         ta(i+6,j+6) = ta(i,j)
         ta(i+9,j+9) = ta(i,j)
         ta(i+12,j+12) = ta(i,j)
         ta(i+15,j+15) = ta(i,j)
         ta(i+18,j+18) = ta(i,j)
         ta(i+21,j+21) = ta(i,j)
         tb(i+3,j+3) = tb(i,j)
         tb(i+6,j+6) = tb(i,j)
         tb(i+9,j+9) = tb(i,j)
         tb(i+12,j+12) = tb(i,j)
         tb(i+15,j+15) = tb(i,j)
         tb(i+18,j+18) = tb(i,j)
         tb(i+21,j+21) = tb(i,j)
         tg(i+3,j+3) = tg(i,j)
         tg(i+6,j+6) = tg(i,j)
         tg(i+9,j+9) = tg(i,j)
         tg(i+12,j+12) = tg(i,j)
         tg(i+15,j+15) = tg(i,j)
         tg(i+18,j+18) = tg(i,j)
         tg(i+21,j+21) = tg(i,j)
   10 continue
c
      call clmx3 (24,dt0)
c
      do 20 i=1,24
      do 20 j=1,24
         dt0(4,i,j) = ta(i,j)
         dt0(5,i,j) = tb(i,j)
         dt0(6,i,j) = tg(i,j)
   20 continue
c
c********************************************************************
c     coefficient of relative displacement vector
c
      call clrmx (6,cf1)
      call clrmx (6,cf2)
      call clrmx (6,cf3)
      call clrmx (6,cf4)
      call clrmx (6,cf5)
      call clrmx (6,cf6)
c
      cf1(1,1) = -1.d0
      cf1(1,5) = a*cg1*sb1
      cf1(1,6) = a*sg1*cb1
      cf1(2,2) = -1.d0
      cf1(2,4) = a*sg1*sa1-a*cg1*sb1*ca1
      cf1(2,5) = -a*cg1*cb1*sa1
      cf1(2,6) = -a*cg1*ca1+a*sg1*sb1*sa1
      cf1(3,3) = -1.d0
      cf1(3,4) = -a*sg1*ca1-a*cg1*sb1*sa1
      cf1(3,5) = a*cg1*cb1*ca1
      cf1(3,6) = -a*cg1*sa1-a*sg1*sb1*ca1
      cf1(4,4) = -1.d0
      cf1(4,5) = (g2-g1)/2.d0*cbb2
      cf1(4,6) = -sbb2
      cf1(5,4) = -(g2-g1)/2.d0*cbb2*caa2-(b2-b1)/2.d0*saa2
      cf1(5,5) = (g2-g1)/2.d0*sbb2*saa2-caa2
      cf1(5,6) = cbb2*saa2
      cf1(6,4) = -(g2-g1)/2.d0*cbb2*saa2+(b2-b1)/2.d0*caa2
      cf1(6,5) = -saa2-(g2-g1)/2.d0*sbb2*caa2
      cf1(6,6) = -cbb2*caa2
c
      cf2(1,1) = 1.d0
      cf2(2,2) = 1.d0
      cf2(3,3) = 1.d0
      cf2(4,4) = 1.d0
      cf2(4,5) = cf1(4,5)
      cf2(4,6) = -cf1(4,6)
      cf2(5,4) = cf1(5,4)
      cf2(5,5) = (g2-g1)/2.d0*sbb2*saa2+caa2
      cf2(5,6) = -cf1(5,6)
      cf2(6,4) = cf1(6,4)
      cf2(6,5) = -(g2-g1)/2.d0*sbb2*caa2+saa2
      cf2(6,6) = -cf1(6,6)
c
      cf3(1,1) = -1.d0
      cf3(1,5) = a*cg1*sb1-b*sg1*sb1
      cf3(1,6) = a*sg1*cb1+b*cg1*cb1
      cf3(2,2) = -1.d0
      cf3(2,4) = a*sg1*sa1-a*cg1*sb1*ca1+b*cg1*sa1+b*sg1*sb1*ca1
      cf3(2,5) = -a*cg1*cb1*sa1+b*sg1*cb1*sa1
      cf3(2,6) = -a*cg1*ca1+a*sg1*sb1*sa1+b*sg1*ca1+b*cg1*sb1*sa1
      cf3(3,3) = -1.d0
      cf3(3,4) = -a*sg1*ca1-a*cg1*sb1*sa1-b*cg1*ca1+b*sg1*sb1*sa1
      cf3(3,5) = a*cg1*cb1*ca1-b*sg1*cb1*ca1
      cf3(3,6) = -a*cg1*sa1-a*sg1*sb1*ca1+b*sg1*sa1-b*cg1*sb1*ca1
      cf3(4,4) = -1.d0
      cf3(4,5) = (g3-g1)/2.d0*cbb3
      cf3(4,6) = -sbb3
      cf3(5,4) = -(g3-g1)/2.d0*cbb3*caa3-(b3-b1)/2.d0*saa3
      cf3(5,5) = (g3-g1)/2.d0*sbb3*saa3-caa3
      cf3(5,6) = cbb3*saa3
      cf3(6,4) = -(g3-g1)/2.d0*cbb3*saa3+(b3-b1)/2.d0*caa3
      cf3(6,5) = -saa3-(g3-g1)/2.d0*sbb3*caa3
      cf3(6,6) = -cbb3*caa3
c
      cf4(1,1) = 1.d0
      cf4(2,2) = 1.d0
      cf4(3,3) = 1.d0
      cf4(4,4) = 1.d0
      cf4(4,5) = cf3(4,5)
      cf4(4,6) = -cf3(4,6)
      cf4(5,4) = cf3(5,4)
      cf4(5,5) = (g3-g1)/2.d0*sbb3*saa3+caa3
      cf4(5,6) = -cf3(5,6)
      cf4(6,4) = cf3(6,4)
      cf4(6,5) = -(g3-g1)/2.d0*sbb3*caa3+saa3
      cf4(6,6) = -cf3(6,6)
c
      cf5(1,1) = -1.d0
      cf5(1,5) = -b*sg1*sb1
      cf5(1,6) = b*cg1*cb1
      cf5(2,2) = -1.d0
      cf5(2,4) = b*cg1*sa1+b*sg1*sb1*ca1
      cf5(2,5) = b*sg1*cb1*sa1
      cf5(2,6) = b*sg1*ca1+b*cg1*sb1*sa1
      cf5(3,3) = -1.d0
      cf5(3,4) = -b*cg1*ca1+b*sg1*sb1*sa1
      cf5(3,5) = -b*sg1*cb1*ca1
      cf5(3,6) = b*sg1*sa1-b*cg1*sb1*ca1
      cf5(4,4) = -1.d0
      cf5(4,5) = (g4-g1)/2.d0*cbb4
      cf5(4,6) = -sbb4
      cf5(5,4) = -(g4-g1)/2.d0*cbb4*caa4-(b4-b1)/2.d0*saa4
      cf5(5,5) = (g4-g1)/2.d0*sbb4*saa4-caa4
      cf5(5,6) = cbb4*saa4
      cf5(6,4) = -(g4-g1)/2.d0*cbb4*saa4+(b4-b1)/2.d0*caa4
      cf5(6,5) = -saa4-(g4-g1)/2.d0*sbb4*caa4
      cf5(6,6) = -cbb4*caa4
c
      cf6(1,1) = 1.d0
      cf6(2,2) = 1.d0
      cf6(3,3) = 1.d0
      cf6(4,4) = 1.d0
      cf6(4,5) = cf5(4,5)
      cf6(4,6) = -cf5(4,6)
      cf6(5,4) = cf5(5,4)
      cf6(5,5) = (g4-g1)/2.d0*sbb4*saa4+caa4
      cf6(5,6) = -cf5(5,6)
      cf6(6,4) = cf5(6,4)
      cf6(6,5) = -(g4-g1)/2.d0*sbb4*caa4+saa4
      cf6(6,6) = -cf5(6,6)
c
      call clrmx (24,cf)
      do 1000 i=1,6
      do 1000 j=1,6
         cf(i+6,j) = cf1(i,j)
         cf(i+6,j+6) = cf2(i,j)
         cf(i+12,j) = cf3(i,j)
         cf(i+12,j+12) = cf4(i,j)
         cf(i+18,j) = cf5(i,j)
         cf(i+18,j+18) = cf6(i,j)
 1000 continue
c
      call trns(tl,d(4),d(5),d(6))
c
      call clrmx (24,tngnt1)
      call clrmx (24,tngnt2)
      call clrmx (24,tngnt3)
      call clrv (24,dmy1)
      call clrv (24,dmy2)
      call clrmx (24,dmy3)
      call clrmx (24,dmy4)
      call clrmx (24,dmy5)
      call clrmx (24,dmy6)
c
c********************************************************************
c
      call mxtv (24,tl,dl,dmy1)
      call mxv (24,s,dmy1,dmy2)
      do 30 k=1,24
      do 30 i=1,24
      do 30 j=1,24
         tngnt1(i,k) = tngnt1(i,k) + dt0(k,i,j) * dmy2(j)
   30 continue
c
c********************************************************************
c
      do 40 k=1,24
      do 40 i=1,24
      do 40 j=1,24
         dmy3(i,k) = dmy3(i,k) + dt0(k,j,i) * dl(j)
   40 continue
c
      call mxmx (24,s,dmy3,dmy4)
      call mxmx (24,tl,dmy4,tngnt2)
c
c********************************************************************
c
      call mxtmx (24,tl,cf,dmy5)
      call mxmx (24,s,dmy5,dmy6)
      call mxmx(24,tl,dmy6,tngnt3)
c
c********************************************************************
c
      call clrmx (24,sk1)
      call clrmx (24,sk2)
      call clrmx (24,sk3)
      call clrmx (24,sk)
c
      do 80 i=1,24
      do 80 j=1,24
         sk1(i,j) = tngnt1(i,j) + tngnt2(i,j) + tngnt3(i,j)
   80 continue
c
      call euler (te,d(4),d(5),d(10),d(11),d(16),d(17),d(22),d(23))
      call mxmx (24,ti,sk1,sk2)
      call mxmx (24,sk2,te,sk3)
      call mxmxt (24,sk3,ti,sk)
c
c********************************************************************
c
      call clrv (24,dmy1)
      call clrv (24,dmy2)
      call clrv (24,dmy7)
      call clrv (24,dmy8)
      call mxtv (24,tl,dl,dmy1)
      call mxv (24,s,dmy1,dmy2)
      call mxv (24,tl,dmy2,dmy7)
      call mxv (24,ti,dmy7,dmy8)
      do 120 i=1,24
         ba(i) = ba(i) - dmy8(i)
  120 continue
c
      return
      end
c
c====================================================================
       subroutine mxtmx(n,a,b,c)
       implicit real*8(a-h,o-z)
       dimension a(n,n),b(n,n),c(n,n)
       do 10 i=1,n
       do 10 k=1,n
       c(i,k)=0.d0
       do 10 j=1,n
       c(i,k)=c(i,k)+a(j,i)*b(j,k)
10     continue
       return
       end
c
c====================================================================
       subroutine mxmx(n,a,b,c)
       implicit real*8(a-h,o-z)
       dimension a(n,n),b(n,n),c(n,n)
       do 10 i=1,n
       do 10 k=1,n
       c(i,k)=0.d0
       do 10 j=1,n
       c(i,k)=c(i,k)+a(i,j)*b(j,k)
10     continue
       return
       end
c
c====================================================================
       subroutine mxtv(n,a,b,c)
       implicit real*8(a-h,o-z)
       dimension a(n,n),b(n),c(n)
       do 10 i=1,n
       c(i)=0.d0
       do 10 j=1,n
       c(i)=c(i)+a(j,i)*b(j)
10     continue
       return
       end
c
c====================================================================
       subroutine mxv(n,a,b,c)
       implicit real*8(a-h,o-z)
       dimension a(n,n),b(n),c(n)
       do 10 i=1,n
       c(i)=0.d0
       do 10 j=1,n
       c(i)=c(i)+a(i,j)*b(j)
10     continue
       return
       end
c
c====================================================================
       subroutine mxmxt(n,a,b,c)
       implicit real*8(a-h,o-z)
       dimension a(n,n),b(n,n),c(n,n)
       do 10 i=1,n
       do 10 k=1,n
       c(i,k)=0.d0
       do 10 j=1,n
       c(i,k)=c(i,k)+a(i,j)*b(k,j)
10     continue
       return
       end
c
c====================================================================
      subroutine smx(s,th,e,um,a,b)
      implicit real*8(a-h,o-z)
      dimension s(24,24)
c
      call clrmx (24,s)
c
      p1 = ( b / ( 3.d0 * a ) + ( 1.d0 - um ) * a / ( 6.d0 * b ) )
     &     * e * th / ( 1.d0 - um * um )
      p2 = ( um / 4.d0 + ( 1 - um ) / 8.d0 )
     &     * e * th / ( 1.d0 - um * um )
      p3 = ( a / ( 3.d0 * b ) + ( 1.d0 - um ) * b / ( 6.d0 * a ) )
     &     * e * th / ( 1.d0 - um * um )
      p4 = ( - b / ( 3.d0 * a ) + ( 1.d0 - um ) * a / ( 12.d0 * b ) )
     &     * e * th / ( 1 - um * um )
      p5 = ( um / 4.d0 - ( 1.d0 - um ) / 8.d0 )
     &     * e * th / ( 1.d0 - um * um )
      p6 = ( a / ( 6.d0 * b ) - ( 1.d0 - um ) * b / ( 6.d0 
     &     * a ) ) * e * th / ( 1.d0 - um * um )
      p7 = ( b / ( 6.d0 * a ) - ( 1.d0 - um ) * a / ( 6.d0 * b ) )
     &     * e * th / ( 1.d0 - um * um )
      p8 = ( - a / ( 3.d0 * b ) + ( 1.d0 - um ) * b / ( 12.d0 * a ) )
     &     * e * th / ( 1.d0 - um * um )
c
      s(1,1) = p1
      s(1,7) = p4
      s(1,13) = -p1 / 2.d0
      s(1,19) = p7
      s(1,2) = p2
      s(1,8) = p5
      s(1,14) = -p2
      s(1,20) = -p5
      s(7,1) = p4
      s(7,7) = p1
      s(7,13) = p7
      s(7,19) = -p1 / 2.d0
      s(7,2) = -p5
      s(7,8) = -p2
      s(7,14) = p5
      s(7,20) = p2
      s(13,1) = -p1 / 2.d0
      s(13,7) = p7
      s(13,13) = p1
      s(13,19) = p4
      s(13,2) = -p2
      s(13,8) = -p5
      s(13,14) = p2
      s(13,20) = p5
      s(19,1) = p7
      s(19,7) = -p1 / 2.d0
      s(19,13) = p4
      s(19,19) = p1
      s(19,2) = p5
      s(19,8) = p2
      s(19,14) = -p5
      s(19,20) = -p2
      s(2,1) = p2
      s(2,7) = -p5
      s(2,13) = -p2
      s(2,19) = p5
      s(2,2) = p3
      s(2,8) = p6
      s(2,14) = -p3 / 2.d0
      s(2,20) = p8
      s(8,1) = p5
      s(8,7) = -p2
      s(8,13) = -p5
      s(8,19) = p2
      s(8,2) = p6
      s(8,8) = p3
      s(8,14) = p8
      s(8,20) = -p3 / 2.d0
      s(14,1) = -p2
      s(14,7) = p5
      s(14,13) = p2
      s(14,19) = -p5
      s(14,2) = -p3 / 2.d0
      s(14,8) = p8
      s(14,14) = p3
      s(14,20) = p6
      s(20,1) = -p5
      s(20,7) = p2
      s(20,13) = p5
      s(20,19) = -p2
      s(20,2) = p8
      s(20,8) = -p3 / 2.d0
      s(20,14) = p6
      s(20,20) = p3
c
      p15 = e * th * th * th / ( 360.d0 * ( 1.d0 - um ** 2 ) * a * b )
      s(3,3) = 120.d0 * b * b / a / a + 120.d0 * a * a / b / b
     &       + 84.d0 - 24.d0 * um
      s(9,3) = - 120.d0 * b * b / a / a + 60.d0 * a * a / b / b
     &       - 84.d0 + 24.d0 * um
      s(15,3) = - 60.d0 * b * b / a / a - 60.d0 * a * a / b / b
     &        + 84.d0 - 24.d0 * um
      s(21,3) = 60.d0 * b * b / a / a - 120.d0 * a * a / b / b
     &        - 84.d0 + 24.d0 * um
      s(3,3) = s(3,3) * p15
      s(9,3) = s(9,3) * p15
      s(15,3) = s(15,3) * p15
      s(21,3) = s(21,3) * p15
      s(9,9) = s(3,3)
      s(15,15) = s(3,3)
      s(21,21) = s(3,3)
      s(3,9) = s(9,3)
      s(21,15) = s(9,3)
      s(15,21) = s(9,3)
      s(21,9) = s(15,3)
      s(3,15) = s(15,3)
      s(9,21) = s(15,3)
      s(15,9) = s(21,3)
      s(9,15) = s(21,3)
      s(3,21) = s(21,3)
c
      s(4,3) = 60.d0 * a * a / b + 6.d0 * b + 24.d0 * b * um
      s(10,3) = 30.d0 * a * a / b - 6.d0 * b - 24.d0 * b * um
      s(16,3) = 30.d0 * a * a / b - 6.d0 * b + 6.d0 * b * um
      s(22,3) = 60.d0 * a * a / b + 6.d0 * b - 6.d0 * b * um
      s(4,3) = s(4,3) * p15
      s(10,3) = s(10,3) * p15
      s(16,3) = s(16,3) * p15
      s(22,3) = s(22,3) * p15
      s(10,9) = s(4,3)
      s(16,15) = - s(4,3)
      s(22,21) = - s(4,3)
      s(4,9) = s(10,3)
      s(22,15) = - s(10,3)
      s(16,21) = - s(10,3)
      s(22,9) = s(16,3)
      s(4,15) = - s(16,3)
      s(10,21) = - s(16,3)
      s(16,9) = s(22,3)
      s(10,15) = - s(22,3)
      s(4,21) = - s(22,3)
      s(3,4) = s(4,3)
      s(3,10) = s(10,3)
      s(3,16) = s(16,3)
      s(3,22) = s(22,3)
      s(9,10) = s(10,9)
      s(15,16) = s(16,15)
      s(21,22) = s(22,21)
      s(9,4) = s(4,9)
      s(15,22) = s(22,15)
      s(21,16) = s(16,21)
      s(9,22) = s(22,9)
      s(15,4) = s(4,15)
      s(21,10) = s(10,21)
      s(9,16) = s(16,9)
      s(15,10) = s(10,15)
      s(21,4) = s(4,21)
c
      s(5,3) = - 60.d0 * b * b / a - 6.d0 * a - 24.d0 * um * a
      s(11,3) = - 60.d0 * b * b / a - 6.d0 * a + 6.d0 * um * a
      s(17,3) = - 30.d0 * b * b / a + 6.d0 * a - 6.d0 * um * a
      s(23,3) = - 30.d0 * b * b / a + 6.d0 * a + 24.d0 * um * a
      s(5,3) = s(5,3) * p15
      s(11,3) = s(11,3) * p15
      s(17,3) = s(17,3) * p15
      s(23,3) = s(23,3) * p15
      s(11,9) = - s(5,3)
      s(17,15) = - s(5,3)
      s(23,21) = s(5,3)
      s(5,9) = - s(11,3)
      s(23,15) = - s(11,3)
      s(17,21) = s(11,3)
      s(23,9) = - s(17,3)
      s(5,15) = - s(17,3)
      s(11,21) = s(17,3)
      s(17,9) = - s(23,3)
      s(11,15) = - s(23,3)
      s(5,21) = s(23,3)
      s(3,5) = s(5,3)
      s(3,11) = s(11,3)
      s(3,17) = s(17,3)
      s(3,23) = s(23,3)
      s(9,11) = s(11,9)
      s(15,17) = s(17,15)
      s(21,23) = s(23,21)
      s(9,5) = s(5,9)
      s(15,23) = s(23,15)
      s(21,17) = s(17,21)
      s(9,23) = s(23,9)
      s(15,5) = s(5,15)
      s(21,11) = s(11,21)
      s(9,17) = s(17,9)
      s(15,11) = s(11,15)
      s(21,5) = s(5,21)
c
      s(4,4) = 8.d0 * b * b + 40.d0 * a * a - 8.d0 * b * b * um
      s(10,4) = - 8.d0 * b * b + 20.d0 * a * a + 8.d0 * b * b * um
      s(16,4) = 2.d0 * b * b + 10.d0 * a * a - 2.d0 * b * b * um
      s(22,4) = - 2.d0 * b * b + 20.d0 * a * a + 2.d0 * b * b * um
      s(4,4) = s(4,4) * p15
      s(10,4) = s(10,4) * p15
      s(16,4) = s(16,4) * p15
      s(22,4) = s(22,4) * p15
      s(10,10) = s(4,4)
      s(16,16) = s(4,4)
      s(22,22) = s(4,4)
      s(4,10) = s(10,4)
      s(22,16) = s(10,4)
      s(16,22) = s(10,4)
      s(22,10) = s(16,4)
      s(4,16) = s(16,4)
      s(10,22) = s(16,4)
      s(16,10) = s(22,4)
      s(10,16) = s(22,4)
      s(4,22) = s(22,4)
c
      s(5,4) = - 30.d0 * um * a * b
      s(11,4) = 0.d0
      s(17,4) = 0.d0
      s(23,4) = 0.d0
      s(5,4) = s(5,4) * p15
      s(11,4) = s(11,4) * p15
      s(17,4) = s(17,4) * p15
      s(23,4) = s(23,4) * p15
      s(11,10) = - s(5,4)
      s(17,16) = s(5,4)
      s(23,22) = - s(5,4)
      s(5,10) = - s(11,4)
      s(23,16) = s(11,4)
      s(17,22) = - s(11,4)
      s(23,10) = - s(17,4)
      s(5,16) = s(17,4)
      s(11,22) = - s(17,4)
      s(17,10) = - s(23,4)
      s(11,16) = s(23,4)
      s(5,22) = - s(23,4)
      s(4,5) = s(5,4)
      s(4,11) = s(11,4)
      s(4,17) = s(17,4)
      s(4,23) = s(23,4)
      s(10,11) = s(11,10)
      s(16,17) = s(17,16)
      s(22,23) = s(23,22)
      s(10,5) = s(5,10)
      s(16,23) = s(23,16)
      s(22,17) = s(17,22)
      s(10,23) = s(23,10)
      s(16,5) = s(5,16)
      s(22,11) = s(11,22)
      s(10,17) = s(17,10)
      s(16,11) = s(11,16)
      s(22,5) = s(5,22)
c
      s(5,5) = 40.d0 * b * b + 8.d0 * a * a - 8.d0 * a * a * um
      s(11,5) = 20.d0 * b * b - 2.d0 * a * a + 2.d0 * a * a * um
      s(17,5) = 10.d0 * b * b + 2.d0 * a * a - 2.d0 * a * a * um
      s(23,5) = 20.d0 * b * b - 8.d0 * a * a + 8.d0 * a * a * um
      s(5,5) = s(5,5) * p15
      s(11,5) = s(11,5) * p15
      s(17,5) = s(17,5) * p15
      s(23,5) = s(23,5) * p15
      s(11,11) = s(5,5)
      s(17,17) = s(5,5)
      s(23,23) = s(5,5)
      s(5,11) = s(11,5)
      s(23,17) = s(11,5)
      s(17,23) = s(11,5)
      s(23,11) = s(17,5)
      s(5,17) = s(17,5)
      s(11,23) = s(17,5)
      s(17,11) = s(23,5)
      s(11,17) = s(17,11)
      s(5,23) = s(17,11)
c
      p9 = e * th * a * b * 0.03d0
         s(6, 6) = p9
         s(6, 12) = -p9 / 3.d0
         s(6, 18) = -p9 / 3.d0
         s(6, 24) = -p9 / 3.d0
         s(12, 6) = -p9 / 3.d0
         s(12, 12) = p9
         s(12, 18) = -p9 / 3.d0
         s(12, 24) = -p9 / 3.d0
         s(18, 6) = -p9 / 3.d0
         s(18, 12) = -p9 / 3.d0
         s(18, 18) = p9
         s(18, 24) = -p9 / 3.d0
         s(24, 6) = -p9 / 3.d0
         s(24, 12) = -p9 / 3.d0
         s(24, 18) = -p9 / 3.d0
         s(24, 24) = p9
c
      return
      end
c====================================================================
      subroutine trns (tl,al,bt,gm)
      implicit real*8(a-h,o-z)
      dimension t(3,3), tl(24,24)
c
      do 110 i=1,24
      do 110 j=1,24
         tl(i,j) = 0.d0
  110 continue
c check!!
      t(1,1) = cos(gm)*cos(bt)
      t(1,2) = -sin(gm)*cos(bt)
      t(1,3) = sin(bt)
      t(2,1) = sin(gm)*cos(al)+cos(gm)*sin(bt)*sin(al)
      t(2,2) = cos(gm)*cos(al)-sin(gm)*sin(bt)*sin(al)
      t(2,3) = -cos(bt)*sin(al)
      t(3,1) = sin(gm)*sin(al)-cos(gm)*sin(bt)*cos(al)
      t(3,2) = cos(gm)*sin(al)+sin(gm)*sin(bt)*cos(al)
      t(3,3) = cos(bt)*cos(al)
c
      do 10 i=1,3
      do 10 j=1,3
         tl(i,j) = t(i,j)
         tl(i+3,j+3) = t(i,j)
         tl(i+6,j+6) = t(i,j)
         tl(i+9,j+9) = t(i,j)
         tl(i+12,j+12) = t(i,j)
         tl(i+15,j+15) = t(i,j)
         tl(i+18,j+18) = t(i,j)
         tl(i+21,j+21) = t(i,j)
   10 continue
c
      return
      end
c
c====================================================================
      subroutine euler(te,a1,b1,a2,b2,a3,b3,a4,b4)
      implicit real*8 (a-h,o-z)
      dimension te(24,24)
c
      do 10 i=1,24
      do 10 j=1,24
         te(i,j) = 0.d0
   10 continue
c
      do 20 i=1,24
         te(i,i)=1.d0
   20 continue
c
      te(4,5) = sin(a1)*tan(b1)
      te(4,6) = -cos(a1)*tan(b1)
      te(5,5) = cos(a1)
      te(5,6) = sin(a1)
      te(6,5) = -sin(a1)/cos(b1)
      te(6,6) = cos(a1)/cos(b1)
      te(10,11) = sin(a2)*tan(b2)
      te(10,12) = -cos(a2)*tan(b2)
      te(11,11) = cos(a2)
      te(11,12) = sin(a2)
      te(12,11) = -sin(a2)/cos(b2)
      te(12,12) = cos(a2)/cos(b2)
      te(16,17) = sin(a3)*tan(b3)
      te(16,18) = -cos(a3)*tan(b3)
      te(17,17) = cos(a3)
      te(17,18) = -sin(a3)
      te(18,17) = -sin(a3)/cos(b3)
      te(18,18) = cos(a3)/cos(b3)
      te(22,23) = sin(a4)*tan(b4)
      te(22,24) = -cos(a4)*tan(b4)
      te(23,23) = cos(a4)
      te(23,24) = sin(a4)
      te(24,23) = -sin(a4)/cos(b4)
      te(24,24) = cos(a4)/cos(b4)
      return
      end
c
c====================================================================
        subroutine clrv(n,a)
        implicit real*8(a-h,o-z)
         dimension a(n)
          do 10 i=1,n
             a(i)=0.d0
10        continue
         return
        end
c
c====================================================================
      subroutine clmx3 (n,a)
      implicit real*8(a-h,o-z)
      dimension a(n,n,n)
      do 10 k=1,n
      do 10 i=1,n
      do 10 j=1,n
         a(k,i,j) = 0.d0
   10 continue
      return
      end
c
c====================================================================
      subroutine clrmx (n,a)
      implicit real*8(a-h,o-z)
c
      dimension a(n,n)
      do 10 i=1,n
      do 10 j=1,n
         a(i,j)=0.d0
   10 continue
c
      return
      end
c====================================================================
      subroutine solve (sinc,dd,df,dp,det,idt,n0,n01,jdt)
c
c  Gaussian Elimination Method
c
      implicit real*8 (a-h,o-z)
      dimension sinc(n01,n01),dd(n01),df(n01)
c
      det=1.d0
      idt=0
      jdt=0
c                                 Forward Elimination
      do 100 i=1,n0
         inext=i+1
         if(sinc(i,i).le.0.d0) jdt=jdt+1
         det=det*sinc(i,i)
         j=idint(dlog10(dabs(det)))
         idt=idt+j
         det=det/10.d0**(j)
c
         sinc(i,i)=1.d0/sinc(i,i)
         df(i)=df(i)*sinc(i,i)
         do 110 j=inext,n01
            sinc(i,j)=sinc(i,j)*sinc(i,i)
            df(j)=df(j)-sinc(j,i)*df(i)
            do 110 k=inext,n01
               sinc(k,j)=sinc(k,j)-sinc(k,i)*sinc(i,j)
  110    continue
  100 continue
c                                 Backward Substitution
      dd(n01)=df(n01)/sinc(n01,n01)
      dp=dd(n01)
      do 120 i=1,n0
         k=n01-i
         dd(k)=df(k)
         knext=k+1
         do 120 j=knext,n01
            dd(k)=dd(k)-sinc(k,j)*dd(j)
  120 continue
c
      return
      end
