c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c+                                                                +
c+              4-node plate bending element program              +
c+                                           by T. Nishihara      +
c+                                                                +
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      program kasou
      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), alp(ndof)
      dimension  cor(ndof,2)
      dimension  nodn(ndof,4), mds(nout), mmds(nout), jcon(4)
      dimension  nbc(nbcc)
c
      character*20 file1
c
      open(7,file='input.dat')
      open(8,file='outkasou.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 = (nodn(ie, 1) - 1) * 6
         ib2 = (nodn(ie, 2) - 1) * 6
         ib3 = (nodn(ie, 3) - 1) * 6
         ib4 = (nodn(ie, 4) - 1) * 6
      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)
  150 continue
      do 151 i=1,24
         d24(i) = d(ie,i)
  151 continue
      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
      call stmx (s,sk,a,b,d24,dd24,ba24,alp3)
c
      do 160 i=1,24
         d(ie,i) = d24(i)
  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      write(*,*)'start Newton-Raphson'
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-6) 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 441 i = 1,4
         jcon(i) = nodn(ie, i)
  441 continue
         ib1 = (nodn(ie, 1) - 1) * 6
         ib2 = (nodn(ie, 2) - 1) * 6
         ib3 = (nodn(ie, 3) - 1) * 6
         ib4 = (nodn(ie, 4) - 1) * 6
      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)
         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
      do 251 i=1,24
         d24(i) = d(ie,i)
  251 continue
      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
      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
c
      call solve (sinc,dd,ba,dp,det,idt,ndof,ndof+1,jdt)
c
      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
      pi = 4.d0 * atan ( 1.d0 )
      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.d0) then
c         stop
c      end if
c
      go to 1000
      end
c====================================================================
      subroutine stmx (s,sk,a,b,d,dd,ba,alp)
      implicit real*8(a-h,o-z)
c
      dimension d(24), dl(24), dt0(24,24,24), sk(24,24), te(24,24),
     &          tl(24,24), s(24,24), alp(3), R(24,24),
     &          dmy1(24), dmy2(24), dmy3(24,24), dmy4(24,24),
     &          sk1(24,24), ba(24), dd(24), ti(24,24),
     &          dte(24,24,24), ste(3,3), dR(24,24,24),
     &          dR1(3,3), dR2(3,3)
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     変位ベクトル
      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
      caa2 = cos((a1+a2) / 2.d0)
      caa3 = cos((a1+a3) / 2.d0)
      caa4 = cos((a1+a4) / 2.d0)
      cbb2 = cos((b1+b2) / 2.d0)
      cbb3 = cos((b1+b3) / 2.d0)
      cbb4 = cos((b1+b4) / 2.d0)
      cgg2 = cos((g1+g2) / 2.d0)
      cgg3 = cos((g1+g3) / 2.d0)
      cgg4 = cos((g1+g4) / 2.d0)
      saa2 = sin((a1+a2) / 2.d0)
      saa3 = sin((a1+a3) / 2.d0)
      saa4 = sin((a1+a4) / 2.d0)
      sbb2 = sin((b1+b2) / 2.d0)
      sbb3 = sin((b1+b3) / 2.d0)
      sbb4 = sin((b1+b4) / 2.d0)
      sgg2 = sin((g1+g2) / 2.d0)
      sgg3 = sin((g1+g3) / 2.d0)
      sgg4 = sin((g1+g4) / 2.d0)
      so2=sin((a1-a2)/2.d0)
      so3=sin((a1-a3)/2.d0)
      so4=sin((a1-a4)/2.d0)
      co2=cos((a1-a2)/2.d0)
      co3=cos((a1-a3)/2.d0)
      co4=cos((a1-a4)/2.d0)
      cabn2=cos((a1-a2-b1-b2)/2.d0)
      cabn3=cos((a1-a3-b1-b3)/2.d0)
      cabn4=cos((a1-a4-b1-b4)/2.d0)
      caba2=cos((a1-a2+b1+b2)/2.d0)
      caba3=cos((a1-a3+b1+b3)/2.d0)
      caba4=cos((a1-a4+b1+b4)/2.d0)
      s2b1=sin(2.d0*b1)
      c2b1=cos(2.d0*b1)
      s2g1=sin(2.d0*g1)
      c2g1=cos(2.d0*g1)
      sab11=sin(a1-b1)
      cab11=cos(a1-b1)
      saba11=sin(a1+b1)
      caba11=cos(a1+b1)
c
c     相対変位ベクトル
      call clrv (24,dl)
c check!!  ok!
      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 check ok!
      call clrmx (24,R)
      R(7,1)=-(cb1*cg1)
      R(7,2)=-(cg1*sa1*sb1) - ca1*sg1
      R(7,3)=ca1*cg1*sb1 - sa1*sg1
      R(7,4)=sa1*(-((w1 - w2)*cg1*sb1) + (v1 - v2)*sg1)
     &+ ca1*(-((v1 - v2)*cg1*sb1) + (-w1 + w2)*sg1)
      R(7,5)=-(cg1*(-((w1 - w2)*ca1*cb1) + (v1 - v2)*cb1*sa1
     &+ a*sb1 - u1*sb1 + u2*sb1 - a*cg1*s2b1
     &+ a*ca1**2*cg1*s2b1 + a*cg1*sa1**2*s2b1))
      R(7,6)=-((w1 - w2)*cg1*sa1) - a*cb1*sg1 + u1*cb1*sg1 - u2*cb1*sg1
     &+ v1*sa1*sb1*sg1 - v2*sa1*sb1*sg1 + ca1*((-v1 + v2)*cg1
     &+ (-w1 + w2)*sb1*sg1) + a*cb1**2*s2g1
     &- a*ca1**2*cb1**2*s2g1 - a*sa1**2*s2g1 + a*sa1**2*sb1**2*s2g1
      R(8,1)=cb1*sg1
      R(8,2)=-(ca1*cg1) + sa1*sb1*sg1
      R(8,3)=-(cg1*sa1) - ca1*sb1*sg1
      R(8,4)=ca1*((-w1 + w2)*cg1 + (v1 - v2)*sb1*sg1)
     &+ sa1*((v1 - v2)*cg1 + (w1 - w2)*sb1*sg1)
      R(8,5)=(-((w1 - w2)*ca1*cb1) + (v1 - v2)*cb1*sa1
     &+ a*sb1 - u1*sb1 + u2*sb1 - a*cg1*s2b1 + a*ca1**2*cg1*s2b1
     &+ a*cg1*sa1**2*s2b1)*sg1
      R(8,6)=-((a - u1 + u2)*cb1*cg1) - a*cb1**2*cg1**2*sa1**2
     &+ a*cb1**2*c2g1*sa1**2 - cg1*(w1*ca1 - w2*ca1 - v1*sa1 + v2*sa1)
     &*sb1 + sg1*(v1*ca1 - v2*ca1 + w1*sa1 - w2*sa1 + a*sa1**2*sg1
     &- a*sa1**2*sb1**2*sg1)
      R(9,1)=-sb1
      R(9,2)=cb1*sa1
      R(9,3)=-(ca1*cb1)
      R(9,4)=cb1*((v1 - v2)*ca1 + (w1 - w2)*sa1)
      R(9,5)=((-v1 + v2)*cab11 + 2.d0*(a - u1 + u2)*cb1 + v1*caba11
     &- v2*caba11 - w1*sab11 + w2*sab11 + w1*saba11 - w2*saba11)/2.d0
      R(10,4)=-(cb1*cg1) + (((g1 - g2)*caa2*cbb2 + (b1 - b2)*saa2)
     &*(cg1*sa1*sb1 + ca1*sg1))/2.d0 + ((-g1 + g2)*caa2*cbb2
     &+ (-b1 + b2)*saa2)*(cg1*sa1*sb1 + ca1*sg1) + ((-b1 + b2)*caa2
     &+ (g1 - g2)*cbb2*saa2)*(ca1*cg1*sb1 - sa1*sg1) + (((b1 - b2)*caa2
     &+ (-g1 + g2)*cbb2*saa2)*(ca1*cg1*sb1 - sa1*sg1))/2.d0
      R(10,5)=((-g1 + g2)*cb1*cbb2*cg1)/2.d0 - ca1*cb1*cg1*((-g1 + g2)
     &*caa2*cbb2 + (-b1 + b2)*saa2) + cb1*cg1*sa1*((-b1 + b2)*caa2
     &+ (g1 - g2)*cbb2*saa2) - cg1*sb1*(-a1 + a2 + (-g1 + g2)*sbb2)
     &+ (-caa2 + ((-g1 + g2)*saa2*sbb2)/2.d0)*(cg1*sa1*sb1 + ca1*sg1)
     &+ ((2.d0*saa2 + (-g1 + g2)*caa2*sbb2)
     &*(ca1*cg1*sb1 - sa1*sg1))/2.d0
      R(10,6)=-(cb1*cg1*sbb2) - cb1*(-a1 + a2 + (-g1 + g2)*sbb2)*sg1
     &+ cbb2*saa2*(cg1*sa1*sb1 + ca1*sg1) + caa2*cbb2*(ca1*cg1*sb1
     &- sa1*sg1) + ((-g1 + g2)*caa2*cbb2 + (-b1 + b2)*saa2)*(cg1*sa1
     &+ ca1*sb1*sg1) + ((-b1 + b2)*caa2 + (g1 - g2)*cbb2*saa2)
     &*(ca1*cg1 - sa1*sb1*sg1)
      R(10,10)=(2.d0*cb1*cg1 + ((g1 - g2)*caa2*cbb2 + (b1 - b2)*saa2)
     &*(cg1*sa1*sb1 + ca1*sg1) + ((b1 - b2)*caa2
     &+ (-g1 + g2)*cbb2*saa2)*(ca1*cg1*sb1 - sa1*sg1))/2.d0
      R(10,11)=((-g1 + g2)*cb1*cbb2*cg1)/2.d0 + (caa2 + ((-g1 + g2)
     &*saa2*sbb2)/2.d0)*(cg1*sa1*sb1 + ca1*sg1) + (saa2 + ((g1 - g2)
     &*caa2*sbb2)/2.d0)*(-(ca1*cg1*sb1) + sa1*sg1)
      R(10,12)=cb1*cg1*sbb2 - cbb2*saa2*(cg1*sa1*sb1 + ca1*sg1)
     &+ caa2*cbb2*(-(ca1*cg1*sb1) + sa1*sg1)
      R(11,4)=cb1*sg1 + ((-b1 + b2)*caa2 + (g1 - g2)*cbb2*saa2)
     &*(-(cg1*sa1) - ca1*sb1*sg1) - (((b1 - b2)*caa2 + (-g1 + g2)
     &*cbb2*saa2)*(cg1*sa1 + ca1*sb1*sg1))/2.d0 + (((g1 - g2)*caa2*cbb2
     &+ (b1 - b2)*saa2)*(ca1*cg1 - sa1*sb1*sg1))/2.d0 + ((-g1 + g2)
     &*caa2*cbb2 + (-b1 + b2)*saa2)*(ca1*cg1 - sa1*sb1*sg1)
      R(11,5)=((g1 - g2)*cb1*cbb2*sg1)/2.d0 + ca1*cb1*((-g1 + g2)
     &*caa2*cbb2 + (-b1 + b2)*saa2)*sg1 - cb1*sa1*((-b1 + b2)*caa2
     &+ (g1 - g2)*cbb2*saa2)*sg1 + sb1*(-a1 + a2 + (-g1 + g2)*sbb2)*sg1
     &- ((2.d0*saa2 + (-g1 + g2)*caa2*sbb2)*(cg1*sa1 + ca1*sb1*sg1))
     &/2.d0 + (-caa2 + ((-g1 + g2)*saa2*sbb2)/2.d0)
     &*(ca1*cg1 - sa1*sb1*sg1)
      R(11,6)=-(cb1*cg1*(-a1 + a2 + (-g1 + g2)*sbb2)) + cb1*sbb2*sg1
     &+ ((-b1 + b2)*caa2 + (g1 - g2)*cbb2*saa2)*(-(cg1*sa1*sb1)
     &- ca1*sg1) + ((-g1 + g2)*caa2*cbb2 + (-b1 + b2)*saa2)*(ca1*cg1*sb1
     &- sa1*sg1) - caa2*cbb2*(cg1*sa1 + ca1*sb1*sg1)
     &+ cbb2*saa2*(ca1*cg1 - sa1*sb1*sg1)
      R(11,10)=(-2.d0*cb1*sg1 - ((b1 - b2)*caa2 + (-g1 + g2)*cbb2*saa2)
     &*(cg1*sa1 + ca1*sb1*sg1) + ((g1 - g2)*caa2*cbb2
     &+ (b1 - b2)*saa2)*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      R(11,11)=((g1 - g2)*cb1*cbb2*sg1)/2.d0 + (saa2 + ((g1 - g2)
     &*caa2*sbb2)/2.d0)*(cg1*sa1 + ca1*sb1*sg1)
     &+ (caa2 + ((-g1 + g2)*saa2*sbb2)/2.d0)*(ca1*cg1 - sa1*sb1*sg1)
      R(11,12)=-(cb1*sbb2*sg1) + caa2*cbb2*(cg1*sa1 + ca1*sb1*sg1)
     &- cbb2*saa2*(ca1*cg1 - sa1*sb1*sg1)
      R(12,4)=((g1 - g2)*caa2*cb1*cbb2*sa1 + b1*cb1*sa1*saa2
     &- b2*cb1*sa1*saa2 + ca1*cb1*((b1 - b2)*caa2
     &+ (-g1 + g2)*cbb2*saa2) - 2.d0*sb1)/2.d0
      R(12,5)=(((g1 - g2)*(-1.d0 + 2.d0*co2)*cbb2 + 2.d0*(-b1 + b2)*so2)
     &*sb1 + cb1*(-2.d0*a1 + 2.d0*a2 - 2.d0*ca1*saa2 - 2.d0*g1*sbb2
     &+ 2.d0*g2*sbb2 + g1*sa1*saa2*sbb2 - g2*sa1*saa2*sbb2
     &+ caa2*(2.d0*sa1 + g1*ca1*sbb2 - g2*ca1*sbb2)))/2.d0
      R(12,6)=-(ca1*caa2*cb1*cbb2) - cb1*cbb2*sa1*saa2 - sb1*sbb2
      R(12,10)=-(cb1*sa1*((g1 - g2)*caa2*cbb2 + (b1 - b2)*saa2))/2.d0
     &- (ca1*cb1*((b1 - b2)*caa2 + (-g1 + g2)*cbb2*saa2))/2.d0 + sb1
      R(12,11)=(2.d0*ca1*cb1*saa2 + caa2*cb1*(-2.d0*sa1 + (g1 - g2)
     &*ca1*sbb2) - (g1 - g2)*(cbb2*sb1 - cb1*sa1*saa2*sbb2))/2.d0
      R(12,12)=ca1*caa2*cb1*cbb2 + cb1*cbb2*sa1*saa2 + sb1*sbb2
      R(13,4)=-(sa1*((w1 - w3)*cg1*sb1 + (b - v1 + v3)*sg1))
     &+ ca1*((b - v1 + v3)*cg1*sb1 + (-w1 + w3)*sg1)
      R(13,5)=(cg1*((w1 - w3)*cab11 + (w1 - w3)*caba11 + b*sab11
     &- v1*sab11 + v3*sab11 - 2.d0*a*sb1 + 2.d0*u1*sb1 - 2.d0*u3*sb1
     &+ b*saba11 - v1*saba11 + v3*saba11))/2.d0
      R(13,6)=-((a - u1 + u3)*cb1*sg1) + ca1*((b - v1 + v3)*cg1
     &+ (-w1 + w3)*sb1*sg1) + cb1**2*(b*c2g1 + a*s2g1) - ca1**2*cb1**2
     &*(b*c2g1 + a*s2g1) + sa1*(-(w1*cg1) + w3*cg1 - b*cg1**2*sa1
     &+ b*cg1**2*sa1*sb1**2 - b*sb1*sg1 + v1*sb1*sg1 - v3*sb1*sg1
     &+ b*sa1*sg1**2 - b*sa1*sb1**2*sg1**2 - a*sa1*s2g1
     &+ a*sa1*sb1**2*s2g1)
      R(14,4)=-(ca1*((w1 - w3)*cg1 + (b - v1 + v3)*sb1*sg1))
     &- sa1*((b - v1 + v3)*cg1 + (-w1 + w3)*sb1*sg1)
      R(14,5)=-(((w1 - w3)*cab11 + (w1 - w3)*caba11 + b*sab11
     &- v1*sab11 + v3*sab11 - 2.d0*a*sb1 + 2.d0*u1*sb1 - 2.d0*u3*sb1
     &+ b*saba11 - v1*saba11 + v3*saba11)*sg1)/2.d0
      R(14,6)=-((a - u1 + u3)*cb1*cg1) - w1*ca1*cg1*sb1 + w3*ca1*cg1*sb1
     &- b*cg1*sa1*sb1 + v1*cg1*sa1*sb1 - v3*cg1*sa1*sb1 - b*ca1*sg1
     &+ v1*ca1*sg1 - v3*ca1*sg1 + w1*sa1*sg1 - w3*sa1*sg1
     &+ a*sa1**2*sg1**2 - a*sa1**2*sb1**2*sg1**2 - cb1**2*sa1**2*sg1
     &*(2.d0*b*cg1 + a*sg1) + b*sa1**2*s2g1 - b*sa1**2*sb1**2*s2g1
      R(15,4)=-(cb1*((b - v1 + v3)*ca1 + (-w1 + w3)*sa1))
      R(15,5)=((b - v1 + v3)*cab11 + 2.d0*(a - u1 + u3)*cb1 - b*caba11
     &+ v1*caba11 - v3*caba11 - w1*sab11 + w3*sab11
     &+ w1*saba11 - w3*saba11)/2.d0
      R(16,4)=-(cb1*cg1) + (((g1 - g3)*caa3*cbb3 + (b1 - b3)*saa3)
     &*(cg1*sa1*sb1 + ca1*sg1))/2.d0 + ((-g1 + g3)*caa3*cbb3
     &+ (-b1 + b3)*saa3)*(cg1*sa1*sb1 + ca1*sg1) + ((-b1 + b3)*caa3
     &+ (g1 - g3)*cbb3*saa3)*(ca1*cg1*sb1 - sa1*sg1) + (((b1 - b3)*caa3
     &+ (-g1 + g3)*cbb3*saa3)*(ca1*cg1*sb1 - sa1*sg1))/2.d0
      R(16,5)=((-g1 + g3)*cb1*cbb3*cg1)/2.d0 - ca1*cb1*cg1*((-g1 + g3)
     &*caa3*cbb3 + (-b1 + b3)*saa3) + cb1*cg1*sa1*((-b1 + b3)*caa3
     &+ (g1 - g3)*cbb3*saa3) - cg1*sb1*(-a1 + a3 + (-g1 + g3)*sbb3)
     &+ (-caa3 + ((-g1 + g3)*saa3*sbb3)/2.d0)*(cg1*sa1*sb1 + ca1*sg1)
     &+ ((2.d0*saa3 + (-g1 + g3)*caa3*sbb3)*(ca1*cg1*sb1 - sa1*sg1))
     &/2.d0
      R(16,6)=-(cb1*cg1*sbb3) - cb1*(-a1 + a3 + (-g1 + g3)*sbb3)*sg1
     &+ cbb3*saa3*(cg1*sa1*sb1 + ca1*sg1) + caa3*cbb3*(ca1*cg1*sb1
     &- sa1*sg1) + ((-g1 + g3)*caa3*cbb3 + (-b1 + b3)*saa3)*(cg1*sa1
     &+ ca1*sb1*sg1) + ((-b1 + b3)*caa3 + (g1 - g3)*cbb3*saa3)
     &*(ca1*cg1 - sa1*sb1*sg1)
      R(16,16)=(2.d0*cb1*cg1 + ((g1 - g3)*caa3*cbb3 + (b1 - b3)*saa3)
     &*(cg1*sa1*sb1 + ca1*sg1) + ((b1 - b3)*caa3 + (-g1 + g3)
     &*cbb3*saa3)*(ca1*cg1*sb1 - sa1*sg1))/2.d0
      R(16,17)=((-g1 + g3)*cb1*cbb3*cg1)/2.d0 + (caa3 + ((-g1 + g3)
     &*saa3*sbb3)/2.d0)*(cg1*sa1*sb1 + ca1*sg1) + (saa3 + ((g1 - g3)
     &*caa3*sbb3)/2.d0)*(-(ca1*cg1*sb1) + sa1*sg1)
      R(16,18)=cb1*cg1*sbb3 - cbb3*saa3*(cg1*sa1*sb1 + ca1*sg1)
     &+ caa3*cbb3*(-(ca1*cg1*sb1) + sa1*sg1)
      R(17,4)=cb1*sg1 + ((-b1 + b3)*caa3 + (g1 - g3)*cbb3*saa3)
     &*(-(cg1*sa1) - ca1*sb1*sg1) - (((b1 - b3)*caa3 + (-g1 + g3)
     &*cbb3*saa3)*(cg1*sa1 + ca1*sb1*sg1))/2.d0 + (((g1 - g3)*caa3*cbb3
     &+ (b1 - b3)*saa3)*(ca1*cg1 - sa1*sb1*sg1))/2.d0 + ((-g1 + g3)
     &*caa3*cbb3 + (-b1 + b3)*saa3)*(ca1*cg1 - sa1*sb1*sg1)
      R(17,5)=((g1 - g3)*cb1*cbb3*sg1)/2.d0 + ca1*cb1*((-g1 + g3)
     &*caa3*cbb3 + (-b1 + b3)*saa3)*sg1 - cb1*sa1*((-b1 + b3)*caa3
     &+ (g1 - g3)*cbb3*saa3)*sg1 + sb1*(-a1 + a3 + (-g1 + g3)*sbb3)*sg1
     &- ((2.d0*saa3 + (-g1 + g3)*caa3*sbb3)*(cg1*sa1 + ca1*sb1*sg1))
     &/2.d0 + (-caa3 + ((-g1 + g3)*saa3*sbb3)/2.d0)
     &*(ca1*cg1 - sa1*sb1*sg1)
      R(17,6)=-(cb1*cg1*(-a1 + a3 + (-g1 + g3)*sbb3)) + cb1*sbb3*sg1
     &+ ((-b1 + b3)*caa3 + (g1 - g3)*cbb3*saa3)*(-(cg1*sa1*sb1)
     &- ca1*sg1) + ((-g1 + g3)*caa3*cbb3 + (-b1 + b3)*saa3)
     &*(ca1*cg1*sb1 - sa1*sg1) - caa3*cbb3*(cg1*sa1 + ca1*sb1*sg1)
     &+ cbb3*saa3*(ca1*cg1 - sa1*sb1*sg1)
      R(17,16)=(-2.d0*cb1*sg1 - ((b1 - b3)*caa3 + (-g1 + g3)*cbb3*saa3)
     &*(cg1*sa1 + ca1*sb1*sg1) + ((g1 - g3)*caa3*cbb3
     &+ (b1 - b3)*saa3)*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      R(17,17)=((g1 - g3)*cb1*cbb3*sg1)/2.d0 + (saa3 + ((g1 - g3)
     &*caa3*sbb3)/2.d0)*(cg1*sa1 + ca1*sb1*sg1)
     &+ (caa3 + ((-g1 + g3)*saa3*sbb3)/2.d0)*(ca1*cg1 - sa1*sb1*sg1)
      R(17,18)=-(cb1*sbb3*sg1) + caa3*cbb3*(cg1*sa1 + ca1*sb1*sg1)
     &- cbb3*saa3*(ca1*cg1 - sa1*sb1*sg1)
      R(18,4)=((g1 - g3)*caa3*cb1*cbb3*sa1 + b1*cb1*sa1*saa3
     &- b3*cb1*sa1*saa3 + ca1*cb1*((b1 - b3)*caa3 + (-g1 + g3)
     &*cbb3*saa3) - 2.d0*sb1)/2.d0
      R(18,5)=(((g1 - g3)*(-1.d0 + 2.d0*co3)*cbb3 + 2.d0*(-b1 + b3)*so3)
     &*sb1 + cb1*(-2.d0*a1 + 2.d0*a3 - 2.d0*ca1*saa3 - 2.d0*g1*sbb3
     &+ 2.d0*g3*sbb3 + g1*sa1*saa3*sbb3 - g3*sa1*saa3*sbb3
     &+ caa3*(2.d0*sa1 + g1*ca1*sbb3 - g3*ca1*sbb3)))/2.d0
      R(18,6)=-(ca1*caa3*cb1*cbb3) - cb1*cbb3*sa1*saa3 - sb1*sbb3
      R(18,16)=-(cb1*sa1*((g1 - g3)*caa3*cbb3 + (b1 - b3)*saa3))/2.d0
     &- (ca1*cb1*((b1 - b3)*caa3 + (-g1 + g3)*cbb3*saa3))/2.d0 + sb1
      R(18,17)=(2.d0*ca1*cb1*saa3 + caa3*cb1*(-2.d0*sa1 + (g1 - g3)*ca1
     &*sbb3) - (g1 - g3)*(cbb3*sb1 - cb1*sa1*saa3*sbb3))/2.d0
      R(18,18)=ca1*caa3*cb1*cbb3 + cb1*cbb3*sa1*saa3 + sb1*sbb3
      R(19,4)=-(sa1*((w1 - w4)*cg1*sb1 + (b - v1 + v4)*sg1))
     &+ ca1*((b - v1 + v4)*cg1*sb1 + (-w1 + w4)*sg1)
      R(19,5)=cg1*((w1 - w4)*ca1*cb1 + (b - v1 + v4)*cb1*sa1 + u1*sb1
     &- u4*sb1 - b*s2b1*sg1 + b*ca1**2*s2b1*sg1 + b*sa1**2*s2b1*sg1)
      R(19,6)=b*cb1**2*cg1**2 + b*(cg1*sa1*sb1 + ca1*sg1)**2 - cb1*sg1
     &*(-u1 + u4 + b*cb1*sg1) + b*(ca1*cg1*sb1 - sa1*sg1)**2 - (cg1*sa1
     &+ ca1*sb1*sg1)*(w1 - w4 + b*cg1*sa1 + b*ca1*sb1*sg1) + (ca1*cg1
     &- sa1*sb1*sg1)*(b - v1 + v4 - b*ca1*cg1 + b*sa1*sb1*sg1)
      R(20,4)=-(ca1*((w1 - w4)*cg1 + (b - v1 + v4)*sb1*sg1))
     &- sa1*((b - v1 + v4)*cg1 + (-w1 + w4)*sb1*sg1)
      R(20,5)=-(sg1*((w1 - w4)*ca1*cb1 + (b - v1 + v4)*cb1*sa1 + u1*sb1
     &- u4*sb1 - b*s2b1*sg1 + b*ca1**2*s2b1*sg1 + b*sa1**2*s2b1*sg1))
      R(20,6)=(u1 - u4)*cb1*cg1 - ca1*((w1 - w4)*cg1*sb1 + (b - v1 + v4)
     &*sg1) - b*cb1**2*sa1**2*s2g1 - sa1*(b*cg1*sb1 - v1*cg1*sb1
     &+ v4*cg1*sb1 - w1*sg1 + w4*sg1 - b*sa1*s2g1 + b*sa1*sb1**2*s2g1)
      R(21,4)=-(cb1*((b - v1 + v4)*ca1 + (-w1 + w4)*sa1))
      R(21,5)=((b - v1 + v4)*cab11 - 2.d0*(u1 - u4)*cb1 - b*caba11
     &+ v1*caba11 - v4*caba11 - w1*sab11 + w4*sab11
     &+ w1*saba11 - w4*saba11)/2.d0
      R(22,4)=-(cb1*cg1) + (((g1 - g4)*caa4*cbb4 + (b1 - b4)*saa4)
     &*(cg1*sa1*sb1 + ca1*sg1))/2.d0 + ((-g1 + g4)*caa4*cbb4
     &+ (-b1 + b4)*saa4)*(cg1*sa1*sb1 + ca1*sg1) + ((-b1 + b4)*caa4
     &+ (g1 - g4)*cbb4*saa4)*(ca1*cg1*sb1 - sa1*sg1) + (((b1 - b4)*caa4
     &+ (-g1 + g4)*cbb4*saa4)*(ca1*cg1*sb1 - sa1*sg1))/2.d0
      R(22,5)=((-g1 + g4)*cb1*cbb4*cg1)/2.d0 - ca1*cb1*cg1*((-g1 + g4)
     &*caa4*cbb4 + (-b1 + b4)*saa4) + cb1*cg1*sa1*((-b1 + b4)*caa4
     &+ (g1 - g4)*cbb4*saa4) - cg1*sb1*(-a1 + a4 + (-g1 + g4)*sbb4)
     &+ (-caa4 + ((-g1 + g4)*saa4*sbb4)/2.d0)*(cg1*sa1*sb1 + ca1*sg1)
     &+ ((2.d0*saa4 + (-g1 + g4)*caa4*sbb4)
     &*(ca1*cg1*sb1 - sa1*sg1))/2.d0
      R(22,6)=-(cb1*cg1*sbb4) - cb1*(-a1 + a4 + (-g1 + g4)*sbb4)*sg1
     &+ cbb4*saa4*(cg1*sa1*sb1 + ca1*sg1) + caa4*cbb4*(ca1*cg1*sb1
     &- sa1*sg1) + ((-g1 + g4)*caa4*cbb4 + (-b1 + b4)*saa4)*(cg1*sa1
     &+ ca1*sb1*sg1) + ((-b1 + b4)*caa4
     &+ (g1 - g4)*cbb4*saa4)*(ca1*cg1 - sa1*sb1*sg1)
      R(22,22)=(2.d0*cb1*cg1 + ((g1 - g4)*caa4*cbb4 + (b1 - b4)*saa4)
     &*(cg1*sa1*sb1 + ca1*sg1) + ((b1 - b4)*caa4 + (-g1 + g4)*cbb4*saa4)
     &*(ca1*cg1*sb1 - sa1*sg1))/2.d0
      R(22,23)=((-g1 + g4)*cb1*cbb4*cg1)/2.d0 + (caa4 + ((-g1 + g4)
     &*saa4*sbb4)/2.d0)*(cg1*sa1*sb1 + ca1*sg1) + (saa4 + ((g1 - g4)
     &*caa4*sbb4)/2.d0)*(-(ca1*cg1*sb1) + sa1*sg1)
      R(22,24)=cb1*cg1*sbb4 - cbb4*saa4*(cg1*sa1*sb1 + ca1*sg1)
     &+ caa4*cbb4*(-(ca1*cg1*sb1) + sa1*sg1)
      R(23,4)=cb1*sg1 + ((-b1 + b4)*caa4 + (g1 - g4)*cbb4*saa4)
     &*(-(cg1*sa1) - ca1*sb1*sg1) - (((b1 - b4)*caa4 + (-g1 + g4)
     &*cbb4*saa4)*(cg1*sa1 + ca1*sb1*sg1))/2.d0 + (((g1 - g4)*caa4*cbb4
     &+ (b1 - b4)*saa4)*(ca1*cg1 - sa1*sb1*sg1))/2.d0 + ((-g1 + g4)
     &*caa4*cbb4 + (-b1 + b4)*saa4)*(ca1*cg1 - sa1*sb1*sg1)
      R(23,5)=((g1 - g4)*cb1*cbb4*sg1)/2.d0 + ca1*cb1*((-g1 + g4)
     &*caa4*cbb4 + (-b1 + b4)*saa4)*sg1 - cb1*sa1*((-b1 + b4)*caa4
     &+ (g1 - g4)*cbb4*saa4)*sg1 + sb1*(-a1 + a4 + (-g1 + g4)*sbb4)*sg1
     &- ((2.d0*saa4 + (-g1 + g4)*caa4*sbb4)*(cg1*sa1 + ca1*sb1*sg1))
     &/2.d0 + (-caa4 + ((-g1 + g4)*saa4*sbb4)/2.d0)
     &*(ca1*cg1 - sa1*sb1*sg1)
      R(23,6)=-(cb1*cg1*(-a1 + a4 + (-g1 + g4)*sbb4)) + cb1*sbb4*sg1
     &+ ((-b1 + b4)*caa4 + (g1 - g4)*cbb4*saa4)*(-(cg1*sa1*sb1)
     &- ca1*sg1) + ((-g1 + g4)*caa4*cbb4 + (-b1 + b4)*saa4)*(ca1*cg1*sb1
     &- sa1*sg1) - caa4*cbb4*(cg1*sa1 + ca1*sb1*sg1)
     &+ cbb4*saa4*(ca1*cg1 - sa1*sb1*sg1)
      R(23,22)=(-2.d0*cb1*sg1 - ((b1 - b4)*caa4 + (-g1 + g4)*cbb4*saa4)
     &*(cg1*sa1 + ca1*sb1*sg1) + ((g1 - g4)*caa4*cbb4 + (b1 - b4)*saa4)
     &*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      R(23,23)=((g1 - g4)*cb1*cbb4*sg1)/2.d0 + (saa4 + ((g1 - g4)*caa4
     &*sbb4)/2.d0)*(cg1*sa1 + ca1*sb1*sg1) + (caa4 + ((-g1 + g4)
     &*saa4*sbb4)/2.d0)*(ca1*cg1 - sa1*sb1*sg1)
      R(23,24)=-(cb1*sbb4*sg1) + caa4*cbb4*(cg1*sa1 + ca1*sb1*sg1)
     &- cbb4*saa4*(ca1*cg1 - sa1*sb1*sg1)
      R(24,4)=((g1 - g4)*caa4*cb1*cbb4*sa1 + b1*cb1*sa1*saa4
     &- b4*cb1*sa1*saa4 + ca1*cb1*((b1 - b4)*caa4 + (-g1 + g4)
     &*cbb4*saa4) - 2.d0*sb1)/2.d0
      R(24,5)=(((g1 - g4)*(-1.d0 + 2.d0*co4)*cbb4 + 2.d0*(-b1 + b4)
     &*so4)*sb1 + cb1*(-2.d0*a1 + 2.d0*a4 - 2.d0*ca1*saa4
     &- 2.d0*g1*sbb4 + 2.d0*g4*sbb4 + g1*sa1*saa4*sbb4
     &- g4*sa1*saa4*sbb4 + caa4*(2.d0*sa1 + g1*ca1*sbb4
     &- g4*ca1*sbb4)))/2.d0
      R(24,6)=-(ca1*caa4*cb1*cbb4) - cb1*cbb4*sa1*saa4 - sb1*sbb4
      R(24,22)=-(cb1*sa1*((g1 - g4)*caa4*cbb4 + (b1 - b4)*saa4))/2.d0
     &- (ca1*cb1*((b1 - b4)*caa4 + (-g1 + g4)*cbb4*saa4))/2.d0 + sb1
      R(24,23)=(2.d0*ca1*cb1*saa4 + caa4*cb1*(-2.d0*sa1 + (g1 - g4)
     &*ca1*sbb4) - (g1 - g4)*(cbb4*sb1 - cb1*sa1*saa4*sbb4))/2.d0
      R(24,24)=ca1*caa4*cb1*cbb4 + cb1*cbb4*sa1*saa4 + sb1*sbb4
      do 801 i=1,3
         do 801 j=1,3
            R(12+i,j) = R(6+i,j)
            R(18+i,j) = R(6+i,j)
            R(6+i,6+j) = -R(6+i,j)
            R(12+i,12+j) = -R(6+i,j)
            R(18+i,18+j) = -R(6+i,j)
  801 continue
c
      call clmx3 (24,dte)
      call clrmx(3,ste)
      do 13 k=1,2
         call dseuler(ste,a1,b1,k)
         do 13 i=1,3
            do 13 j=1,3
               dte(k+3,i+3,j+3) = ste(i,j)
   13 continue
      do 23 k=1,2
         call dseuler(ste,a2,b2,k)
         do 23 i=1,3
            do 23 j=1,3
               dte(k+3,i+9,j+9) = ste(i,j)
   23 continue
      do 33 k=1,2
         call dseuler(ste,a3,b3,k)
         do 33 i=1,3
            do 33 j=1,3
               dte(k+3,i+15,j+15) = ste(i,j)
   33 continue
      do 43 k=1,2
         call dseuler(ste,a4,b4,k)
         do 43 i=1,3
            do 43 j=1,3
               dte(k+3,i+21,j+21) = ste(i,j)
   43 continue
c
      call clmx3 (24,dR)
      dR(1,7,5)=cg1*sb1
      dR(1,7,6)=cb1*sg1
      dR(1,8,5)=-(sb1*sg1)
      dR(1,8,6)=cb1*cg1
      dR(1,9,5)=-cb1
      do 901 i=1,3
         do 901 j=1,3
            dR(1,12+i,3+j) = dR(1,i+6,3+j)
            dR(1,18+i,3+j) = dR(1,i+6,3+j)
  901 continue
c
      dR(2,7,4)=-(ca1*cg1*sb1) + sa1*sg1
      dR(2,7,5)=-(cb1*cg1*sa1)
      dR(2,7,6)=-(ca1*cg1) + sa1*sb1*sg1
      dR(2,8,4)=cg1*sa1 + ca1*sb1*sg1
      dR(2,8,5)=cb1*sa1*sg1
      dR(2,8,6)=cg1*sa1*sb1 + ca1*sg1
      dR(2,9,4)=ca1*cb1
      dR(2,9,5)=-(sa1*sb1)
      do 902 i=1,3
         do 902 j=1,3
            dR(2,12+i,3+j) = dR(2,i+6,3+j)
            dR(2,18+i,3+j) = dR(2,i+6,3+j)
  902 continue
c
      dR(3,7,4)=-(cg1*sa1*sb1) - ca1*sg1
      dR(3,7,5)=ca1*cb1*cg1
      dR(3,7,6)=-(cg1*sa1) - ca1*sb1*sg1
      dR(3,8,4)=-(ca1*cg1) + sa1*sb1*sg1
      dR(3,8,5)=-(ca1*cb1*sg1)
      dR(3,8,6)=-(ca1*cg1*sb1) + sa1*sg1
      dR(3,9,4)=cb1*sa1
      dR(3,9,5)=ca1*sb1
      do 903 i=1,3
         do 903 j=1,3
            dR(3,12+i,3+j) = dR(3,i+6,3+j)
            dR(3,18+i,3+j) = dR(3,i+6,3+j)
  903 continue
c
      dR(4,7,2)=-(ca1*cg1*sb1) + sa1*sg1
      dR(4,7,3)=-(cg1*sa1*sb1) - ca1*sg1
      dR(4,8,2)=cg1*sa1 + ca1*sb1*sg1
      dR(4,8,3)=-(ca1*cg1) + sa1*sb1*sg1
      dR(4,9,2)=ca1*cb1
      dR(4,9,3)=cb1*sa1
c
      dR(4,7,4)=ca1*(-((w1 - w2)*cg1*sb1) + (v1 - v2)*sg1) +
     &sa1*((v1 - v2)*cg1*sb1 + (w1 - w2)*sg1)
      dR(4,7,5)=cb1*cg1*((-v1 + v2)*ca1 + (-w1 + w2)*sa1)
      dR(4,7,6)=ca1*((-w1 + w2)*cg1 + (v1 - v2)*sb1*sg1)
     &+ sa1*((v1 - v2)*cg1 + (w1 - w2)*sb1*sg1)
      dR(4,8,4)=sa1*((w1 - w2)*cg1 + (-v1 + v2)*sb1*sg1)
     &+ ca1*((v1 - v2)*cg1 + (w1 - w2)*sb1*sg1)
      dR(4,8,5)=cb1*((v1 - v2)*ca1 + (w1 - w2)*sa1)*sg1
      dR(4,8,6)=sa1*((w1 - w2)*cg1*sb1 + (-v1 + v2)*sg1)
     &+ ca1*((v1 - v2)*cg1*sb1 + (w1 - w2)*sg1)
      dR(4,9,4)=cb1*((w1 - w2)*ca1 + (-v1 + v2)*sa1)
      dR(4,9,5)=((-v1 + v2)*ca1 + (-w1 + w2)*sa1)*sb1
c
      dR(4,13,4)=-(ca1*((w1 - w3)*cg1*sb1 + (b - v1 + v3)*sg1))
     &- sa1*((b - v1 + v3)*cg1*sb1 + (-w1 + w3)*sg1)
      dR(4,13,5)=cb1*cg1*((b - v1 + v3)*ca1 + (-w1 + w3)*sa1)
      dR(4,13,6)=-(ca1*((w1 - w3)*cg1 + (b - v1 + v3)*sb1*sg1))
     &- sa1*((b - v1 + v3)*cg1 + (-w1 + w3)*sb1*sg1)
      dR(4,14,4)=sa1*((w1 - w3)*cg1 + (b - v1 + v3)*sb1*sg1)
     &- ca1*((b - v1 + v3)*cg1 + (-w1 + w3)*sb1*sg1)
      dR(4,14,5)=-(cb1*((b - v1 + v3)*ca1 + (-w1 + w3)*sa1)*sg1)
      dR(4,14,6)=sa1*((w1 - w3)*cg1*sb1 + (b - v1 + v3)*sg1)
     &- ca1*((b - v1 + v3)*cg1*sb1 + (-w1 + w3)*sg1)
      dR(4,15,4)=cb1*((w1 - w3)*ca1 + (b - v1 + v3)*sa1)
      dR(4,15,5)=((b - v1 + v3)*ca1 + (-w1 + w3)*sa1)*sb1
c
      dR(4,19,4)=-(ca1*((w1 - w4)*cg1*sb1 + (b - v1 + v4)*sg1))
     &- sa1*((b - v1 + v4)*cg1*sb1 + (-w1 + w4)*sg1)
      dR(4,19,5)=cb1*cg1*((b - v1 + v4)*ca1 + (-w1 + w4)*sa1)
      dR(4,19,6)=-(ca1*((w1 - w4)*cg1 + (b - v1 + v4)*sb1*sg1))
     &- sa1*((b - v1 + v4)*cg1 + (-w1 + w4)*sb1*sg1)
      dR(4,20,4)=sa1*((w1 - w4)*cg1 + (b - v1 + v4)*sb1*sg1)
     &- ca1*((b - v1 + v4)*cg1 + (-w1 + w4)*sb1*sg1)
      dR(4,20,5)=-(cb1*((b - v1 + v4)*ca1 + (-w1 + w4)*sa1)*sg1)
      dR(4,20,6)=sa1*((w1 - w4)*cg1*sb1 + (b - v1 + v4)*sg1)
     &- ca1*((b - v1 + v4)*cg1*sb1 + (-w1 + w4)*sg1)
      dR(4,21,4)=cb1*((w1 - w4)*ca1 + (b - v1 + v4)*sa1)
      dR(4,21,5)=((b - v1 + v4)*ca1 + (-w1 + w4)*sa1)*sb1
c
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr4 (sa1,sb1,sg1,ca1,cb1,cg1,saa2,sbb2,caa2,cbb2,
     &         so2,co2,cabn2,caba2,b1,g1,b2,g2)
      do 904 i=1,3
         do 904 j=1,3
            dR(4,9+i,3+j) = dR1(i,j)
            dR(4,9+i,9+j) = dR2(i,j)
  904 continue
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr4 (sa1,sb1,sg1,ca1,cb1,cg1,saa3,sbb3,caa3,cbb3,
     &         so3,co3,cabn3,caba3,b1,g1,b3,g3)
      do 905 i=1,3
         do 905 j=1,3
            dR(4,15+i,3+j) = dR1(i,j)
            dR(4,15+i,15+j) = dR2(i,j)
  905 continue
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr4 (sa1,sb1,sg1,ca1,cb1,cg1,saa4,sbb4,caa4,cbb4,
     &         so4,co4,cabn4,caba4,b1,g1,b4,g4)
      do 906 i=1,3
         do 906 j=1,3
            dR(4,21+i,3+j) = dR1(i,j)
            dR(4,21+i,21+j) = dR2(i,j)
  906 continue
c
      do 907 i=1,3
         do 907 j=1,3
            dR(4,12+i,j) = dR(4,6+i,j)
            dR(4,18+i,j) = dR(4,6+i,j)
            dR(4,6+i,6+j) = -dR(4,6+i,j)
            dR(4,12+i,12+j) = -dR(4,6+i,j)
            dR(4,18+i,18+j) = -dR(4,6+i,j)
  907 continue
c
      dR(5,7,1)=cg1*sb1
      dR(5,7,2)=-(cb1*cg1*sa1)
      dR(5,7,3)=ca1*cb1*cg1
      dR(5,8,1)=-(sb1*sg1)
      dR(5,8,2)=cb1*sa1*sg1
      dR(5,8,3)=-(ca1*cb1*sg1)
      dR(5,9,1)=-cb1
      dR(5,9,2)=-(sa1*sb1)
      dR(5,9,3)=ca1*sb1
c
      dR(5,7,4)=cb1*cg1*((-v1 + v2)*ca1 + (-w1 + w2)*sa1)
      dR(5,7,5)=-(cg1*((-v1 + v2)*cab11 + 2.d0*(a - u1 + u2)*cb1
     &+ v1*caba11 - v2*caba11 - w1*sab11 + w2*sab11 + w1*saba11
     &- w2*saba11))/2.d0
      dR(5,7,6)=(-((w1 - w2)*ca1*cb1) + (a - u1 + u2)*sb1 + cb1*(v1*sa1
     &- v2*sa1 - 4.d0*a*cg1*sb1 + 4.d0*a*cg1*sa1**2*sb1)
     &+ 2.d0*a*ca1**2*cg1*s2b1)*sg1
      dR(5,8,4)=cb1*((v1 - v2)*ca1 + (w1 - w2)*sa1)*sg1
      dR(5,8,5)=(((-v1 + v2)*cab11 + 2.d0*(a - u1 + u2)*cb1 + v1*caba11
     &- v2*caba11 - w1*sab11 + w2*sab11 + w1*saba11
     &- w2*saba11)*sg1)/2.d0
      dR(5,8,6)=-(cg1*((w1 - w2)*cab11 + (w1 - w2)*caba11
     &- v1*sab11 + v2*sab11 - 2.d0*a*sb1 + 2.d0*u1*sb1 - 2.d0*u2*sb1
     &- v1*saba11 + v2*saba11))/2.d0
      dR(5,9,4)=((-v1 + v2)*ca1 + (-w1 + w2)*sa1)*sb1
      dR(5,9,5)=(w1 - w2)*ca1*cb1 - (a - u1 + u2)*sb1 + cb1*(-(v1*sa1)
     &+ v2*sa1 + 4.d0*a*cg1*sb1 - 4.d0*a*cg1*sa1**2*sb1)
     &- 2.d0*a*ca1**2*cg1*s2b1
c
      dR(5,13,4)=cb1*cg1*((b - v1 + v3)*ca1 + (-w1 + w3)*sa1)
      dR(5,13,5)=(cg1*(-((b - v1 + v3)*cab11) - 2.d0*(a - u1 + u3)*cb1
     &+ b*caba11 - v1*caba11 + v3*caba11 + w1*sab11 - w3*sab11
     &- w1*saba11 + w3*saba11))/2.d0
      dR(5,13,6)=-(((w1 - w3)*cab11 + (w1 - w3)*caba11 + b*sab11
     &- v1*sab11 + v3*sab11 - 2.d0*a*sb1 + 2.d0*u1*sb1 - 2.d0*u3*sb1
     &+ b*saba11 - v1*saba11 + v3*saba11)*sg1)/2.d0
      dR(5,14,4)=-(cb1*((b - v1 + v3)*ca1 + (-w1 + w3)*sa1)*sg1)
      dR(5,14,5)=(((b - v1 + v3)*cab11 + 2.d0*(a - u1 + u3)*cb1
     &- b*caba11 + v1*caba11 - v3*caba11 - w1*sab11 + w3*sab11
     &+ w1*saba11 - w3*saba11)*sg1)/2
      dR(5,14,6)=-(cg1*((w1 - w3)*cab11 + (w1 - w3)*caba11 + b*sab11
     &- v1*sab11 + v3*sab11 - 2.d0*a*sb1 + 2.d0*u1*sb1 - 2.d0*u3*sb1
     &+ b*saba11 - v1*saba11 + v3*saba11))/2.d0
      dR(5,15,4)=((b - v1 + v3)*ca1 + (-w1 + w3)*sa1)*sb1
      dR(5,15,5)=((w1 - w3)*cab11 + (w1 - w3)*caba11 + b*sab11
     &- v1*sab11 + v3*sab11 - 2.d0*a*sb1 + 2.d0*u1*sb1 - 2.d0*u3*sb1
     &+ b*saba11 - v1*saba11 + v3*saba11)/2.d0
c
      dR(5,19,4)=cb1*cg1*((b - v1 + v4)*ca1 + (-w1 + w4)*sa1)
      dR(5,19,5)=(cg1*(-((b - v1 + v4)*cab11) + 2.d0*(u1 - u4)*cb1
     &+ b*caba11 - v1*caba11 + v4*caba11 + w1*sab11 - w4*sab11
     &- w1*saba11 + w4*saba11))/2.d0
      dR(5,19,6)=-(((w1 - w4)*cab11 + (w1 - w4)*caba11 + b*sab11
     &- v1*sab11 + v4*sab11 + 2.d0*u1*sb1 - 2.d0*u4*sb1 + b*saba11
     &- v1*saba11 + v4*saba11)*sg1)/2.d0
      dR(5,20,4)=-(cb1*((b - v1 + v4)*ca1 + (-w1 + w4)*sa1)*sg1)
      dR(5,20,5)=(((b - v1 + v4)*cab11 - 2.d0*(u1 - u4)*cb1 - b*caba11
     &+ v1*caba11 - v4*caba11 - w1*sab11 + w4*sab11
     &+ w1*saba11 - w4*saba11)*sg1)/2.d0
      dR(5,20,6)=-(cg1*((w1 - w4)*ca1*cb1 + (u1 - u4)*sb1
     &+ 2.d0*b*ca1**2*s2b1*sg1 + cb1*(b*sa1 - v1*sa1 + v4*sa1
     &- 4.d0*b*sb1*sg1 + 4.d0*b*sa1**2*sb1*sg1)))
      dR(5,21,4)=((b - v1 + v4)*ca1 + (-w1 + w4)*sa1)*sb1
      dR(5,21,5)=(w1 - w4)*ca1*cb1 + (u1 - u4)*sb1
     &+ 2.d0*b*ca1**2*s2b1*sg1 + cb1*(b*sa1 - v1*sa1 + v4*sa1
     &- 4.d0*b*sb1*sg1 + 4.d0*b*sa1**2*sb1*sg1)
c
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr5 (sa1,sb1,sg1,ca1,cb1,cg1,saa2,sbb2,caa2,cbb2,
     &               so2,co2,b1,g1,a1,b2,g2,a2)
      do 908 i=1,3
         do 908 j=1,3
            dR(5,9+i,3+j) = dR1(i,j)
            dR(5,9+i,9+j) = dR2(i,j)
  908 continue
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr5 (sa1,sb1,sg1,ca1,cb1,cg1,saa3,sbb3,caa3,cbb3,
     &               so3,co3,b1,g1,a1,b3,g3,a3)
      do 909 i=1,3
         do 909 j=1,3
            dR(5,15+i,3+j) = dR1(i,j)
            dR(5,15+i,15+j) = dR2(i,j)
  909 continue
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr5 (sa1,sb1,sg1,ca1,cb1,cg1,saa4,sbb4,caa4,cbb4,
     &               so4,co4,b1,g1,a1,b4,g4,a4)
      do 910 i=1,3
         do 910 j=1,3
            dR(5,21+i,3+j) = dR1(i,j)
            dR(5,21+i,21+j) = dR2(i,j)
  910 continue
c
      do 911 i=1,3
         do 911 j=1,3
            dR(5,12+i,j) = dR(5,6+i,j)
            dR(5,18+i,j) = dR(5,6+i,j)
            dR(5,6+i,6+j) = -dR(5,6+i,j)
            dR(5,12+i,12+j) = -dR(5,6+i,j)
            dR(5,18+i,18+j) = -dR(5,6+i,j)
  911 continue
c
      dR(6,7,1)=cb1*sg1
      dR(6,7,2)=-(ca1*cg1) + sa1*sb1*sg1
      dR(6,7,3)=-(cg1*sa1) - ca1*sb1*sg1
      dR(6,8,1)=cb1*cg1
      dR(6,8,2)=cg1*sa1*sb1 + ca1*sg1
      dR(6,8,3)=-(ca1*cg1*sb1) + sa1*sg1
c
      dR(6,7,4)=ca1*((-w1 + w2)*cg1 + (v1 - v2)*sb1*sg1)
     &+ sa1*((v1 - v2)*cg1 + (w1 - w2)*sb1*sg1)
      dR(6,7,5)=(-((w1 - w2)*ca1*cb1) + (a - u1 + u2)*sb1 + cb1*(v1*sa1
     &- v2*sa1 - 4.d0*a*cg1*sb1 + 4.d0*a*cg1*sa1**2*sb1)
     &+ 2.d0*a*ca1**2*cg1*s2b1)*sg1
      dR(6,7,6)=-((a - u1 + u2)*cb1*cg1) - 2.d0*a*cb1**2*cg1**2*sa1**2
     &+ 2.d0*a*cb1**2*c2g1*sa1**2 - cg1*(w1*ca1 - w2*ca1 - v1*sa1
     &+ v2*sa1)*sb1 + sg1*(v1*ca1 - v2*ca1 + w1*sa1 - w2*sa1
     &+ 2.d0*a*sa1**2*sg1 - 2.d0*a*sa1**2*sb1**2*sg1)
      dR(6,8,4)=sa1*((w1 - w2)*cg1*sb1 + (-v1 + v2)*sg1)
     &+ ca1*((v1 - v2)*cg1*sb1 + (w1 - w2)*sg1)
      dR(6,8,5)=-(cg1*((w1 - w2)*cab11 + (w1 - w2)*caba11 - v1*sab11
     &+ v2*sab11 - 2.d0*a*sb1 + 2.d0*u1*sb1 - 2.d0*u2*sb1
     &- v1*saba11 + v2*saba11))/2.d0
      dR(6,8,6)=(a*cb1 - u1*cb1 + u2*cb1 - v1*sa1*sb1 + v2*sa1*sb1)*sg1
     &+ ca1*((v1 - v2)*cg1 + (w1 - w2)*sb1*sg1) + cg1*(w1*sa1 - w2*sa1
     &- 4.d0*a*cb1**2*sg1 + 4.d0*a*sa1**2*sg1
     &- 4.d0*a*sa1**2*sb1**2*sg1) + 2.d0*a*ca1**2*cb1**2*s2g1
c
      dR(6,13,4)=-(ca1*((w1 - w3)*cg1 + (b - v1 + v3)*sb1*sg1))
     &- sa1*((b - v1 + v3)*cg1 + (-w1 + w3)*sb1*sg1)
      dR(6,13,5)=-(((w1 - w3)*cab11 + (w1 - w3)*caba11 + b*sab11
     &- v1*sab11 + v3*sab11 - 2.d0*a*sb1 + 2.d0*u1*sb1 - 2.d0*u3*sb1
     &+ b*saba11 - v1*saba11 + v3*saba11)*sg1)/2.d0
      dR(6,13,6)=-((a - u1 + u3)*cb1*cg1) - w1*ca1*cg1*sb1
     &+ w3*ca1*cg1*sb1 - b*cg1*sa1*sb1 + v1*cg1*sa1*sb1
     &- v3*cg1*sa1*sb1 - b*ca1*sg1 + v1*ca1*sg1 - v3*ca1*sg1
     &+ w1*sa1*sg1 - w3*sa1*sg1 + 2.d0*a*sa1**2*sg1**2
     &- 2.d0*a*sa1**2*sb1**2*sg1**2 - 2.d0*cb1**2*sa1**2*sg1*(2.d0*b*cg1
     &+ a*sg1) + 2.d0*b*sa1**2*s2g1 - 2.d0*b*sa1**2*sb1**2*s2g1
      dR(6,14,4)=sa1*(w1*cg1*sb1 - w3*cg1*sb1 + b*sg1 - v1*sg1 + v3*sg1)
     &- ca1*((b - v1 + v3)*cg1*sb1 + (-w1 + w3)*sg1)
      dR(6,14,5)=-(cg1*((w1 - w3)*cab11 + (w1 - w3)*caba11 + b*sab11
     &- v1*sab11 + v3*sab11 - 2.d0*a*sb1 + 2.d0*u1*sb1 - 2.d0*u3*sb1
     &+ b*saba11 - v1*saba11 + v3*saba11))/2.d0
      dR(6,14,6)=(a - u1 + u3)*cb1*sg1 - ca1*((b - v1 + v3)*cg1
     &+ (-w1 + w3)*sb1*sg1) - 2.d0*cb1**2*(b*c2g1 + a*s2g1)
     &+ 2.d0*ca1**2*cb1**2*(b*c2g1 + a*s2g1) + sa1*(w1*cg1 - w3*cg1
     &+ 2.d0*b*cg1**2*sa1 - 2.d0*b*cg1**2*sa1*sb1**2 + b*sb1*sg1
     &- v1*sb1*sg1 + v3*sb1*sg1 - 2.d0*b*sa1*sg1**2 + 2.d0*b*sa1
     &*sb1**2*sg1**2 + 2.d0*a*sa1*s2g1 - 2.d0*a*sa1*sb1**2*s2g1)
c
      dR(6,19,4)=-(ca1*((w1 - w4)*cg1 + (b - v1 + v4)*sb1*sg1))
     &- sa1*((b - v1 + v4)*cg1 + (-w1 + w4)*sb1*sg1)
      dR(6,19,5)=-(((w1 - w4)*cab11 + (w1 - w4)*caba11 + b*sab11
     &- v1*sab11 + v4*sab11 + 2.d0*u1*sb1 - 2.d0*u4*sb1
     &+ b*saba11 - v1*saba11 + v4*saba11)*sg1)/2.d0
      dR(6,19,6)=(u1 - u4)*cb1*cg1 - ca1*(w1*cg1*sb1 - w4*cg1*sb1
     &+ b*sg1 - v1*sg1 + v4*sg1) - 2.d0*b*cb1**2*sa1**2*s2g1
     &- sa1*(b*cg1*sb1 - v1*cg1*sb1 + v4*cg1*sb1 - w1*sg1 + w4*sg1
     &- 2.d0*b*sa1*s2g1 + 2.d0*b*sa1*sb1**2*s2g1)
      dR(6,20,4)=sa1*(w1*cg1*sb1 - w4*cg1*sb1 + b*sg1 - v1*sg1 + v4*sg1)
     &- ca1*((b - v1 + v4)*cg1*sb1 + (-w1 + w4)*sg1)
      dR(6,20,5)=-(cg1*((w1 - w4)*ca1*cb1 + (u1 - u4)*sb1
     &+ 2.d0*b*ca1**2*s2b1*sg1 + cb1*(b*sa1 - v1*sa1 + v4*sa1
     &- 4.d0*b*sb1*sg1 + 4.d0*b*sa1**2*sb1*sg1)))
      dR(6,20,6)=-2.d0*b*cb1**2*c2g1 + 2.d0*b*ca1**2*cb1**2*c2g1
     &- (u1 - u4)*cb1*sg1 - ca1*((b - v1 + v4)*cg1 + (-w1 + w4)*sb1*sg1)
     &+ sa1*(w1*cg1 - w4*cg1 + 2.d0*b*cg1**2*sa1
     &- 2.d0*b*cg1**2*sa1*sb1**2 + b*sb1*sg1 - v1*sb1*sg1 + v4*sb1*sg1
     &- 2.d0*b*sa1*sg1**2 + 2.d0*b*sa1*sb1**2*sg1**2)
c
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr6 (sa1,sb1,sg1,ca1,cb1,cg1,saa2,sbb2,caa2,cbb2,
     &               so2,co2,b1,g1,a1,b2,g2,a2)
      do 912 i=1,3
         do 912 j=1,3
            dR(6,9+i,3+j) = dR1(i,j)
            dR(6,9+i,9+j) = dR2(i,j)
  912 continue
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr6 (sa1,sb1,sg1,ca1,cb1,cg1,saa3,sbb3,caa3,cbb3,
     &               so3,co3,b1,g1,a1,b3,g3,a3)
      do 913 i=1,3
         do 913 j=1,3
            dR(6,15+i,3+j) = dR1(i,j)
            dR(6,15+i,15+j) = dR2(i,j)
  913 continue
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr6 (sa1,sb1,sg1,ca1,cb1,cg1,saa4,sbb4,caa4,cbb4,
     &               so4,co4,b1,g1,a1,b4,g4,a4)
      do 914 i=1,3
         do 914 j=1,3
            dR(6,21+i,3+j) = dR1(i,j)
            dR(6,21+i,21+j) = dR2(i,j)
  914 continue
c
      do 915 i=1,3
         do 915 j=1,3
            dR(6,12+i,j) = dR(6,6+i,j)
            dR(6,18+i,j) = dR(6,6+i,j)
            dR(6,6+i,6+j) = -dR(6,6+i,j)
            dR(6,12+i,12+j) = -dR(6,6+i,j)
            dR(6,18+i,18+j) = -dR(6,6+i,j)
  915 continue
c
      do 916 i=1,3
         do 916 j=1,3
            dR(7,6+i,3+j) = -dR(1,i+6,3+j)
  916 continue
c
      do 917 i=1,3
         do 917 j=1,3
            dR(8,6+i,3+j) = -dR(2,i+6,3+j)
  917 continue
c
      do 918 i=1,3
         do 918 j=1,3
            dR(9,6+i,3+j) = -dR(3,i+6,3+j)
  918 continue
c
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr10 (sa1,sb1,sg1,ca1,cb1,cg1,saa2,sbb2,caa2,cbb2,so2,co2,
     &           cabn2,caba2,b1,g1,b2,g2,dR1,dR2)
      do 919 i=1,3
         do 919 j=1,3
            dR(10,9+i,3+j) = dR1(i,j)
            dR(10,9+i,9+j) = dR2(i,j)
  919 continue
c
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr11 (sa1,sb1,sg1,ca1,cb1,cg1,saa2,sbb2,caa2,cbb2,
     &           co2,b1,g1,b2,g2,dR1,dR2)
      do 920 i=1,3
         do 920 j=1,3
            dR(11,9+i,3+j) = dR1(i,j)
            dR(11,9+i,9+j) = dR2(i,j)
  920 continue
c
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr12 (sa1,sb1,sg1,ca1,cb1,cg1,saa2,sbb2,caa2,cbb2,
     &           co2,so2,b1,g1,b2,g2,dR1,dR2)
      do 921 i=1,3
         do 921 j=1,3
            dR(12,9+i,3+j) = dR1(i,j)
            dR(12,9+i,9+j) = dR2(i,j)
  921 continue
c
      do 922 i=1,3
         do 922 j=1,3
            dR(13,12+i,3+j) = -dR(1,i+6,3+j)
  922 continue
c
      do 923 i=1,3
         do 923 j=1,3
            dR(14,12+i,3+j) = -dR(2,i+6,3+j)
  923 continue
c
      do 924 i=1,3
         do 924 j=1,3
            dR(15,12+i,3+j) = -dR(3,i+6,3+j)
  924 continue
c
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr10 (sa1,sb1,sg1,ca1,cb1,cg1,saa3,sbb3,caa3,cbb3,so3,co3,
     &           cabn3,caba3,b1,g1,b3,g3,dR1,dR2)
      do 925 i=1,3
         do 925 j=1,3
            dR(16,15+i,3+j) = dR1(i,j)
            dR(16,15+i,15+j) = dR2(i,j)
  925 continue
c
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr11 (sa1,sb1,sg1,ca1,cb1,cg1,saa3,sbb3,caa3,cbb3,
     &           co3,b1,g1,b3,g3,dR1,dR2)
      do 926 i=1,3
         do 926 j=1,3
            dR(17,15+i,3+j) = dR1(i,j)
            dR(17,15+i,15+j) = dR2(i,j)
  926 continue
c
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr12 (sa1,sb1,sg1,ca1,cb1,cg1,saa3,sbb3,caa3,cbb3,
     &           co3,so3,b1,g1,b3,g3,dR1,dR2)
      do 927 i=1,3
         do 927 j=1,3
            dR(18,15+i,3+j) = dR1(i,j)
            dR(18,15+i,15+j) = dR2(i,j)
  927 continue
c
      do 928 i=1,3
         do 928 j=1,3
            dR(19,18+i,3+j) = -dR(1,i+6,3+j)
  928 continue
c
      do 929 i=1,3
         do 929 j=1,3
            dR(20,18+i,3+j) = -dR(2,i+6,3+j)
  929 continue
c
      do 930 i=1,3
         do 930 j=1,3
            dR(21,18+i,3+j) = -dR(3,i+6,3+j)
  930 continue
c
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr10 (sa1,sb1,sg1,ca1,cb1,cg1,saa4,sbb4,caa4,cbb4,so4,co4,
     &           cabn4,caba4,b1,g1,b4,g4,dR1,dR2)
      do 931 i=1,3
         do 931 j=1,3
            dR(22,21+i,3+j) = dR1(i,j)
            dR(22,21+i,21+j) = dR2(i,j)
  931 continue
c
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr11 (sa1,sb1,sg1,ca1,cb1,cg1,saa3,sbb3,caa3,cbb3,
     &           co3,b1,g1,b3,g3,dR1,dR2)
      do 932 i=1,3
         do 932 j=1,3
            dR(23,21+i,3+j) = dR1(i,j)
            dR(23,21+i,21+j) = dR2(i,j)
  932 continue
c
      call clrmx (3,dR1)
      call clrmx (3,dR2)
      call dr12 (sa1,sb1,sg1,ca1,cb1,cg1,saa4,sbb4,caa4,cbb4,
     &           co4,so4,b1,g1,b4,g4,dR1,dR2)
      do 933 i=1,3
         do 933 j=1,3
            dR(24,21+i,3+j) = dR1(i,j)
            dR(24,21+i,21+j) = dR2(i,j)
  933 continue
c
      call clrmx (24,tl)
      call clrmx (24,te)
      call trns(tl,d(4),d(5),d(6))
      call euler (te,d(4),d(5),d(10),d(11),d(16),d(17),d(22),d(23))
c
c********************************************************************
c  dE(transpose) R(transpose) K T(transpose) r の計算
c
      call clrv (24,dmy1)
      call clrv (24,dmy2)
      call clrmx (24,sk1)
      call mxtv (24,tl,dl,dmy1)
      call mxv (24,s,dmy1,dmy2)
      call mxtv (24,R,dmy2,dmy1)
      do 30 k=1,24
         do 30 i=1,24
            do 30 j=1,24
               sk1(i,k) = sk1(i,k) + dte(k,j,i) * dmy1(j)
   30 continue
c
c********************************************************************
c  E(transpose) dR(transpose) k T(transpose) r の計算
c
      call clrv (24,dmy1)
      call clrv (24,dmy2)
      call clrmx (24,dmy3)
      call clrmx (24,dmy4)
      call mxtv (24,tl,dl,dmy1)
      call mxv (24,s,dmy1,dmy2)
      do 40 k=1,24
         do 40 i=1,24
            do 40 j=1,24
               dmy3(i,k) = dmy3(i,k) + dR(k,j,i) * dmy2(j)
   40 continue
c
      call mxtmx (24,te,dmy3,dmy4)
      do 700 i=1,24
         do 700 j=1,24
            sk1(i,j) = sk1(i,j) + dmy4(i,j)
  700 continue
c
c********************************************************************
c  E(transpose) R(transpose) k R の計算
c
      call clrmx (24,dmy3)
      call clrmx (24,dmy4)
      call mxmx (24,s,r,dmy3)
      call mxtmx (24,r,dmy3,dmy4)
      call mxtmx(24,te,dmy4,dmy3)
      do 710 i=1,24
         do 710 j=1,24
            sk1(i,j) = sk1(i,j) + dmy3(i,j)
  710 continue
c
c********************************************************************
c
      call clrmx (24,sk)
      call mxmx (24,ti,sk1,sk)
      call mxmx (24,sk,te,sk1)
      call mxmxt (24,sk1,ti,sk)
c
c********************************************************************
c     不釣合い力の計算
c
      call clrv (24,dmy1)
      call clrv (24,dmy2)
      call mxtv (24,tl,dl,dmy1)
      call mxv (24,s,dmy1,dmy2)
      call mxtv (24,R,dmy2,dmy1)
      call mxtv (24,te,dmy1,dmy2)
      call mxv (24,ti,dmy2,dmy1)
      do 120 i=1,24
         ba(i) = ba(i) - dmy1(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 check!!
      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(7,1) = p4
      s(7,7) = p1
      s(7,2) = -p5
      s(13,1) = -p1 / 2.d0
      s(13,7) = p7
      s(13,13) = p1
      s(13,2) = -p2
      s(13,8) = -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(2,1) = p2
      s(2,2) = p3
      s(8,1) = p5
      s(8,7) = -p2
      s(8,2) = p6
      s(8,8) = p3
      s(14,1) = -p2
      s(14,7) = p5
      s(14,13) = p2
      s(14,2) = -p3 / 2.d0
      s(14,8) = p8
      s(14,14) = p3
      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(21,15) = s(9,3)
      s(21,9) = s(15,3)
      s(15,9) = 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(22,15) = - s(10,3)
      s(22,9) = s(16,3)
      s(16,9) = s(22,3)
      s(9,4) = s(10,3)
      s(21,16) = - s(10,3)
      s(15,4) = - s(16,3)
      s(21,10) = - s(16,3)
      s(15,10) = - s(22,3)
      s(21,4) = - s(22,3)
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(23,15) = - s(11,3)
      s(23,9) = - s(17,3)
      s(17,9) = - s(23,3)
      s(9,5) = - s(11,3)
      s(21,17) = s(11,3)
      s(15,5) = - s(17,3)
      s(21,11) = s(17,3)
      s(15,11) = - s(23,3)
      s(21,5) = s(23,3)
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(22,16) = s(10,4)
      s(22,10) = s(16,4)
      s(16,10) = s(22,4)
c
      s(5,4) = - 30.d0 * um * a * b
      s(5,4) = s(5,4) * p15
      s(11,10) = - s(5,4)
      s(17,16) = s(5,4)
      s(23,22) = - s(5,4)
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(23,17) = s(11,5)
      s(23,11) = s(17,5)
      s(17,11) = s(23,5)
c
      p9 = e * th * a * b * 0.03d0
         s(6, 6) = p9
         s(12, 6) = -p9 / 3.d0
         s(12, 12) = p9
         s(18, 6) = -p9 / 3.d0
         s(18, 12) = -p9 / 3.d0
         s(18, 18) = p9
         s(24, 6) = -p9 / 3.d0
         s(24, 12) = -p9 / 3.d0
         s(24, 18) = -p9 / 3.d0
         s(24, 24) = p9
c
      do 100 i=1,24
         do 100 j=1,24
            if( i.gt.j ) then
               s(j,i) = s(i,j)
            end if
  100 continue
      return
      end
c====================================================================
      subroutine trns (tl,al,bt,gm)
      implicit real*8(a-h,o-z)
      dimension tl(24,24)
c
      call clrmx (24,tl)
c check!!
      tl(1,1) = cos(gm)*cos(bt)
      tl(1,2) = -sin(gm)*cos(bt)
      tl(1,3) = sin(bt)
      tl(2,1) = sin(gm)*cos(al)+cos(gm)*sin(bt)*sin(al)
      tl(2,2) = cos(gm)*cos(al)-sin(gm)*sin(bt)*sin(al)
      tl(2,3) = -cos(bt)*sin(al)
      tl(3,1) = sin(gm)*sin(al)-cos(gm)*sin(bt)*cos(al)
      tl(3,2) = cos(gm)*sin(al)+sin(gm)*sin(bt)*cos(al)
      tl(3,3) = cos(bt)*cos(al)
c
      do 10 i=1,3
         do 10 j=1,3
            tl(i+3,j+3) = tl(i,j)
            tl(i+6,j+6) = tl(i,j)
            tl(i+9,j+9) = tl(i,j)
            tl(i+12,j+12) = tl(i,j)
            tl(i+15,j+15) = tl(i,j)
            tl(i+18,j+18) = tl(i,j)
            tl(i+21,j+21) = tl(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
      call clrmx (24,te)
      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 dseuler(ste,a,b,n)
      implicit real*8 (a-h,o-z)
      dimension ste(3,3)
c
      call clrmx(3,ste)
      if( n.eq.1 ) then
         ste(1,2) = cos(a)*tan(b)
         ste(1,3) = sin(a)*tan(b)
         ste(2,2) = -sin(a)
         ste(2,3) = cos(a)
         ste(3,2) = -cos(a)/cos(b)
         ste(3,3) = -sin(a)/cos(b)
      else if( n.eq.2 ) then
         ste(1,2) = sin(a)/cos(b)/cos(b)
         ste(1,3) = -cos(a)/cos(b)/cos(b)
         ste(3,2) = -sin(a)*tan(b)/cos(b)
         ste(3,3) = cos(a)*tan(b)/cos(b)
      end if
c
      return
      end
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
c
c====================================================================
      subroutine dr10 (sa1,sb1,sg1,ca1,cb1,cg1,saa,sbb,caa,cbb,so,co,
     &                cabn,caba,b1,g1,b,g,dR1,dR2)
      implicit real*8 (a-h,o-z)
c
      dimension  dR1(3,3), dR2(3,3)
      call clrmx (3,dR1)
      call clrmx (3,dR2)
c
      dR1(1,1)=(-(((b1 - b)*caa + (-g1 + g)*cbb*saa)
     &*(cg1*sa1*sb1 + ca1*sg1)) + ((g1 - g)*caa*cbb
     &+ (b1 - b)*saa)*(ca1*cg1*sb1 - sa1*sg1))/4.d0
      dR1(1,2)=(2.d0*cb1*cg1*sa1*((g1 - g)*caa*cbb
     &+ (b1 - b)*saa) + 2.d0*ca1*cb1*cg1*((b1 - b)*caa
     &+ (-g1 + g)*cbb*saa) - 4.d0*cg1*sb1 + (2.d0*saa + (-g1 + g)
     &*caa*sbb)*(cg1*sa1*sb1 + ca1*sg1) + (2.d0*caa
     &+ (g1 - g)*saa*sbb)*(ca1*cg1*sb1 - sa1*sg1))/4.d0
      dR1(1,3)=(-2.d0*cb1*sg1 + caa*cbb*(cg1*sa1*sb1 + ca1*sg1)
     &+ cbb*saa*(-(ca1*cg1*sb1) + sa1*sg1) - ((b1 - b)*caa
     &+ (-g1 + g)*cbb*saa)*(cg1*sa1 + ca1*sb1*sg1)
     &+ ((g1 - g)*caa*cbb + (b1 - b)*saa)*(ca1*cg1
     &- sa1*sb1*sg1))/2.d0
      dR1(2,1)=-(((g1 - g)*caa*cbb + (b1 - b)*saa)*(cg1*sa1
     &+ ca1*sb1*sg1))/4.d0 - (((b1 - b)*caa + (-g1 + g)*cbb*saa)
     &*(ca1*cg1 - sa1*sb1*sg1))/4.d0
      dR1(2,2)=(-2.d0*cb1*sa1*((g1 - g)*caa*cbb + (b1 - b)*saa)
     &*sg1 - 2.d0*ca1*cb1*((b1 - b)*caa + (-g1 + g)*cbb*saa)*sg1
     &+ 4.d0*sb1*sg1 - (2.d0*caa + (g1 - g)*saa*sbb)*(cg1*sa1
     &+ ca1*sb1*sg1) + (2.d0*saa + (-g1 + g)*caa*sbb)
     &*(ca1*cg1 - sa1*sb1*sg1))/4.d0
      dR1(2,3)=(-2.d0*cb1*cg1 - ((g1 - g)*caa*cbb + (b1 - b)
     &*saa)*(cg1*sa1*sb1 + ca1*sg1) - ((b1 - b)*caa
     &+ (-g1 + g)*cbb*saa)*(ca1*cg1*sb1 - sa1*sg1)
     &+ cbb*saa*(cg1*sa1 + ca1*sb1*sg1)
     &+ caa*cbb*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      dR1(3,1)=-(cb1*((g1 - g)*cabn
     &+ (g1 - g)*caba
     &+ 2.d0*(-b1 + b)*so))/8.d0
      dR1(3,2)=(2.d0*((b1 - b)*co + (g1 - g)*cbb
     &*so)*sb1 + cb1*(4.d0 - 2.d0*sa1*saa
     &+ g1*caa*sa1*sbb - g*caa*sa1*sbb + ca1*(-2.d0*caa
     &- g1*saa*sbb + g*saa*sbb)))/4.d0
      dR1(3,3)=-(cb1*cbb*so)/2.d0
c
      dR2(1,1)=(((b1 - b)*caa + (-g1 + g)*cbb*saa)
     &*(cg1*sa1*sb1 + ca1*sg1) - ((g1 - g)*caa*cbb + (b1 - b)*saa)
     &*(ca1*cg1*sb1 - sa1*sg1))/4.d0
      dR2(1,2)=-((2.d0*saa + (g1 - g)*caa*sbb)*(cg1*sa1*sb1
     &+ ca1*sg1))/4.d0 - ((2.d0*caa + (-g1 + g)*saa*sbb)
     &*(ca1*cg1*sb1 - sa1*sg1))/4.d0
      dR2(1,3)=-(cbb*(caa*(cg1*sa1*sb1 + ca1*sg1)
     &+ saa*(-(ca1*cg1*sb1) + sa1*sg1)))/2.d0
      dR2(2,1)=(((g1 - g)*caa*cbb + (b1 - b)*saa)*(cg1*sa1
     &+ ca1*sb1*sg1) + ((b1 - b)*caa + (-g1 + g)*cbb*saa)*(ca1*cg1
     &- sa1*sb1*sg1))/4.d0
      dR2(2,2)=((2.d0*caa + (-g1 + g)*saa*sbb)*(cg1*sa1
     &+ ca1*sb1*sg1) - (2.d0*saa + (g1 - g)*caa*sbb)
     &*(ca1*cg1 - sa1*sb1*sg1))/4.d0
      dR2(2,3)=(cbb*(-(saa*(cg1*sa1 + ca1*sb1*sg1))
     &- caa*(ca1*cg1 - sa1*sb1*sg1)))/2.d0
      dR2(3,1)=(cb1*(ca1*((g1 - g)*caa*cbb + (b1 - b)*saa)
     &+ sa1*((-b1 + b)*caa + (g1 - g)*cbb*saa)))/4.d0
      dR2(3,2)=(cb1*(sa1*(2.d0*saa + (g1 - g)*caa*sbb)
     &+ ca1*(2.d0*caa + (-g1 + g)*saa*sbb)))/4.d0
      dR2(3,3)=(cb1*cbb*so)/2.d0
c
      return
      end
c
c====================================================================
      subroutine dr11 (sa1,sb1,sg1,ca1,cb1,cg1,saa,sbb,caa,cbb,
     &           co,b1,g1,b,g,dR1,dR2)
      implicit real*8 (a-h,o-z)
c
      dimension  dR1(3,3), dR2(3,3)
      call clrmx (3,dR1)
      call clrmx (3,dR2)
c
      dR1(1,1)=(sa1*(cg1*sb1*(2.d0*saa + (g1 - g)*caa*sbb)
     &+ (-2.d0*caa + g1*saa*sbb - g*saa*sbb)*sg1) + ca1*(saa
     &*(-(g1*cg1*sb1*sbb) + g*cg1*sb1*sbb + 2.d0*sg1)
     &+ caa*(2.d0*cg1*sb1 + g1*sbb*sg1 - g*sbb*sg1)))/4.d0
      dR1(1,2)=(-((g1 - g)*cg1*(-2.d0*cbb*sb1 + cbb*sa1*saa*sb1
     &- cb1*sbb + 2.d0*cb1*sa1*saa*sbb)) + ca1*saa*(-4.d0*cb1*cg1
     &- g1*cbb*sg1 + g*cbb*sg1) + caa*(2.d0*cb1*cg1*(2.d0*sa1
     &- g1*ca1*sbb + g*ca1*sbb)
     &- (g1 - g)*cbb*(ca1*cg1*sb1 - sa1*sg1)))/4.d0
      dR1(1,3)=(-(cb1*cbb*cg1) + (g1 - g)*cb1*cbb*sg1 - saa*sbb
     &*(cg1*sa1*sb1 + ca1*sg1) + caa*sbb*(-(ca1*cg1*sb1) + sa1*sg1)
     &+ 2.d0*(saa + ((g1 - g)*caa*sbb)/2.d0)*(cg1*sa1 + ca1*sb1*sg1)
     &+ 2.d0*(caa + ((-g1 + g)*saa*sbb)/2.d0)*(ca1*cg1
     &- sa1*sb1*sg1))/2.d0
      dR1(2,1)=-((2.d0*caa + (-g1 + g)*saa*sbb)*(cg1*sa1
     &+ ca1*sb1*sg1))/4.d0 + (saa + ((g1 - g)*caa*sbb)/2.d0)*(ca1*cg1
     &- sa1*sb1*sg1) - ((2.d0*saa + (g1 - g)*caa*sbb)
     &*(ca1*cg1 - sa1*sb1*sg1))/4.d0
      dR1(2,2)=((g1 - g)*(-2.d0*cbb*sb1 + cbb*sa1*saa*sb1 - cb1*sbb
     &+ 2.d0*cb1*sa1*saa*sbb)*sg1 + ca1*saa*(-(g1*cbb*cg1) + g*cbb*cg1
     &+ 4.d0*cb1*sg1) + caa*(-2.d0*cb1*(2.d0*sa1 - g1*ca1*sbb
     &+ g*ca1*sbb)*sg1 + (g1 - g)*cbb*(cg1*sa1 + ca1*sb1*sg1)))/4.d0
      dR1(2,3)=((g1 - g)*cb1*cbb*cg1 + cb1*cbb*sg1 + 2.d0*(caa
     &+ ((-g1 + g)*saa*sbb)/2.d0)*(-(cg1*sa1*sb1) - ca1*sg1)
     &+ 2.d0*(saa + ((g1 - g)*caa*sbb)/2.d0)*(ca1*cg1*sb1 - sa1*sg1)
     &+ caa*sbb*(cg1*sa1 + ca1*sb1*sg1)
     &- saa*sbb*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      dR1(3,1)=-(cb1*(sa1*(2.d0*saa + (g1 - g)*caa*sbb)
     &+ ca1*(2.d0*caa + (-g1 + g)*saa*sbb)))/4.d0
      dR1(3,2)=((g1 - g)*(-2.d0 + co)*cb1*cbb + sb1*(-4.d0*ca1*saa
     &+ g1*sbb - g*sbb - 2.d0*g1*sa1*saa*sbb + 2.d0*g*sa1*saa*sbb
     &+ 2.d0*caa*(2.d0*sa1 - g1*ca1*sbb + g*ca1*sbb)))/4.d0
      dR1(3,3)=(-(cbb*sb1) + co*cb1*sbb)/2.d0
      dR2(1,1)=-((2*saa + (g1 - g)*caa*sbb)*(cg1*sa1*sb1 + ca1*sg1))
     &/4.d0 - ((2.d0*caa + (-g1 + g)*saa*sbb)
     &*(ca1*cg1*sb1 - sa1*sg1))/4.d0
      dR2(1,2)=((-g1 + g)*(-(cb1*cg1*sbb) + cbb*saa*(cg1*sa1*sb1
     &+ ca1*sg1) + caa*cbb*(ca1*cg1*sb1 - sa1*sg1)))/4.d0
      dR2(1,3)=(cb1*cbb*cg1 + saa*sbb*(cg1*sa1*sb1 + ca1*sg1)
     &+ caa*sbb*(ca1*cg1*sb1 - sa1*sg1))/2.d0
      dR2(2,1)=((2.d0*caa + (-g1 + g)*saa*sbb)*(cg1*sa1 + ca1*sb1*sg1)
     &- (2.d0*saa + (g1 - g)*caa*sbb)*(ca1*cg1 - sa1*sb1*sg1))/4.d0
      dR2(2,2)=((-g1 + g)*(cb1*sbb*sg1 - caa*cbb*(cg1*sa1 + ca1*sb1*sg1)
     &+ cbb*saa*(ca1*cg1 - sa1*sb1*sg1)))/4.d0
      dR2(2,3)=(-(cb1*cbb*sg1) - caa*sbb*(cg1*sa1 + ca1*sb1*sg1)
     &+ saa*sbb*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      dR2(3,1)=(cb1*(sa1*(2.d0*saa + (g1 - g)*caa*sbb)
     &+ ca1*(2.d0*caa + (-g1 + g)*saa*sbb)))/4.d0
      dR2(3,2)=((g1 - g)*(ca1*caa*cb1*cbb + cb1*cbb*sa1*saa
     &+ sb1*sbb))/4.d0
      dR2(3,3)=(cbb*sb1 - co*cb1*sbb)/2.d0
c
      return
      end
c====================================================================
      subroutine dr12 (sa1,sb1,sg1,ca1,cb1,cg1,saa,sbb,caa,cbb,
     &           co,so,b1,g1,b,g,dR1,dR2)
      implicit real*8 (a-h,o-z)
c
      dimension  dR1(3,3), dR2(3,3)
      call clrmx (3,dR1)
      call clrmx (3,dR2)
c
      dR1(1,1)=(cbb*(caa*(cg1*sa1*sb1 + ca1*sg1)
     &+ saa*(-(ca1*cg1*sb1) + sa1*sg1)))/2.d0
      dR1(1,2)=(-((-1.d0 + 2.d0*co)*cb1*cbb*cg1)
     &+ sbb*((-2.d0 + co)*cg1*sb1 - so*sg1))/2.d0
      dR1(1,3)=-(cb1*sbb*sg1) + caa*cbb*(cg1*sa1 + ca1*sb1*sg1)
     &- cbb*saa*(ca1*cg1 - sa1*sb1*sg1)
      dR1(2,1)=(cbb*(sa1*(cg1*saa - caa*sb1*sg1)
     &+ ca1*(caa*cg1 + saa*sb1*sg1)))/2.d0
      dR1(2,2)=-(cb1*cbb*sg1)/2.d0 + ca1*caa*cb1*cbb*sg1
     &+ cb1*cbb*sa1*saa*sg1 + sb1*sbb*sg1 - (caa*sbb*(cg1*sa1
     &+ ca1*sb1*sg1))/2.d0 + (saa*sbb*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      dR1(2,3)=-(cb1*cg1*sbb) + cbb*saa*(cg1*sa1*sb1 + ca1*sg1)
     &+ caa*cbb*(ca1*cg1*sb1 - sa1*sg1)
      dR1(3,1)=-(cb1*cbb*so)/2.d0
      dR1(3,2)=((1.d0 - 2.d0*co)*cbb*sb1 - (-2.d0 + co)*cb1*sbb)/2.d0
      dR2(1,1)=-(cbb*(caa*(cg1*sa1*sb1 + ca1*sg1)
     &+ saa*(-(ca1*cg1*sb1) + sa1*sg1)))/2.d0
      dR2(1,2)=(cb1*cbb*cg1 + saa*sbb*(cg1*sa1*sb1 + ca1*sg1)
     &+ caa*sbb*(ca1*cg1*sb1 - sa1*sg1))/2.d0
      dR2(2,1)=(cbb*(-(saa*(cg1*sa1 + ca1*sb1*sg1)) - caa*(ca1*cg1
     &- sa1*sb1*sg1)))/2.d0
      dR2(2,2)=(-(cb1*cbb*sg1) - caa*sbb*(cg1*sa1 + ca1*sb1*sg1)
     &+ saa*sbb*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      dR2(3,1)=(cb1*cbb*so)/2.d0
      dR2(3,2)=(cbb*sb1 - co*cb1*sbb)/2.d0
c
      return
      end
c====================================================================
c
      subroutine dr4(sa1,sb1,sg1,ca1,cb1,cg1,saa,sbb,caa,cbb,
     &               so,co,cabn,caba,b1,g1,b,g)
      implicit real*8 (a-h,o-z)
c
      dimension dR1(3,3),dR2(3,3)
      call clrmx (3,dR1)
      call clrmx (3,dR2)
c
      dR1(1,1)=(sa1*(saa*(-(g1*cbb*cg1*sb1) + g*cbb*cg1*sb1 + b1*sg1
     &- b*sg1) + caa*(b1*cg1*sb1 - b*cg1*sb1 + g1*cbb*sg1 - g*cbb*sg1))
     &+ ca1*(caa*(-(g1*cbb*cg1*sb1) + g*cbb*cg1*sb1 + b1*sg1 - b*sg1)
     &+ saa*(-(b1*cg1*sb1) + b*cg1*sb1 - g1*cbb*sg1 + g*cbb*sg1)))/4.d0
      dR1(1,2)=(-2.d0*b1*cb1*cg1*sa1*saa + 2.d0*b*cb1*cg1*sa1*saa
     &+ 4.d0*cg1*sb1 - 2.d0*cg1*sa1*saa*sb1 + g1*sa1*saa*sbb*sg1
     &- g*sa1*saa*sbb*sg1 + caa*sa1*(-2.d0*g1*cb1*cbb*cg1
     &+ 2.d0*g*cb1*cbb*cg1 + g1*cg1*sb1*sbb - g*cg1*sb1*sbb + 2.d0*sg1)
     &- ca1*(saa*(-2.d0*g1*cb1*cbb*cg1 + 2.d0*g*cb1*cbb*cg1
     &+ g1*cg1*sb1*sbb - g*cg1*sb1*sbb + 2.d0*sg1) 
     &+ caa*(2.d0*b1*cb1*cg1 - 2.d0*b*cb1*cg1 + 2.d0*cg1*sb1
     &- g1*sbb*sg1 + g*sbb*sg1)))/4.d0
      dR1(1,3)=(-(g1*cbb*cg1*sa1*saa) + g*cbb*cg1*sa1*saa + 2.d0*cb1*sg1
     &- cbb*sa1*saa*sg1 + b1*sa1*saa*sb1*sg1 - b*sa1*saa*sb1*sg1
     &+ caa*sa1*(b1*cg1 - b*cg1 - cbb*cg1*sb1 + g1*cbb*sb1*sg1
     &- g*cbb*sb1*sg1) + ca1*(-(caa*(g1*cbb*cg1 - g*cbb*cg1 + cbb*sg1
     &- b1*sb1*sg1 + b*sb1*sg1)) + saa*(-(b1*cg1) + b*cg1
     &+ cbb*cg1*sb1 - g1*cbb*sb1*sg1 + g*cbb*sb1*sg1)))/2.d0
      dR1(2,1)=(ca1*(saa*(-(g1*cbb*cg1) + g*cbb*cg1 + b1*sb1*sg1
     &- b*sb1*sg1) + caa*(b1*cg1 - b*cg1 + g1*cbb*sb1*sg1
     &- g*cbb*sb1*sg1)) + sa1*(caa*(g1*cbb*cg1 - g*cbb*cg1 - b1*sb1*sg1
     &+ b*sb1*sg1) + saa*(b1*cg1 - b*cg1 + g1*cbb*sb1*sg1
     &- g*cbb*sb1*sg1)))/4.d0
      dR1(2,2)=(-2.d0*cb1*sa1*((g1 - g)*caa*cbb + (b1 - b)*saa)*sg1
     &- 4.d0*cb1*sa1*((-g1 + g)*caa*cbb + (-b1 + b)*saa)*sg1
     &- 4.d0*ca1*cb1*((-b1 + b)*caa + (g1 - g)*cbb*saa)*sg1
     &- 2.d0*ca1*cb1*((b1 - b)*caa + (-g1 + g)*cbb*saa)*sg1
     &- 4.d0*sb1*sg1 + (2.d0*caa + (g1 - g)*saa*sbb)*(cg1*sa1
     &+ ca1*sb1*sg1) - (2.d0*saa + (-g1 + g)*caa*sbb)
     &*(ca1*cg1 - sa1*sb1*sg1))/4.d0
      dR1(2,3)=(2.d0*cb1*cg1 + 2.d0*((-g1 + g)*caa*cbb + (-b1 + b)*saa)
     &*(-(cg1*sa1*sb1) - ca1*sg1) - ((g1 - g)*caa*cbb + (b1 - b)*saa)
     &*(cg1*sa1*sb1 + ca1*sg1) - ((b1 - b)*caa + (-g1 + g)*cbb*saa)
     &*(ca1*cg1*sb1 - sa1*sg1) + 2.d0*((-b1 + b)*caa + (g1 - g)*cbb*saa)
     &*(-(ca1*cg1*sb1) + sa1*sg1) - cbb*saa*(cg1*sa1
     &+ ca1*sb1*sg1) - caa*cbb*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      dR1(3,1)=(cb1*((g1 - g)*cabn + (g1 - g)*caba
     &+ 2.d0*(-b1 + b)*so))/8.d0
      dR1(3,2)=(2.d0*((-b1 + b)*co + (-g1 + g)*cbb*so)*sb1
     &+ cb1*(-4.d0 + 2.d0*sa1*saa - g1*caa*sa1*sbb + g*caa*sa1*sbb
     &+ ca1*(2.d0*caa + (g1 - g)*saa*sbb)))/4.d0
      dR1(3,3)=(cb1*cbb*so)/2.d0
      dR2(1,1)=(-(((b1 - b)*caa + (-g1 + g)*cbb*saa)*(cg1*sa1*sb1
     &+ ca1*sg1)) + ((g1 - g)*caa*cbb + (b1 - b)*saa)
     &*(ca1*cg1*sb1 - sa1*sg1))/4.d0
      dR2(1,2)=(sa1*(cg1*sb1*(2.d0*saa + (g1 - g)*caa*sbb) + (-2.d0*caa
     &+ g1*saa*sbb - g*saa*sbb)*sg1) + ca1*(saa*(-(g1*cg1*sb1*sbb)
     &+ g*cg1*sb1*sbb + 2.d0*sg1) + caa*(2.d0*cg1*sb1
     &+ g1*sbb*sg1 - g*sbb*sg1)))/4.d0
      dR2(1,3)=(cbb*(caa*(cg1*sa1*sb1 + ca1*sg1)
     &+ saa*(-(ca1*cg1*sb1) + sa1*sg1)))/2.d0
      dR2(2,1)=-(((g1 - g)*caa*cbb + (b1 - b)*saa)*(cg1*sa1
     &+ ca1*sb1*sg1))/4.d0 - (((b1 - b)*caa + (-g1 + g)*cbb*saa)
     &*(ca1*cg1 - sa1*sb1*sg1))/4.d0
      dR2(2,2)=-((2.d0*caa + (-g1 + g)*saa*sbb)*(cg1*sa1
     &+ ca1*sb1*sg1))/4.d0 + (saa + ((g1 - g)*caa*sbb)/2.d0)*(ca1*cg1
     &- sa1*sb1*sg1) - ((2.d0*saa + (g1 - g)*caa*sbb)
     &*(ca1*cg1 - sa1*sb1*sg1))/4.d0
      dR2(2,3)=(cbb*(sa1*(cg1*saa - caa*sb1*sg1)
     &+ ca1*(caa*cg1 + saa*sb1*sg1)))/2
      dR2(3,1)=-(cb1*((g1 - g)*cabn + (g1 - g)*caba
     &+ 2.d0*(-b1 + b)*so))/8.d0
      dR2(3,2)=-(cb1*(sa1*(2.d0*saa + (g1 - g)*caa*sbb)
     &+ ca1*(2.d0*caa + (-g1 + g)*saa*sbb)))/4.d0
      dR2(3,3)=-(cb1*cbb*so)/2.d0
c
      return
      end
c====================================================================
c
      subroutine dr5(sa1,sb1,sg1,ca1,cb1,cg1,saa,sbb,caa,cbb,
     &               so,co,b1,g1,a1,b,g,a)
      implicit real*8 (a-h,o-z)
c
      dimension dR1(3,3),dR2(3,3)
      call clrmx (3,dR1)
      call clrmx (3,dR2)
c
      dR1(1,1)=(-2.d0*b1*cb1*cg1*sa1*saa + 2.d0*b*cb1*cg1*sa1*saa
     &+ 4.d0*cg1*sb1 - 2.d0*cg1*sa1*saa*sb1 + g1*sa1*saa*sbb*sg1
     &- g*sa1*saa*sbb*sg1 + caa*sa1*(-2.d0*g1*cb1*cbb*cg1
     &+ 2.d0*g*cb1*cbb*cg1 + g1*cg1*sb1*sbb - g*cg1*sb1*sbb
     &+ 2.d0*sg1) - ca1*(saa*(-2.d0*g1*cb1*cbb*cg1
     &+ 2.d0*g*cb1*cbb*cg1 + g1*cg1*sb1*sbb - g*cg1*sb1*sbb
     &+ 2.d0*sg1) + caa*(2.d0*b1*cb1*cg1 - 2.d0*b*cb1*cg1
     &+ 2.d0*cg1*sb1 - g1*sbb*sg1 + g*sbb*sg1)))/4.d0
      dR1(1,2)=(4.d0*(b1 - b)*cg1*so*sb1 + cb1*cg1*(4.d0*a1 - 4.d0*a
     &+ 8.d0*ca1*saa + 5.d0*g1*sbb - 5.d0*g*sbb - 4.d0*g1*sa1*saa*sbb
     &+ 4.d0*g*sa1*saa*sbb - 4.d0*caa*(2.d0*sa1 + g1*ca1*sbb
     &- g*ca1*sbb)) - (g1 - g)*cbb*(-4.d0*cg1*sb1
     &+ 5.d0*ca1*caa*cg1*sb1 + 5.d0*cg1*sa1*saa*sb1
     &- caa*sa1*sg1 + ca1*saa*sg1))/4.d0
      dR1(1,3)=(-(cb1*cbb*cg1) + 2.d0*ca1*caa*cb1*cbb*cg1
     &+ 2.d0*cb1*cbb*cg1*sa1*saa + 2.d0*cg1*sb1*sbb + (g1 - g)
     &*cb1*cbb*sg1 + 2.d0*ca1*cb1*((-g1 + g)*caa*cbb + (-b1 + b)*saa)
     &*sg1 - 2.d0*cb1*sa1*((-b1 + b)*caa + (g1 - g)*cbb*saa)*sg1
     &+ 2.d0*sb1*(-a1 + a + (-g1 + g)*sbb)*sg1 - saa*sbb*(cg1*sa1*sb1
     &+ ca1*sg1) + caa*sbb*(-(ca1*cg1*sb1) + sa1*sg1) - (2.d0*saa
     &+ (-g1 + g)*caa*sbb)*(cg1*sa1 + ca1*sb1*sg1) + 2.d0*(-caa
     &+ ((-g1 + g)*saa*sbb)/2.d0)*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      dR1(2,1)=(-2.d0*cb1*sa1*((g1 - g)*caa*cbb + (b1 - b)*saa)*sg1
     &- 4.d0*cb1*sa1*((-g1 + g)*caa*cbb + (-b1 + b)*saa)*sg1
     &- 4.d0*ca1*cb1*((-b1 + b)*caa + (g1 - g)*cbb*saa)*sg1
     &- 2.d0*ca1*cb1*((b1 - b)*caa + (-g1 + g)*cbb*saa)*sg1
     &- 4.d0*sb1*sg1 + (2.d0*caa + (g1 - g)*saa*sbb)*(cg1*sa1
     &+ ca1*sb1*sg1) - (2.d0*saa + (-g1 + g)*caa*sbb)*(ca1*cg1
     &- sa1*sb1*sg1))/4.d0
      dR1(2,2)=((-4.d0*a1*cb1 + 4.d0*a*cb1 - 4.d0*g1*cbb*sb1
     &+ 4.d0*g*cbb*sb1 + 5.d0*g1*cbb*sa1*saa*sb1 - 5.d0*g*cbb*sa1*saa
     &*sb1 - 5.d0*g1*cb1*sbb + 5.d0*g*cb1*sbb + 4.d0*g1*cb1*sa1*saa*sbb
     &- 4.d0*g*cb1*sa1*saa*sbb)*sg1 + ca1*saa*(-(g1*cbb*cg1)
     &+ g*cbb*cg1 - 8.d0*cb1*sg1 + 4.d0*b1*sb1*sg1 - 4.d0*b*sb1*sg1)
     &+ caa*(4.d0*(2.d0*cb1*sa1 - b1*sa1*sb1 + b*sa1*sb1
     &+ g1*ca1*cb1*sbb - g*ca1*cb1*sbb)*sg1 + (g1 - g)*cbb
     &*(cg1*sa1 + 5.d0*ca1*sb1*sg1)))/4.d0
      dR1(2,3)=((g1 - g)*cb1*cbb*cg1 + 2.d0*ca1*cb1*cg1*((-g1 + g)
     &*caa*cbb + (-b1 + b)*saa) - 2.d0*cb1*cg1*sa1*((-b1 + b)*caa
     &+ (g1 - g)*cbb*saa) + 2.d0*cg1*sb1*(-a1 + a + (-g1 + g)*sbb)
     &+ cb1*cbb*sg1 - 2.d0*ca1*caa*cb1*cbb*sg1 - 2.d0*cb1*cbb*sa1*saa
     &*sg1 - 2.d0*sb1*sbb*sg1 + 2.d0*(-caa + ((-g1 + g)*saa*sbb)/2.d0)
     &*(-(cg1*sa1*sb1) - ca1*sg1) + 2.d0*(-saa + ((g1 - g)*caa*sbb)
     &/2.d0)*(ca1*cg1*sb1 - sa1*sg1) + caa*sbb*(cg1*sa1 + ca1*sb1*sg1)
     &- saa*sbb*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      dR1(3,1)=(2.d0*((-b1 + b)*co + (-g1 + g)*cbb*so)*sb1 + cb1
     &*(-4.d0 + 2.d0*sa1*saa - g1*caa*sa1*sbb + g*caa*sa1*sbb + ca1
     &*(2.d0*caa + (g1 - g)*saa*sbb)))/4.d0
      dR1(3,2)=(cb1*((g1 - g)*(-4.d0 + 5.d0*co)*cbb + 4.d0*(-b1 + b)
     &*so) + sb1*(4.d0*a1 - 4.d0*a + 8.d0*ca1*saa + 5.d0*g1*sbb
     &- 5.d0*g*sbb - 4.d0*g1*sa1*saa*sbb + 4.d0*g*sa1*saa*sbb
     &- 4.d0*caa*(2.d0*sa1 + g1*ca1*sbb - g*ca1*sbb)))/4.d0
      dR1(3,3)=((-1.d0 + 2.d0*co)*cbb*sb1 + (-2.d0 + co)*cb1*sbb)/2.d0
      dR2(1,1)=(2.d0*cb1*cg1*sa1*((g1 - g)*caa*cbb + (b1 - b)*saa)
     &+ 2.d0*ca1*cb1*cg1*((b1 - b)*caa + (-g1 + g)*cbb*saa)
     &- 4.d0*cg1*sb1 + (2.d0*saa + (-g1 + g)*caa*sbb)*(cg1*sa1*sb1
     &+ ca1*sg1) + (2.d0*caa + (g1 - g)*saa*sbb)*(ca1*cg1*sb1
     &- sa1*sg1))/4.d0
      dR2(1,2)=(-((g1 - g)*cg1*(-2.d0*cbb*sb1 + cbb*sa1*saa*sb1
     &- cb1*sbb + 2.d0*cb1*sa1*saa*sbb)) + ca1*saa*(-4.d0*cb1*cg1
     &- g1*cbb*sg1 + g*cbb*sg1) + caa*(2.d0*cb1*cg1*(2.d0*sa1
     &- g1*ca1*sbb + g*ca1*sbb) - (g1 - g)*cbb
     &*(ca1*cg1*sb1 - sa1*sg1)))/4.d0
      dR2(1,3)=(-((-1.d0 + 2.d0*co)*cb1*cbb*cg1)
     &+ sbb*((-2.d0 + co)*cg1*sb1 - so*sg1))/2.d0
      dR2(2,1)=(-2.d0*cb1*sa1*((g1 - g)*caa*cbb + (b1 - b)*saa)*sg1
     &- 2.d0*ca1*cb1*((b1 - b)*caa + (-g1 + g)*cbb*saa)*sg1
     &+ 4.d0*sb1*sg1 - (2.d0*caa + (g1 - g)*saa*sbb)*(cg1*sa1
     &+ ca1*sb1*sg1) + (2.d0*saa + (-g1 + g)*caa*sbb)
     &*(ca1*cg1 - sa1*sb1*sg1))/4.d0
      dR2(2,2)=((g1 - g)*(-2.d0*cbb*sb1 + cbb*sa1*saa*sb1 - cb1*sbb
     &+ 2.d0*cb1*sa1*saa*sbb)*sg1 + ca1*saa*(-(g1*cbb*cg1) + g*cbb*cg1
     &+ 4.d0*cb1*sg1) + caa*(-2.d0*cb1*(2.d0*sa1 - g1*ca1*sbb
     &+ g*ca1*sbb)*sg1 + (g1 - g)*cbb*(cg1*sa1 + ca1*sb1*sg1)))/4.d0
      dR2(2,3)=-(cb1*cbb*sg1)/2.d0 + ca1*caa*cb1*cbb*sg1
     &+ cb1*cbb*sa1*saa*sg1 + sb1*sbb*sg1 - (caa*sbb*(cg1*sa1
     &+ ca1*sb1*sg1))/2.d0 + (saa*sbb*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      dR2(3,1)=(2.d0*((b1 - b)*co + (g1 - g)*cbb*so)*sb1 + cb1*(4.d0
     &- 2.d0*sa1*saa + g1*caa*sa1*sbb - g*caa*sa1*sbb + ca1
     &*(-2.d0*caa - g1*saa*sbb + g*saa*sbb)))/4.d0
      dR2(3,2)=((g1 - g)*(-2.d0 + co)*cb1*cbb + sb1*(-4.d0*ca1*saa
     &+ g1*sbb - g*sbb - 2.d0*g1*sa1*saa*sbb + 2.d0*g*sa1*saa*sbb
     &+ 2.d0*caa*(2.d0*sa1 - g1*ca1*sbb + g*ca1*sbb)))/4.d0
      dR2(3,3)=((1.d0 - 2.d0*co)*cbb*sb1 - (-2.d0 + co)*cb1*sbb)/2.d0
c
      return
      end
c====================================================================
c
      subroutine dr6(sa1,sb1,sg1,ca1,cb1,cg1,saa,sbb,caa,cbb,
     &               so,co,b1,g1,a1,b,g,a)
      implicit real*8 (a-h,o-z)
c
      dimension dR1(3,3),dR2(3,3)
      call clrmx (3,dR1)
      call clrmx (3,dR2)
c
      dR1(1,1)=(-(g1*cbb*cg1*sa1*saa) + g*cbb*cg1*sa1*saa + 2.d0*cb1*sg1
     &- cbb*sa1*saa*sg1 + b1*sa1*saa*sb1*sg1 - b*sa1*saa*sb1*sg1
     &+ caa*sa1*(b1*cg1 - b*cg1 - cbb*cg1*sb1 + g1*cbb*sb1*sg1
     &- g*cbb*sb1*sg1) + ca1*(-(caa*(g1*cbb*cg1 - g*cbb*cg1 + cbb*sg1
     &- b1*sb1*sg1 + b*sb1*sg1)) + saa*(-(b1*cg1) + b*cg1 + cbb*cg1*sb1
     &- g1*cbb*sb1*sg1 + g*cbb*sb1*sg1)))/2.d0
      dR1(1,2)=(-(cb1*cbb*cg1) + 2.d0*ca1*caa*cb1*cbb*cg1
     &+ 2.d0*cb1*cbb*cg1*sa1*saa + 2.d0*cg1*sb1*sbb + (g1 - g)
     &*cb1*cbb*sg1 + 2.d0*ca1*cb1*((-g1 + g)*caa*cbb + (-b1 + b)*saa)
     &*sg1 - 2.d0*cb1*sa1*((-b1 + b)*caa + (g1 - g)*cbb*saa)*sg1
     &+ 2.d0*sb1*(-a1 + a + (-g1 + g)*sbb)*sg1 - saa*sbb*(cg1*sa1*sb1
     &+ ca1*sg1) + caa*sbb*(-(ca1*cg1*sb1) + sa1*sg1) - (2.d0*saa
     &+ (-g1 + g)*caa*sbb)*(cg1*sa1 + ca1*sb1*sg1) + 2.d0*(-caa
     &+ ((-g1 + g)*saa*sbb)/2.d0)*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      dR1(1,3)=-(cb1*cg1*(-a1 + a + (-g1 + g)*sbb)) + 2.d0*cb1*sbb*sg1
     &+ ((-b1 + b)*caa + (g1 - g)*cbb*saa)*(-(cg1*sa1*sb1) - ca1*sg1)
     &+ ((-g1 + g)*caa*cbb + (-b1 + b)*saa)*(ca1*cg1*sb1 - sa1*sg1)
     &- 2.d0*caa*cbb*(cg1*sa1 + ca1*sb1*sg1)
     &+ 2.d0*cbb*saa*(ca1*cg1 - sa1*sb1*sg1)
      dR1(2,1)=(2.d0*cb1*cg1 + 2.d0*((-g1 + g)*caa*cbb + (-b1 + b)*saa)
     &*(-(cg1*sa1*sb1) - ca1*sg1) - ((g1 - g)*caa*cbb + (b1 - b)*saa)
     &*(cg1*sa1*sb1 + ca1*sg1) - ((b1 - b)*caa + (-g1 + g)*cbb*saa)
     &*(ca1*cg1*sb1 - sa1*sg1) + 2.d0*((-b1 + b)*caa + (g1 - g)*cbb*saa)
     &*(-(ca1*cg1*sb1) + sa1*sg1) - cbb*saa*(cg1*sa1 + ca1*sb1*sg1)
     &- caa*cbb*(ca1*cg1 - sa1*sb1*sg1))/2
      dR1(2,2)=((g1 - g)*cb1*cbb*cg1 + 2.d0*ca1*cb1*cg1*((-g1 + g)
     &*caa*cbb + (-b1 + b)*saa) - 2.d0*cb1*cg1*sa1*((-b1 + b)*caa
     &+ (g1 - g)*cbb*saa) + 2.d0*cg1*sb1*(-a1 + a + (-g1 + g)*sbb)
     &+ cb1*cbb*sg1 - 2.d0*ca1*caa*cb1*cbb*sg1 - 2.d0*cb1*cbb*sa1*saa
     &*sg1 - 2.d0*sb1*sbb*sg1 + 2.d0*(-caa + ((-g1 + g)*saa*sbb)/2.d0)
     &*(-(cg1*sa1*sb1) - ca1*sg1) + 2.d0*(-saa + ((g1 - g)*caa*sbb)
     &/2.d0)*(ca1*cg1*sb1 - sa1*sg1) + caa*sbb*(cg1*sa1
     &+ ca1*sb1*sg1) - saa*sbb*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      dR1(2,3)=2.d0*cb1*cg1*sbb + cb1*(-a1 + a + (-g1 + g)*sbb)*sg1
     &- 2.d0*cbb*saa*(cg1*sa1*sb1 + ca1*sg1) - 2.d0*caa*cbb*(ca1*cg1*sb1
     &- sa1*sg1) + ((-g1 + g)*caa*cbb + (-b1 + b)*saa)*(-(cg1*sa1)
     &- ca1*sb1*sg1) + ((-b1 + b)*caa + (g1 - g)*cbb*saa)
     &*(-(ca1*cg1) + sa1*sb1*sg1)
      dR1(3,1)=(cb1*cbb*so)/2.d0
      dR1(3,2)=((-1.d0 + 2.d0*co)*cbb*sb1 + (-2.d0 + co)*cb1*sbb)/2.d0
      dR2(1,1)=(-2.d0*cb1*sg1 + caa*cbb*(cg1*sa1*sb1 + ca1*sg1)
     &+ cbb*saa*(-(ca1*cg1*sb1) + sa1*sg1) - ((b1 - b)*caa + (-g1 + g)
     &*cbb*saa)*(cg1*sa1 + ca1*sb1*sg1) + ((g1 - g)*caa*cbb
     &+ (b1 - b)*saa)*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      dR2(1,2)=(-(cb1*cbb*cg1) + (g1 - g)*cb1*cbb*sg1 - saa*sbb
     &*(cg1*sa1*sb1 + ca1*sg1) + caa*sbb*(-(ca1*cg1*sb1) + sa1*sg1)
     &+ 2.d0*(saa + ((g1 - g)*caa*sbb)/2.d0)*(cg1*sa1 + ca1*sb1*sg1)
     &+ 2.d0*(caa + ((-g1 + g)*saa*sbb)/2.d0)
     &*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      dR2(1,3)=-(cb1*sbb*sg1) + caa*cbb*(cg1*sa1 + ca1*sb1*sg1)
     &- cbb*saa*(ca1*cg1 - sa1*sb1*sg1)
      dR2(2,1)=(-2.d0*cb1*cg1 - ((g1 - g)*caa*cbb + (b1 - b)*saa)
     &*(cg1*sa1*sb1 + ca1*sg1) - ((b1 - b)*caa + (-g1 + g)*cbb*saa)
     &*(ca1*cg1*sb1 - sa1*sg1) + cbb*saa*(cg1*sa1 + ca1*sb1*sg1)
     &+ caa*cbb*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      dR2(2,2)=((g1 - g)*cb1*cbb*cg1 + cb1*cbb*sg1 + 2.d0*(caa
     &+ ((-g1 + g)*saa*sbb)/2.d0)*(-(cg1*sa1*sb1) - ca1*sg1)
     &+ 2.d0*(saa + ((g1 - g)*caa*sbb)/2.d0)*(ca1*cg1*sb1 - sa1*sg1)
     &+ caa*sbb*(cg1*sa1 + ca1*sb1*sg1)
     &- saa*sbb*(ca1*cg1 - sa1*sb1*sg1))/2.d0
      dR2(2,3)=-(cb1*cg1*sbb) + cbb*saa*(cg1*sa1*sb1 + ca1*sg1)
     &+ caa*cbb*(ca1*cg1*sb1 - sa1*sg1)
      dR2(3,1)=-(cb1*cbb*so)/2.d0
      dR2(3,2)=(-(cbb*sb1) + co*cb1*sbb)/2.d0
c
      return
      end
