c 弾性(曲がり)梁の座屈や大変位問題を、直線梁要素の
c 有限変位・有限要素法で解くプログラム。
c 定式化は、
c 『後藤文彦,小林 裕,斉木 功,岩熊哲夫:
c   空間固定三軸回りの回転自由度を用いた空間梁解析
c   応用力学論文集, 1, 1998, pp.319-327.』
c の式(48)で座標変換をオイラー角で表した接線剛性方程式を用いた。
c 尚、反りに関して、ねじり率の自由度は、
c  『後藤文彦,小林裕,岩熊哲夫:
c   オイラー角を用いた簡潔な有限変位解析手法
c   構造工学論文集, 43A, 1997, pp.333-338.
c の定式化のように扱った。
c
c 02/6/5版
c プログラムを主に作った人:後藤文彦
c 増分式を Mathematica で計算してくれた人:小林裕さん
c いくつかの代数計算やスカイライン法のサブルーチンの
c 元プログラムの作者:岩熊先生(http://www.civil.tohoku.ac.jp/~bear/)
c 
c 改造などは自由ですが、もとのプログラムが、
c http://gthmhk.virtualave.net/programoj/hari/hari.f
c であることをプログラム中に明記して戴けると嬉しいです。
c バグの報告は、
c http://gthmhk.virtualave.net/denbin.html
c にお願いします。
c
         program  hari
        implicit real*8(a-h,o-z)
        dimension nnl(256+1),nnr(256+1),lorg(256*4),lg4(4),
     &  noutl(256),noutr(256),nuvwl(256),nuvwr(256),
     &  f(7*256+8),bc(7*256+8),d(14*256),dd(7*256+8),
     &  ddl(7*256+8),b(7*256+8),alp(256),gmm(256),phi(256),
     &  d14(14),b14(14),dd14(14),sko(7,7),s(14,14),sk(14,14),
     &  stuh(84*256+36),stlh(84*256+36),imx(7*256+9),
     &  msl(7*256+8),nnegtm(7*256+8),anegtm(7*256+8)
         call cl4ivs(256+1,nnl,256+1,nnr,256*4,lorg,4,lg)
         call cl4ivs(256,noutl,256,noutr,256,nuvwl,256,nuvwr)
         call cl4rvs(7*256+8,f,7*256+8,bc,16*256,d,7*256+8,dd)
         call cl4rvs(7*256+8,ddl,7*256+8,b,256,alp,256,gmm)
         call cl4rvs(256,phi, 14,d14, 14,b14, 14,dd14)
         call cl4rms(7,sko, 14,s, 14,sgz, 14,sgx)
         call clrmx(14,sgy)
         call clrmx(14,sk)
         call clrv(84*256+36,stuh)
         call clrv(84*256+36,stlh)
         call cliv(7*256+9,imx)
         call cliv(7*256+8,msl)
         call cliv(7*256+8,nnegtm)
         call clrv(7*256+8,anegtm)
c 一発目の弧長が大きすぎて座屈してしまった場合に、初期状態から
c 計算しなおす時のために、初期変位(ゼロ)をフォーマットなしで fort.11
c に書き出しておく。
        write(11) dd,d,ddl,dp,dpl,ro,p
        rewind 11
c 荷重−変位(u,v,w)関係を出力するファイル:uvw.out
        open(8,file='uvw.out')
c
c ### まずはデータの入力 ###(変数の説明を兼ねる)
c
c  要素数:nofe
        read(*,*) nofe
        write(*,*)'要素数:',nofe
c  座屈モード:mode
        read(*,*) mode
        write(*,*) mode,'番目に低い不安定点を捜します'
c 前に計算したモードに対する倍率:gmode
c 座屈後などを計算する場合は、座屈点付近の接線剛性行列の
c 固有値解析をして得た変位のモードに、倍率 gmode を掛けたものを
c 外乱や(初期)不整として取り込む場合に指定する。
c 不整を与えずに完全形の変位や座屈点などを求める場合は
c gmode=0 を与える。
        read(*,*) gmode
c モード計算のための出力が必要:modeo=1
c モード計算のための出力が不要:modeo=0
c 上述のようなモード計算のために、計算終了時(収束点)の
c 接線剛性行列を出力する場合は modeo=1
        read(*,*) modeo
        nmax=7*256+8  
        n7=7*nofe+7
        n8=7*nofe+8
        call maxa(nofe,n8,imx)
        imx9=imx(n8+1)
        tol=1.d-10
c 初期状態における全ての要素の方向をオイラー角で与える
        do 8 i=1,nofe
         read(*,510) alp(i),gmm(i),phi(i)
8       continue
510     format(1p3d17.9) 
c
c  断面定数
c ea:伸び剛性(ヤング率×断面積)
          read(*,*) ea
c eix:x軸回りの曲げ剛性(ヤング率×断面二次モーメント)
          read(*,*) eix
c eiy:y軸回りの曲げ剛性(ヤング率×断面二次モーメント)
          read(*,*) eiy
c gj:ねじり剛性:(せん断弾性係数×ねじり定数)
          read(*,*) gj
c eiw:そりねじり剛性(ヤング率×そりねじり定数)
          read(*,*) eiw
c el:一要素の長さ
        read(*,*) el
        ell=el*nofe
        th=-alp(nofe)*2.d0
        pi=asin(1.d0)*2.d0
c 比較対照解の種類:nhikak
c nhikak=1:圧縮単純梁のオイラー座屈の解
c nhikak=2:等曲げ円弧アーチの横ねじれ座屈のウラソフの解
c nhikak=3: 片持ち梁の横ねじれ座屈のトラヘアの解(面内変位無視)
        read(*,*) nhikak
          goto (511,512,513), nhikak
511    continue
c 柱のオイラー座屈公式
         eul=(pi/ell)**2*eiy
         write(*,*) 'Euler=',eul
         goto 599
512    continue
c 等曲げを受ける円弧アーチの横ねじれ座屈のウラソフの解
         eix2=eix*1.d6
         call vlasov(mode,ell,eix2,eiy,gj,eiw,th,pi,vln,vlp)
         write(*,*)'Vlasov(負)=',vln
         write(*,*)'Vlasov(正)=',vlp
c と修正ウラソフの解(座屈前の面内変位の影響考慮)
         call vlasov(mode,ell,eix,eiy,gj,eiw,th,pi,vln,vlp)
         write(*,*)'修正Vlasov(負)=',vln
         write(*,*)'修正Vlasov(正)=',vlp
          goto 599
513    continue
c 片持ち梁の横ねじれ座屈のトラヘアの解(座屈前面内変位無視)
         a1=sqrt(eiy*gj)/ell/ell
         a2=sqrt(eiw/gj)*pi/ell
         tra=a1*(3.95d0+3.52d0*a2)
         write(*,*) 'Trahair=',tra
           goto 599
599    continue
c
c
c  ###  載荷条件  ###
c 載荷されている節点の数:nofln
         read(*,*) nofln
c 載荷節点の番号:nnl
c 載荷節点における変位ベクトル:f(自由度数+1)
        do 10 i=1,nofln
         read(*,*) nnl(i)
        do 20 j=1,7
         nf=7*(nnl(i)-1)+j
         read(*,*) f(nf)
c 荷重は el と eiy で無次元化
c これに対応して変位も無次元量で求まる。
          if(j.lt.4) f(nf)=f(nf)*el*el/eiy
          if(j.gt.3) f(nf)=f(nf)*el/eiy
          if(j.eq.7) f(nf)=f(nf)/eiy
20      continue
10      continue
c
c  ###  境界条件 ####
          do 30 i=1,n8
            bc(i)=1.d0
30        continue
c 拘束されている節点の数:nofrn
         read(*,*) nofrn
c 拘束節点の番号:nnr
c と拘束節点における変位ベクトル:bc(自由度数+1)
        do 40 i=1,nofrn
         read(*,*) nnr(i)
           do 50 j=1,7
            nb=7*(nnr(i)-1)+j
            read(*,*) bc(nb)
50         continue
40      continue
c 境界条件を局所座標系で与える(節点変位をそこだけ局所系で
c 定義して全体系に変換する座標変換を掛けないでおく)要素の数:noflbc
c (円弧アーチの横倒れ座屈を解く際に、端部の面外回転を半径軸回りに
c 与える場合など)
        read(*,*) noflbc
c その要素番号を入れ、次の四項目の境界条件をどう与えるのか
c (局所系1/全体系0)を lorg(要素数*4) に入れる。四項目とは、
c 左の並進、左の回転、右の並進、右の回転
         do 60 i=1,noflbc
          read(*,*) lbcen
           l1=lbcen-1
       read(*,500) lorg(l1*4+1),lorg(l1*4+2),lorg(l1*4+3),lorg(l1*4+4)
60       continue
500    format(4i2)
c
c  ### 出力条件 ###
c
c 左の節点変位を出力したい要素の要素数:noflo, その要素番号:noutl
        read(*,*) noflo
         do 70 i=1,noflo
          read(*,*) noutl(i)
70       continue
c 右の節点変位を出力したい要素の要素数:nofro, その要素番号:noutr
        read(*,*) nofro
         do 80 i=1,nofro
          read(*,*) noutr(i)
80       continue
c 左の節点並進変位をファイル出力したい要素の要素数:noflof
c その要素番号: nuvwl
        read(*,*) noflof
         do 75 i=1,noflof
          read(*,*) nuvwl(i)
75       continue
c 右の節点並進変位をファイル出力したい要素の要素数:nofrof
c その要素番号: nuvwr
        read(*,*) nofrof
         do 85 i=1,nofrof
          read(*,*) nuvwr(i)
85       continue
c
c 一番最初の一発目の増分ステップの弧長:r
        read(*,*) r
c 2ステップ目以降の計算に使う弧長増分:rn
        read(*,*) rn
c
c
c ############# ここまでがデータ入力 #############
c
c 収束精度:収束時の荷重が、N-R反復における一つ前のステップの
c 荷重と有効数字何桁まで合うようにしたいかに応じて cnvr の小ささを
c 決める
        cnvr=5.d-8
c 外乱の大きさに応じて収束精度を甘くしたいときとかは、
cc       if(gmode.gt.1.d-10) cnvr=1.d-3
c みたいなのを入れたりするとか。
c
c 三次元梁の線形剛性行列(反りを含む1節点7自由度)を s(14,14)
c に入れる
        call smx(el,ea,eix,eiy,gj,eiw,s)
c
c ###  初期条件 ### 
c
c 変位増分ベクトル:dd(自由度数+1)
c 変位増分の回転角成分は空間固定3軸回りの微少回転角
c 変位ベクトル:d(節点自由度数*要素数)
c トータルな変位は、要素どうしが繋がる2節点のオイラー角を
c 別々に覚えておく。
c ddl は dd を各増分ステップの繰り返しで総和したもの
c dp は荷重増分(弧長法なので、増分変位と同様に解かれる未知変数)
        do 100 i=1,n7
         dd(i)=0.d0
         d(i)=0.d0
         ddl(i)=0.d0
100     continue
         dp=0.d0
         p=0.d0
         ro=r
         dpl=r
c まだ一度も収束値が求まっていないことを表す:nccl(後述)
         nccl=0
c  ### 各荷重ステップの一発目の弧長増分過程 ###
1000   continue        
           nofcs=0
c 計算開始からの増分繰り返し数:ninc
        ninc=ninc+1
         if(ninc.gt.4000) then
           write(*,*)'(増分繰り返しが 4000 を越えました)'
           close(8)
           stop
         endif
c
c 全変位増分ddと全変位dの中から、1要素ぶんずつをそれぞれ
c dd14とd14に抜き出す
       do 1100 i=0,nofe-1
        do 120 j=1,14
         dd14(j)=dd(i*7+j)
         d14(j)=d(i*14+j)
120     continue
         do 130 j=1,4
          lg4(j)=lorg(i*4+j)
130      continue
c dd14とd14を使って1要素ずつ接線剛性行列を計算していく
c 求まった1要素ぶんの接線剛性行列はskに入る
        call stmx(nofe,i,alp,gmm,dd14,d14,sk,b14,lg4,s)
         do 140 j=1,14
          d(i*14+j)=d14(j)
140      continue
c sk をスカイライン法の1次元ベクトルに書き込む
c stuh に上三角、stlh に下三角が入る
        call skyln(nofe,i,sk,sko,stuh,stlh,imx)
1100   continue
c
        do 150 i=1,n7
         i8=imx(n8)+n8-i
         stuh(i8)=-f(i)
         stlh(i8)=ddl(i)/ro
         ddl(i)=0.d0
         b(i)=0.d0
150     continue
        do 160 i=1,n8
         stlh(imx(i))=b(i)
160     continue
         stuh(imx(n8))=dpl/ro
         stlh(imx(n8))=r
         ro=r
         dpl=0.d0
c 境界条件を入れる
        call bncn(nofe,bc,stuh,stlh,imx)
c 代数方程式をガウスの消去法で解いて変位ベクトルを求め、
c その際に対角項の符号から、行列の負の固有値の数 jdt を求める
        call solve(stuh,stlh,dp,imx,msl,nnegtm,anegtm,
     &                   det,idt,jdt,n8)
        do 170 i=1,n8
         dd(i)=stlh(imx(i))
170     continue
c
c  ### ニュートン−ラフソン 収束過程 ###
c(基本的に上の「一発目の弧長増分過程」の計算と同じ)
2000   continue
c 各増分ステップでの収束までの繰り返し数:nofcs
        nofcs=nofcs+1
        do 180 i=1,n7
         ddl(i)=ddl(i)+dd(i)
         ddi2=ddi2+dd(i)**2
         ddl2=ddl2+ddl(i)**2
180     continue
         dpl=dpl+dp
         p=p+dp
         arc=ddi2/ddl2
         ddi2=0
         ddl2=0
       if(nofcs.gt.6) then
        write(*,*)'反復数が6を越えても収束しないので弧長を半分にする'
c 反復数が6を越えて辛うじて収束するようなことが何回も重なると、
c 答えがおかしくなってたりするので、その回数を nlet で数えて、
c 適当な回数に達したら強制的に終了して答えを出力するように
c 荷重が収束した(p1=p)したことにしてやる。
        nlet=nlet+1
          if(nlet.gt.10) then
           p1=p
           goto 4004
          endif
        goto 4001
       endif
        if(arc.lt.tol) goto 3000
        do 190 i=1,n7
        b(i)=p*f(i)
190     continue
        do 1200 i=0,nofe-1
        do 200 j=1,14
        dd14(j)=dd(i*7+j)
        d14(j)=d(i*14+j)
200     b14(j)=b(i*7+j)
        do 210 j=1,4
         lg4(j)=lorg(i*4+j)
210     continue
c 接線剛性行列を1要素ずつ計算してskに入れる
        call stmx(nofe,i,alp,gmm,dd14,d14,sk,b14,lg4,s)
        do 220 j=1,14
        d(i*14+j)=d14(j)
220     continue
c skをスカイライン法の1次元配列に書き込む
c stuhが上三角行列、stlhが下三角行列
        call skyln(nofe,i,sk,sko,stuh,stlh,imx)
        do 230 j=1,14
        b(i*7+j)=b14(j)
230     continue
1200    continue
        do 240 i=1,n7
        i8=imx(n8)+n8-i
        stuh(i8)=-f(i)
        stlh(i8)=ddl(i)
240     continue
        do 250 i=1,n8
        stlh(imx(i))=b(i)
250     continue
        stuh(imx(n8))=p
        stlh(imx(n8))=0.d0
c 境界条件を入れる
        call bncn(nofe,bc,stuh,stlh,imx)
c 座屈点の探査のとき、荷重が前の増分ステップと誤差範囲に
c になったら終了する。モード計算用の出力が必要な場合(modeo=1)は
c fort.7 に収束時の接線剛性行列をフォーマットなしで書き出す。
       if( (abs((p1-p)/p).lt.cnvr).and.(gmode.lt.1.d-20) ) then
            if(modeo.eq.1) then
              write(7) nofe+1,stuh,stlh,imx,imx9
              rewind 7
            endif
          close(8)
c
c 出力した接線行列の固有値解析をせずに、最後のN-R反復で求まる
c (微少な)変位増分を固有ベクトル(座屈モード)の近似として
c 利用する場合のために、一応、近似モードもここで出力しておく。
c こういう近似ができるということは、
c 久田俊明・野口裕久 著『非線形有限要素法の基礎と応用』
c(丸善株式会社)p.317とかに書いてある。
c
          call modeout(nofe,dd)
c
c 但し、接線剛性方程式の回転自由度は空間固定3軸回りの成分なので、
c 座屈モードの回転成分もその成分に対応して求まる。
c 座屈後へ分岐させる外乱としては、3軸回り成分をオイラー角に与えても
c 十分に用を足すとは思うが。
c オイラー角成分がほしいときは、変換行列をかけるとかして。
c
          stop
        endif
c
c
c 接線剛性方程式を解く
        call solve(stuh,stlh,dp,imx,msl,nnegtm,anegtm,
     &                   det,idt,jdt,n8)
        do 260 i=1,n8
        dd(i)=stlh(imx(i))
260     continue
        goto 2000
c
c  ### 出力 ### 
c
3000   continue
           arc=0.d0
c 収束しないために弧長が小さくなり過ぎたのを、また大きくしたい
c ときは、nlet の大きさに応じて弧長を大きくしてやるとか
c           if(nlet.gt.1) rn=rn*(2.d0)**nlet
c ここを通過しているということは、増分ステップが収束したという
c ことなので増分ステップ内の繰り返し数 nlet を0にする。
            nlet=0
c 座屈点の探査は、座屈点を超えたら、越える前の荷重ステップまで
c 戻り、弧長を半分にして進み、また行き過ぎたらまた戻って弧長を
c 半分にして次の荷重ステップを計算して座屈点に近づいていく。
c この座屈点を跨いだいったりきたりの計算において、同じ点を重複して計算
c しないように、「アキレスと亀」方式で弧長を毎回半分にしていく。
c
c 接線剛性行列の負の固有値数 jdt と固有ベクトルとして得られる
c 座屈モードの次数とはだいたい対応していることから、jdt が
c 指定した座屈モードの次数(mode)を越えたかどうかで座屈点を
c 越えたかどうかを判断する。座屈点を超えたかどうかにも
c 様々な場合があり得るので、その状況をnacl,nbcl,ncclを使って判断する。
c
c 一発目の弧長を大きく取りすぎて、最初から座屈点を超えてしまった場合
c は、正しい答えがまだ求まっていないので、最初から計算をやり直す
c 必要がある。そのため、まだ一つも収束値が求まっていない時は
c nccl=0, 一つでも収束値が(座屈点を超える前に)求まったら
c nccl=1 とする。
        if(jdt.lt.mode) nccl=1
        write(*,*) 'N-R反復数=',nofcs
        write(*,570)       det,    idt,jdt
570     format(' ','det=',f11.3,'d',i4,'  ,  ','負の固有値数:',i4)
          goto (571,572,573), nhikak
571    write(*,581) ' p=',p,'   PL/EA=',p*ell/ea,'    Euler=',eul
       goto 579
572    write(*,581) ' p=',p,'   ML^2/8EI=',p*ell**2/8.d0/eix
       write(*,580) ' 修正Vlasov(正)=',vlp,'   修正Vlasov(負)=',vln
       goto 579
573    write(*,581)' p=',p,'   PL^3/3EI=',p*ell**3/3./eix,'   Tra=',tra
       goto 579
579     continue
c
c        write(*,580)' p=',p,'            Ml/EIy=',p*ell/e/yi
580     format(' ',a,1pd13.5,a,1pd13.5)
581     format(' ',a,1pd13.5,a,1pd13.5,a,1pd13.5)
c 節点変位は、u, v, w, α, γ, φ, λの順に出力。
c α, γ, φ はオイラー角、λはねじれ率
c d は無次元化されているので、出力するときは有次元化する
c
c 指定要素の左節点を出力する場合
        if(noflo.ne.0) then
         do 460 i=1,noflo
          j=noutl(i)-1
       write(*,590) noutl(i),'左',d(j*14+1)*el,d(j*14+2)*el,
     & d(j*14+3)*el,d(j*14+4),d(j*14+5),d(j*14+6),d(j*14+7)/el
460      continue 
        endif
c 指定要素の右節点を出力する場合
        if(nofro.ne.0) then
         do 470 i=1,nofro
          j=noutr(i)-1
       write(*,590) noutr(i),'右',d(j*14+8)*el,d(j*14+9)*el,
     & d(j*14+10)*el,d(j*14+11),d(j*14+12),d(j*14+13),d(j*14+14)/el
470      continue 
        endif
c 指定要素の左節点のuvwをファイル出力
        if(noflof.ne.0) then
         do 480 i=1,noflof
          j=nuvwl(i)-1
       write(8,560) p,d(j*14+1)*el,d(j*14+2)*el,d(j*14+3)*el
480      continue 
        endif
c 指定要素の右節点のuvwをファイル出力
        if(nofrof.ne.0) then
         do 490 i=1,nofrof
          j=nuvwr(i)-1
       write(8,560) p,d(j*14+8)*el,d(j*14+9)*el,d(j*14+10)*el
490      continue 
        endif
560      format(1p4d11.3)
590      format(i3,a,1p6d11.3,1pd10.2)
        if(mode.eq.0) goto 4003
c 接線剛性行列の負の固有値数jdtが指定した座屈モード次数を
c 越えたら、座屈したと見なす。
        if(jdt.ge.mode) then
         write(*,*)'座屈したみでだよ................'
c 最初に座屈点を超えたら、まずnbcl=1,nacl=0としておく。
         nbcl=1
         nacl=0
        endif
4004    continue
c
c ### 外乱の入れ方 ###
c 一発目の弧長増分ステップが収束した後に入れる。
c 予め座屈点を計算しておいて、一発目の弧長を座屈点の直前の
c 分岐しやすい荷重で外乱が入るように調整する。
            if (nccl.eq.1) then
c
c 「アキレスと亀」方式で座屈点に達したら、そこで外乱を読み込む
c 場合は、
c       if(abs((p1-p)/p).lt.1.d-1) then
c とか。
        if(abs(gmode).lt.1.d-20) goto 5002
c 但し、外乱は一度しか読み込まない
        if(iyonda.eq.1) goto 5002
         write(*,*)'外乱読み込み中'
c 外乱を読み込んだ後に収束精度を甘くしたいときは、
c         cnvr=1.d-3
c とか。外乱を1度 読んだらiyonda=1
          iyonda=1
c
c
           call modein(nofe,d,gmode)
c
c
          rn=rn*1.5d0
             endif
5002    continue
        p1=p
c 座屈点を超えたために座屈前に戻って1ステップ以上、座屈点を超えない
c 点を解いているとき(nbcl=0, nacl=1)は
c 弧長は半分にするけど、座屈してない状態だから前のステップには
c 戻らなくていい。rn=rn/2の後、4002まで飛んで1000に戻る。
       if(nacl.eq.1) rn=rn/2.d0
c 座屈点を超えたばかりの状態(nbcl=1, nacl=0)のときは、
c 座屈前に戻って弧長を半分にする(それでも、また座屈した時は
c 上でnbcl=1になってるので、またここを通過する)。
c 次にここを通ときに座屈前に戻っていることを教えるために
c nbcl=0, nacl=1 とする(nbcl=nacl=1の組み合わせはない)
       if(nbcl.eq.1) then
           nbcl=0
           nacl=1
         goto 4001
       endif
c こういうふうに goto が交差した書き方をしたのは、
c 一旦ここでif文を閉じておいて、上から4001へ飛べるようにするため。
         goto 4002
4001    continue
           read(11) dd,d,ddl,dp,dpl,ro,p
           rewind 11
           rn=rn/2.d0
c 一発目の弧長が大きすぎて座屈してしまった場合(nccl=0)は、
c 弧長を一割小さくして初期状態から計算しなおす。
           if(nccl.eq.0) then
            r=r*0.9d0
            ro=r
            dpl=ro
           endif
4002    continue
c 一度でも収束値が求まっている(nccl=1)なら、収束するたびに、
c 変位をfort.11にフォーマットなしで書き出す。但し、
c 座屈点を超えた場合は、上記のread(11)で、座屈前の変位を読み込んだ
c ものが、ここで書き出される。つまり、ここでは常に座屈していない
c 状態での最も後に計算した変位が出力される。そして、次に座屈点を
c 越えたときに、この変位の点まで戻れるようにしておく。
         if(nccl.eq.1) then
            write(11) dd,d,ddl,dp,dpl,ro,p
            rewind 11
            r=rn
         endif
4003    continue
       write(*,*)'----------------------------------------------------'
c ここまでで一つの弧長増分ステップの終わり。次の弧長増分ステップ
c に戻る。
        goto 1000
        end
c
c ものすごく長いけど、ここまでがメインプログラム。
c 以下がサブルーチン
c
c 等曲げを受ける円弧アーチの横ねじれ座屈に対するウラソフの解析解
        subroutine vlasov(n,blen,eix,eiy,gj,eiw,th,pi,vln,vlp)
        implicit real*8(a-h,o-z)
        if(th.lt.1.d-6) th=1.d-6
        r=blen/th
        f=(n*pi/blen)**2*eiw+gj
        a=eiy/eix-1.d0+(1.d0/eix-eiy/eix/eix)*f
        b=(eiy+(1.d0-2.d0*eiy/eix)*f)/r
        c=((n*pi/blen)**2-1.d0/r/r)*eiy*f
        d=b**2-4.d0*a*c
        vln=(-b+sqrt(d))/2.d0/a
        vlp=(-b-sqrt(d))/2.d0/a
        return
        end
c
c 梁の線形剛性行列(無次元化してある)
c 2節点14自由度
        subroutine smx(el,ea,eix,eiy,gj,eiw,s)
        implicit real*8(a-h,o-z)
        dimension s(14,14)
        s(3,3)=   ea*el*el/eiy
        s(10,3)= -ea*el*el/eiy
        s(10,10)= ea*el*el/eiy
        s(2,2)=  12.d0*eix/eiy
        s(4,2)=  -6.d0*eix/eiy
        s(4,4)=   4.d0*eix/eiy
        s(11,2)= -6.d0*eix/eiy
        s(11,4)=  2.d0*eix/eiy
        s(1,1)=  12.d0
        s(5,1)=   6.d0
        s(5,5)=   4.d0
        s(12,1)=  6.d0
        s(12,5)=  2.d0
        s(6,6)=  12.d0*eiw/eiy/el/el+gj/eiy*6.d0/5.d0
        s(13,6)=-12.d0*eiw/eiy/el/el-gj/eiy*6.d0/5.d0
        s(7,6)=  -6.d0*eiw/eiy/el/el-gj/eiy/10.d0
        s(14,6)= -6.d0*eiw/eiy/el/el-gj/eiy/10.d0
        s(7,7)=   4.d0*eiw/eiy/el/el+gj/eiy*2.d0/15.d0
        s(14,7)=  2.d0*eiw/eiy/el/el-gj/eiy/30.d0
        s(14,14)= 4.d0*eiw/eiy/el/el+gj/eiy*2.d0/15.d0
        s(9,2)=-s(2,2)
        s(9,9)=s(2,2)
        s(9,4)=-s(4,2)
        s(11,9)=s(9,4)
        s(11,11)=s(4,4)
        s(8,1)=-s(1,1)
        s(8,5)=-s(5,1)
        s(8,8)=s(1,1)
        s(12,8)=-s(5,1)
        s(12,12)=s(5,5)
        s(13,7)=-s(7,6)
        s(13,13)=s(6,6)
        s(14,13)=s(13,7)
c
        do 100 i=1,13
        do 100 j=i+1,14
        s(i,j)=s(j,i)
100     continue
c
        return
        end
c
c 接線剛性行列(変位に関する増分式)の計算
c
      SUBROUTINE STMX(nofe,is,alp,gmm,dd,d,sl,b,lg,ss)
        implicit real*8 (a-h,o-z)
        DIMENSION ss(14,14),sk(14,14),sl(14,14),d(14),dl(14),dd(14),
     & t(14,14),ti(14,14),t14(14,14),eti(14,14),alp(256),gmm(256),
     & cf1(7,7),cf2(7,7),t1(7,7),ta(7,7),tg(7,7),tp(7,7),b(15),lg(4),
     & deti(14,14,14),r7x14(7,14),dr7x14(14,7,14),dss(14,14,14),
     & dt0(7,14,14),dmy1(14),dmy2(14),
     & tngnt1(14,14),tngnt3(14,14),tngnt4(14,14),sngnt4(14,14)
       call cl4rms(14,t14, 14,sl, 14,sk, 14,ti)
       call cl4rms(7,ta, 7,tg, 7,tp, 7,cf1)
       call cl4rms(7,cf2, 14,eti, 7,t1, 14,t)
       call clrv(14,dl)
       call cl4rms(14,tngnt1,14,tngnt3,14,tngnt4,14,sngnt4)
       call cl5rms(deti,r7x14,dr7x14,dss,dt0)
       call clrmx(14,sngnt4)
c
c 曲がり梁の要素の初期方向を与える座標変換:ti
        call trns(ti,alp(is+1),gmm(is+1),0.d0)
c 3軸回りの回転成分をオイラー角成分に変換する行列:eti
        call euler(eti,d(4),d(5),d(11),d(12))
c 端部の境界条件が局所系で与えられていた時は、その変位に
c 掛けられるtiを外す(単位行列に変える)
c 左端の並進変位
       if(lg(1).eq.1) then
         do 1 i=1,3
         do 1 j=1,3
           ti(i,j)=0.d0
           ti(i,i)=1.d0
1        continue
       endif
c 左端の回転変位
       if(lg(2).eq.1) then
         do 2 i=4,6
         do 2 j=4,6
           ti(i,j)=0.d0
           ti(i,i)=1.d0
2        continue
       endif
c 右端の並進変位
       if(lg(3).eq.1) then
         do 3 i=8,10
         do 3 j=8,10
           ti(i,j)=0.d0
           ti(i,i)=1.d0
3        continue
       endif
c 右端の回転変位
       if(lg(4).eq.1) then
         do 4 i=11,13
         do 4 j=11,13
           ti(i,j)=0.d0
           ti(i,i)=1.d0
4        continue
       endif
c 変位増分ddをti系のオイラー角成分に変換する
c dmy1=ti^T * dd    (^T は転置)
       call mxtv(14,ti,dd,dmy1)
c dd=eti * dmy1 = eti * ti^T * dd
       call mxv(14,eti,dmy1,dd)
c オイラー角成分にしてから足し合わせる
       do 80 i=1,14
       d(i)=d(i)+dd(i)
80     continue
c オイラー角が更新されたところで、座標変換行列を再定義:t
       call trns(t,d(4),d(5),d(6))
c 3軸成分をオイラー角成分に変換する行列も再定義:eti
       call euler(eti,d(4),d(5),d(11),d(12))
c 変位ベクトル
        u1=d(1)
        v1=d(2)
        w1=d(3)
        a1=d(4)
        g1=d(5)
        p1=d(6)
        q1=d(7)
        u2=d(8)
        v2=d(9)
        w2=d(10)
        a2=d(11)
        g2=d(12)
        p2=d(13)
        q2=d(14)
        ca1=cos(a1)
        cg1=cos(g1)
        cp1=cos(p1)
        ca2=cos(a2)
        cg2=cos(g2)
        cp2=cos(p2)
        sa1=sin(a1)
        sg1=sin(g1)
        sp1=sin(p1)
        sa2=sin(a2)
        sg2=sin(g2)
        sp2=sin(p2)
c
c 相対変位
        dl(7)=q1
        dl(8)=u2-u1-ca1*sg1
        dl(9)=v2-v1-sa1
        dl(10)=w2-w1-(ca1*cg1-1.d0)
        gg=(g1+g2)/2.d0
        aa=(a1+a2)/2.d0
        pp=(p1+p2)/2.d0
        cpp=cos(pp)
        caa=cos(aa)
        cgg=cos(gg)
        spp=sin(pp)
        saa=sin(aa)
        sgg=sin(gg)
        dl(11)=-cgg*(a2-a1)+sgg*caa*(p2-p1)
        dl(12)=(g2-g1)+saa*(p2-p1)
        dl(13)=sgg*(a2-a1)+cgg*caa*(p2-p1)
        dl(14)=q2
c
c t の増分
        ta(1,1)=-sp1*ca1*sg1
        ta(1,2)=-cp1*ca1*sg1
        ta(1,3)=-sa1*sg1
        ta(2,1)=-sp1*sa1
        ta(2,2)=-cp1*sa1
        ta(2,3)=ca1
        ta(3,1)=-sp1*ca1*cg1
        ta(3,2)=-cp1*ca1*cg1
        ta(3,3)=-sa1*cg1
        tg(1,1)=-cp1*sg1-sp1*sa1*cg1
        tg(1,2)=sp1*sg1-cp1*sa1*cg1
        tg(1,3)=ca1*cg1
        tg(3,1)=-cp1*cg1+sp1*sa1*sg1
        tg(3,2)=sp1*cg1+cp1*sa1*sg1
        tg(3,3)=-ca1*sg1
        tp(1,1)=-sp1*cg1-cp1*sa1*sg1
        tp(1,2)=-cp1*cg1+sp1*sa1*sg1
        tp(2,1)=cp1*ca1
        tp(2,2)=-sp1*ca1
        tp(3,1)=sp1*sg1-cp1*sa1*cg1
        tp(3,2)=cp1*sg1+sp1*sa1*cg1
        do 162 i=1,3
        do 162 j=1,3
        ta(3+i,3+j)=ta(i,j)
        tg(3+i,3+j)=tg(i,j)
        tp(3+i,3+j)=tp(i,j)
162     continue
c
        do 100 i=1,6
        do 100 j=1,6
           dt0(4,i,j)=ta(i,j)
           dt0(4,i+7,j+7)=ta(i,j)
           dt0(5,i,j)=tg(i,j)
           dt0(5,i+7,j+7)=tg(i,j)
           dt0(6,i,j)=tp(i,j)
           dt0(6,i+7,j+7)=tp(i,j)
100     continue
c
c 相対変位の増分
        cf1(1,1)=-1.d0
        cf1(1,4)=sa1*sg1
        cf1(1,5)=-ca1*cg1
        cf1(2,2)=-1.d0
        cf1(2,4)=-ca1
        cf1(3,3)=-1.d0
        cf1(3,4)=sa1*cg1
        cf1(3,5)=ca1*sg1
        cf1(4,4)=cgg-(p2-p1)/2.d0*sgg*saa
        cf1(4,5)=(a2-a1)/2.d0*sgg+(p2-p1)/2.d0*cgg*caa
        cf1(4,6)=-sgg*caa
        cf1(5,4)=(p2-p1)/2.d0*caa
        cf1(5,5)=-1.d0
        cf1(5,6)=-saa
        cf1(6,4)=-sgg-(p2-p1)/2.d0*cgg*saa
        cf1(6,5)=(a2-a1)/2.d0*cgg-(p2-p1)/2.d0*sgg*caa
        cf1(6,6)=-cgg*caa
        cf2(1,1)=1.d0
        cf2(2,2)=1.d0
        cf2(3,3)=1.d0
        cf2(4,4)=-cgg-(p2-p1)/2.d0*sgg*saa
        cf2(4,5)=cf1(4,5)
        cf2(4,6)=sgg*caa
        cf2(5,4)=cf1(5,4)
        cf2(5,5)=1.d0
        cf2(5,6)=saa
        cf2(6,4)=cf1(6,4)+2.d0*sgg
        cf2(6,5)=cf1(6,5)
        cf2(6,6)=cgg*caa
        cf2(7,7)=1.d0
c
c その他、煩雑な増分式は、farejo の中で作る
       call farejo(d,dl,t,cf1,cf2,dt0,r7x14,dr7x14,deti)
c
        call mxtv8(t,dl,dmy1)
        call mxv8(ss,dmy1,dmy2)
       do 140 i=1,14
       do 140 j=1,14
          tngnt1(i,j)=0.d0
       do 140 k=1,7
c       dr7x17 no naka no sori no zoubun ha 0 
          tngnt1(i,j)=tngnt1(i,j)+dr7x14(j,k,i)*dmy2(k+7)
140    continue
c
        call mxm14(ss,r7x14,sngnt4)
        call mtmx14(r7x14,sngnt4,tngnt3)
c 
        do 260 i=1,14
        do 260 j=1,14
           sl(i,j)=tngnt1(i,j)+tngnt3(i,j)
260     continue
c
       call mxmx(14,sl,eti,sk)
       call mxtmx(14,eti,sk,sl)
c
       call mtv14(r7x14,dmy2,dmy1)
c
       do 300 i=1,14
       do 300 j=1,14
          tngnt4(i,j)=0.d0
       do 300 k=1,14
          tngnt4(i,j)=tngnt4(i,j)+deti(j,k,i)*dmy1(k)
300    continue
c
       call mxmx(14,tngnt4,eti,sngnt4)
c
        do 290 i=1,14
        do 290 j=1,14
           sk(i,j)=sngnt4(i,j)+sl(i,j)
290     continue
c
       call mxmxt(14,sk,ti,sngnt4)
       call mxmx(14,ti,sngnt4,sl)
c
c 不釣り合い力の計算
c
       call mxtv(14,t,dl,dmy2)
       call mxv(14,ss,dmy2,dmy1)
       call mtv14(r7x14,dmy1,dmy2)
       call mxtv(14,eti,dmy2,dmy1)
       call mxv(14,ti,dmy1,dmy2)
       do 550 i=1,14
       b(i)=b(i)-dmy2(i)
550    continue
       return
c888    format(14f5.1)
       end
c
c オイラー角で座標変換行列を定義する
        subroutine trns(t,al,gm,ph)
        IMPLICIT REAL*8 (A-H,O-Z)
        dimension t(14,14)
c       write(*,*)'I am in sub trns'
        t(1,1)=cos(ph)*cos(gm)-sin(ph)*sin(al)*sin(gm)
        t(1,2)=-sin(ph)*cos(gm)-cos(ph)*sin(al)*sin(gm)
        t(1,3)=cos(al)*sin(gm)
        t(2,1)=sin(ph)*cos(al)
        t(2,2)=cos(ph)*cos(al)
        t(2,3)=sin(al)
        t(3,1)=-cos(ph)*sin(gm)-sin(ph)*sin(al)*cos(gm)
        t(3,2)=sin(ph)*sin(gm)-cos(ph)*sin(al)*cos(gm)
        t(3,3)=cos(al)*cos(gm)
        do 10 i=1,3
        do 10 j=1,3
        t(i+3,j+3)=t(i,j)
        t(i+7,j+7)=t(i,j)
        t(i+10,j+10)=t(i,j)
10      continue
        t(7,7)=1.d0
        t(14,14)=1.d0
        return
        end
c
c 3軸回りの回転角成分をオイラー角成分に変換する
        subroutine euler(t,a1,g1,a2,g2)
        IMPLICIT REAL*8 (A-H,O-Z)
        dimension t(14,14)
        do 10 i=1,14
        t(i,i)=1.d0
10      continue
        t(4,4)=-cos(g1)
        t(4,6)=sin(g1)
        t(5,4)=-sin(g1)*tan(a1)
        t(5,6)=-cos(g1)*tan(a1)
        t(6,4)=sin(g1)/cos(a1)
        t(6,6)=cos(g1)/cos(a1)
        t(11,11)=-cos(g2)
        t(11,13)=sin(g2)
        t(12,11)=-sin(g2)*tan(a2)
        t(12,13)=-cos(g2)*tan(a2)
        t(13,11)=sin(g2)/cos(a2)
        t(13,13)=cos(g2)/cos(a2)
        return
        end
c
c
c 以下の mx...とかmt...という一連のサブルーチンは行列の
c かけ算を行うためのもの。
c
       subroutine mtmx14(a,b,c)
       implicit real*8(a-h,o-z)
       dimension a(7,14),b(14,14),c(14,14)
       do 10 i=1,7
       do 10 k=1,7
        c(i,k)=0.d0
       do 10 j=1,7
        c(i,k)=c(i,k)+a(j,i)*b(j+7,k)
        c(i,k+7)=c(i,k+7)+a(j,i)*b(j+7,k+7)
        c(i+7,k)=c(i+7,k)+a(j,i+7)*b(j+7,k)
        c(i+7,k+7)=c(i+7,k+7)+a(j,i+7)*b(j+7,k+7)
10     continue
       do 20 j=1,14
          c(7,j)=c(7,j)+b(7,j)
20    continue
       return
       end
c
c
       subroutine mxm14(a,b,c)
       implicit real*8(a-h,o-z)
       dimension a(14,14),b(7,14),c(14,14)
       do 10 i=1,7
       do 10 k=1,7
        c(i,k)=0.d0
       do 10 j=1,7
        c(i,k)=c(i,k)+a(i,j+7)*b(j,k)
        c(i,k+7)=c(i,k+7)+a(i,j+7)*b(j,k+7)
        c(i+7,k)=c(i+7,k)+a(i+7,j+7)*b(j,k)
        c(i+7,k+7)=c(i+7,k+7)+a(i+7,j+7)*b(j,k+7)
10     continue
       do 20 i=1,14
          c(i,7)=c(i,7)+a(i,7)
20    continue
       return
       end
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
       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
       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
       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
       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
       subroutine mxv8(a,b,c)
       implicit real*8(a-h,o-z)
       dimension a(14,14),b(14),c(14)
       do 10 i=1,14
       c(i)=0.d0
       do 10 j=1,7
       c(i)=c(i)+a(i,j+7)*b(j+7)
10     continue
       do 20 i=1,14
        c(i)=c(i)+a(i,7)*b(7)
20     continue
       return
       end
c
       subroutine mxtv8(a,b,c)
       implicit real*8(a-h,o-z)
       dimension a(14,14),b(14),c(14)
       do 10 i=1,14
       c(i)=0.d0
       do 10 j=1,7
       c(i)=c(i)+a(j+7,i)*b(j+7)
10     continue
       do 20 i=1,14
        c(i)=c(i)+a(7,i)*b(7)
20     continue
       return
       end
c
       subroutine mtv14(a,b,c)
       implicit real*8(a-h,o-z)
       dimension a(7,14),b(14),c(14)
       do 10 i=1,14
       c(i)=0.d0
       do 10 j=1,7
       c(i)=c(i)+a(j,i)*b(j+7)
10     continue
        c(7)=c(7)+b(7)
       return
       end
c
        subroutine maxa(nofe,n8,imx)
        dimension imx(n8+1)
c
        imx(1)=1
        do 10 j=2,15
        imx(j)=imx(j-1)+j-1
10      continue
        do 20 nj=1,nofe-1
        do 20 j=9,15
        imx(7*nj+j)=imx(7*nj+j-1)+j-1
20      continue
        imx(n8)=imx(n8-1)+14
        imx(n8+1)=imx(n8)+n8
        return
        end
c
c 各要素ごとの接線剛性行列を、スカイライン法の上三角行列と
c 下三角行列にわけながら全要素について重ね合わせる
        subroutine skyln(nofe,i,s,so,stuh,stlh,imx)
        implicit real*8(a-h,o-z)
        dimension s(14,14),so(7,7),           
     &  stuh(84*nofe+36),stlh(84*nofe+36),imx(7*nofe+9)
c
        do 10 j=1,7
        do 10 k=j,7
        ij=7*i+j
        ik=7*i+k
        im=imx(ik)+ik-ij
        stuh(im)=s(j,k)+so(j,k)
        stlh(im)=s(k,j)+so(k,j)
        so(j,k)=s(j+7,k+7)
        so(k,j)=s(k+7,j+7)
10      continue
        do 15 j=1,7
        do 15 k=1,7
        ij=7*i+j
        ik=7*i+k
        im7=imx(ik+7)+ik+7-ij
        stuh(im7)=s(j,k+7)
        stlh(im7)=s(k+7,j)
15      continue
        if(i.eq.nofe-1) then
        do 20 j=1,7
        do 20 k=j,7
        ij=7*i+j
        ik=7*i+k
        i7=imx(ik+7)+ik-ij
        stuh(i7)=s(j+7,k+7)
        stlh(i7)=s(k+7,j+7)
        so(j,k)=0.d0
        so(k,j)=0.d0
20      continue
        endif 
        return
        end
c
c 境界条件を入れる
        subroutine bncn(nofe,bc,stuh,stlh,imx)
        implicit real*8(a-h,o-z)
        dimension bc(7*nofe+8),                
     &  stuh(84*nofe+36),stlh(84*nofe+36),imx(7*nofe+9)
c
        do 10 i=0,nofe-1
        do 10 j=1,7
        do 10 k=j,14
        ij=7*i+j
        ik=7*i+k
        im=imx(ik)+ik-ij
        stuh(im)=bc(ij)*bc(ik)*stuh(im)
        stlh(im)=bc(ij)*bc(ik)*stlh(im)
10      continue
        do 15 j=7*nofe+1,7*nofe+7
        do 15 k=j,7*nofe+7
        im=imx(k)+k-j
        stuh(im)=bc(j)*bc(k)*stuh(im)
        stlh(im)=bc(j)*bc(k)*stlh(im)
15      continue
        n8=7*nofe+8
        do 20 i=1,n8
        if(bc(i).gt.1.d-6) goto 20
        stuh(imx(i))=1.d0
        stlh(imx(i))=0.d0
        im=imx(n8)+n8-i
        stuh(im)=bc(i)*stuh(im)
        stlh(im)=bc(i)*stlh(im)
20      continue
        return
        end
c
c
c 接線剛性方程式を解いて変位増分を求めるとともに、接線剛性行列の
c 負の固有値数を求める。
c 岩熊先生が、東京大学の大計のライブラリーかどこかの「ガウスの消去法」の
c サブルーチンを2次元骨組をスカイライン法で解くために書き換えたものを
c 拝借した。
c
       subroutine solve(stuh,stlh,dp,imx,msl,nnegtm,anegtm,
     &                   det,idt,jdt,n01)
c
c 方程式は掃き出し法によりここで解かれる
c
      implicit real*8 (a-h,o-z)
      dimension stuh(12*n01-60),stlh(12*n01-60)
      dimension imx(n01+1),msl(n01),nnegtm(n01),anegtm(n01)
c
c  行列式 = det D idt
c jdt : 負の固有値の数
c
        do 100 i=1,n01
        msl(i)=1+i+imx(i)-imx(i+1)
100     continue

      det=stuh(1)
      idt=0
      jdt=0
      if( det.gt.0.d0 ) go to 80
      jdt=1
      nnegtm(jdt)=1
      anegtm(jdt)=stuh(1)
   80 continue
c
c 前進消去
c
      do 10 j=2,n01
      ist=msl(j)+1
      ied=j-1
      if( ist.gt.ied ) go to 20
      do 30 i=ist,ied
      kst=msl(j)
      if( msl(i).gt.kst ) kst=msl(i)
      ked=i-1
      ij=imx(j)+j-i
      do 30 k=kst,ked
      kj=imx(j)+j-k
      ki=imx(i)+i-k
      stuh(ij)=stuh(ij)-stuh(kj)*stlh(ki)
      stlh(ij)=stlh(ij)-stlh(kj)*stuh(ki)
   30 continue
   20 continue
      do 40 i=msl(j),ied
      ij=imx(j)+j-i
      stuh(ij)=stuh(ij)/stuh(imx(i))
      stlh(ij)=stlh(ij)/stuh(imx(i))
   40 continue
      do 50 k=msl(j),ied
      kj=imx(j)+j-k
      stuh(imx(j))=stuh(imx(j))-stuh(kj)*stlh(kj)*stuh(imx(k))
      stlh(imx(j))=stlh(imx(j))-stlh(kj)*stlh(imx(k))
   50 continue
c
      if( j.eq.n01 ) go to 10
      if( stuh(imx(j)).gt.0.d0 ) go to 90
      jdt=jdt+1
      nnegtm(jdt)=j
      anegtm(jdt)=stuh(imx(j))
   90 continue
      det=det*stuh(imx(j))
      jjjdt=idint(dlog10(dabs(det)))
      idt=idt+jjjdt
      det=det/10.d0**(jjjdt)
c
   10 continue
c
c 後進消去
c
      do 60 j=1,n01
      stlh(imx(j))=stlh(imx(j))/stuh(imx(j))
   60 continue
      do 70 j=2,n01
      i=n01+2-j
      ked=i-1
      do 70 k=msl(i),ked
      ii=imx(i)+i-k
      stlh(imx(k))=stlh(imx(k))-stuh(ii)*stlh(imx(i))
   70 continue
      dp=stlh(imx(n01))
c
      return
      end
c
c 以下のcl..というサブルーチン群は配列の中に 0 を書き込んで
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
        subroutine cliv(n,l)
        implicit real*8(a-h,o-z)
         dimension l(n)
          do 10 i=1,n
             l(i)=0
10        continue
         return
        end
c
        subroutine cl4rvs(n1,a1,n2,a2,n3,a3,n4,a4)
        implicit real*8(a-h,o-z)
         dimension a1(n1),a2(n2),a3(n3),a4(n4)
          do 10 i=1,n1
             a1(i)=0.d0
10        continue
          do 20 i=1,n2
             a2(i)=0.d0
20        continue
          do 30 i=1,n3
             a3(i)=0.d0
30        continue
          do 40 i=1,n4
             a4(i)=0.d0
40        continue
         return
        end
c
        subroutine cl4ivs(n1,l1,n2,l2,n3,l3,n4,l4)
        implicit real*8(a-h,o-z)
         dimension l1(n1),l2(n2),l3(n3),l4(n4)
          do 10 i=1,n1
             l1(i)=0
10        continue
          do 20 i=1,n2
             l2(i)=0
20        continue
          do 30 i=1,n3
             l3(i)=0
30        continue
          do 40 i=1,n4
             l4(i)=0
40        continue
         return
        end
c
        subroutine clrmx(n,a)
        implicit real*8(a-h,o-z)
         dimension a(n,n)
          do 10 i=1,n
          do 10 j=1,n
             a(i,j)=0.d0
10        continue
         return
        end
c
        subroutine cl4rms(n1,a1,n2,a2,n3,a3,n4,a4)
        implicit real*8(a-h,o-z)
         dimension a1(n1,n1),a2(n2,n2),a3(n3,n3),a4(n4,n4)
          do 10 i=1,n1
          do 10 j=1,n1
             a1(i,j)=0.d0
10        continue
          do 20 i=1,n2
          do 20 j=1,n2
             a2(i,j)=0.d0
20        continue
          do 30 i=1,n3
          do 30 j=1,n3
             a3(i,j)=0.d0
30        continue
          do 40 i=1,n4
          do 40 j=1,n4
             a4(i,j)=0.d0
40        continue
         return
        end
c
        subroutine cl5rms(deti,r7x14,dr7x14,dss,dt0)
        implicit real*8(a-h,o-z)
       dimension deti(14,14,14),r7x14(7,14),dr7x14(14,7,14),
     &           dss(14,14,14),dt0(7,14,14)
c
          do 10 i=1,14
          do 10 j=1,14
          do 10 k=1,14
             deti(i,j,k)=0.d0 
             dss(i,j,k)=0.d0
10        continue
c
          do 20 i=1,7
          do 20 j=1,14
             r7x14(i,j)=0.d0
          do 20 k=1,14
             dr7x14(j,i,k)=0.d0
             dt0(i,j,k)=0.d0
20        continue
         return
        end
c
c
c
c 増分行列の幾つかがここで作られる
c 尚、このサブルーチン内の増分式は、
c 小林裕さんが Mathematica で計算してくれました。
c
       subroutine farejo(d,dl,t,cf1,cf2,dt0,r7x14,dr7x14,deti)
        IMPLICIT REAL*8 (A-H,O-Z)
       dimension d(14),dl(14),t(14,14),cf1(7,7),cf2(7,7),
     &    dt0(7,14,14),r7x14(7,14),dr7x14(14,7,14),deti(14,14,14),
     &  a(6,3),da(14,6,3),dcf1(14,7,7),dcf2(14,7,7)
       dimension tdb1(14,7,7),tdb2(14,7,7),dtb1(14,7,7),
     &           dtb2(14,7,7),dt1(14,14,14),t1(7,7)
       do 1007 i=1,6
       do 1007 j=1,3
       a(i,j)=0.d0
       do 1007 k=1,14
       da(k,i,j)=0.d0
1007    continue
       do 1013 k=1,14
       do 1013 i=1,14
       do 1013 j=1,14 
       deti(k,i,j)=0.d0
1013   continue
       do 1018 i=1,7
       do 1018 j=1,14
       r7x14(i,j)=0.d0
       do 1018 k=1,14
       dr7x14(k,i,j)=0.d0
1018    continue
       do 1024 k=1,14
       do 1024 i=1,7
       do 1024 j=1,7
       dcf1(k,i,j)=0.d0
       dcf2(k,i,j)=0.d0
1024   continue

       do 460 k=1,7
       do 460 i=1,14
       do 460 j=1,14
       dt1(k,i,j)=dt0(k,i,j)
       dt1(7+k,i,j)=0.d0
460    continue
       do 470 i=1,7
       do 470 j=1,7
        t1(i,j)=t(i,j)
470    continue
        u1=d(1)
        v1=d(2)
        w1=d(3)
        a1=d(4)
        g1=d(5)
        p1=d(6)
        q1=d(7)
        u2=d(8)
        v2=d(9)
        w2=d(10)
        a2=d(11)
        g2=d(12)
        p2=d(13)
        q2=d(14)
c
       ca1=cos(a1)
       cg1=cos(g1)
       cp1=cos(p1)
       ca2=cos(a2)
       cg2=cos(g2)
       cp2=cos(p2)
       sa1=sin(a1)
       sg1=sin(g1)
       sp1=sin(p1)
       sa2=sin(a2)
       sg2=sin(g2)
       sp2=sin(p2)
       ta1=tan(a1)
       ta2=tan(a2)
       tg1=tan(g1)
       tg2=tan(g2)
       tp1=tan(p1)
       tp2=tan(p2)
c
       gmm=(g1+g2)/2.d0
       all=(a1+a2)/2.d0
       phh=(p2-p1)/2.d0
       allp=(a2-a1)/2.d0
c
       cgg=cos(gmm)
       caa=cos(all)
       cph=cos(phh)
       callp=cos(allp)
       sgg=sin(gmm)
       saa=sin(all)
       sph=sin(phh)
       sallp=sin(allp)
c
       ru=dl(8)
       rv=dl(9)
       rw=dl(10)
       rx=dl(11)
       ry=dl(12)
       rz=dl(13)
c
       dcf1(4,1,4)=ca1*sg1
       dcf1(5,1,4)=sa1*cg1
       dcf1(4,1,5)=sa1*cg1
       dcf1(5,1,5)=ca1*sg1
       dcf1(4,2,4)=sa1
       dcf1(4,3,4)=ca1*cg1
       dcf1(5,3,4)=-sa1*sg1
       dcf1(4,3,5)=-sa1*sg1
       dcf1(5,3,5)=ca1*cg1
       dcf1(4,4,4)=-phh*sgg*caa/2.d0
       dcf1(11,4,4)=dcf1(4,4,4)
       dcf1(5,4,4)=(-sgg-phh*cgg*saa)/2.d0
       dcf1(12,4,4)=dcf1(5,4,4)
       dcf1(6,4,4)=sgg*saa/2.d0
       dcf1(13,4,4)=-dcf1(6,4,4)
       dcf1(4,4,5)=(-sgg-phh*cgg*saa)/2.d0
       dcf1(5,4,5)=(allp*cgg-phh*sgg*caa)/2.d0
       dcf1(6,4,5)=-cgg*caa/2.d0
       dcf1(11,4,5)=(sgg-phh*cgg*saa)/2.d0
       dcf1(12,4,5)=dcf1(5,4,5) 
       dcf1(13,4,5)=-dcf1(6,4,5)
       dcf1(4,4,6)=sgg*saa/2.d0
       dcf1(5,4,6)=-cgg*caa/2.d0
       dcf1(11,4,6)=dcf1(4,4,6)
       dcf1(12,4,6)=dcf1(5,4,6)
       dcf2(4,4,4)=-phh*sgg*caa/2.d0
       dcf2(5,4,4)=(sgg-phh*cgg*saa)/2.d0
       dcf2(6,4,4)=sgg*saa/2.d0
       dcf2(11,4,4)=dcf2(4,4,4)
       dcf2(12,4,4)=dcf2(5,4,4)
       dcf2(13,4,4)=-dcf2(6,4,4) 
       dcf1(4,5,4)=-saa*phh/2.d0
       dcf1(11,5,4)=dcf1(4,5,4)
       dcf1(4,5,6)=-caa/2.d0
       dcf1(11,5,6)=dcf1(4,5,6)
       dcf1(6,5,4)=-caa/2.d0
       dcf1(13,5,4)=-dcf1(6,5,4)
       dcf1(4,6,4)=-phh*cgg*caa/2.d0
       dcf1(5,6,4)=(-cgg+phh*sgg*saa)/2.d0
       dcf1(6,6,4)=cgg*saa/2.d0
       dcf1(11,6,4)=dcf1(4,6,4)
       dcf1(12,6,4)=dcf1(5,6,4)
       dcf1(13,6,4)=-dcf1(6,6,4)
       dcf1(4,6,5)=(-cgg+phh*sgg*saa)/2.d0
       dcf1(5,6,5)=(-allp*sgg-phh*cgg*caa)/2.d0
       dcf1(6,6,5)=sgg*caa/2.d0
       dcf1(11,6,5)=(cgg+phh*sgg*saa)/2.d0
       dcf1(12,6,5)=dcf1(5,6,5)
       dcf1(13,6,5)=-dcf1(6,6,5)
       dcf1(4,6,6)=cgg*saa/2.d0
       dcf1(5,6,6)=sgg*caa/2.d0
       dcf1(11,6,6)=dcf1(4,6,6)
       dcf1(12,6,6)=dcf1(5,6,6)
       dcf2(4,6,4)=-phh*cgg*caa/2.d0
       dcf2(5,6,4)=(cgg+phh*sgg*saa)/2.d0
       dcf2(6,6,4)=cgg*saa/2.d0
       dcf2(11,6,4)=dcf2(4,6,4)
       dcf2(12,6,4)=dcf2(5,6,4)
       dcf2(13,6,4)=-dcf2(6,6,4)
       dcf2(4,6,5)=dcf1(4,6,5)
       dcf2(5,6,5)=dcf1(5,6,5)
       dcf2(6,6,5)=dcf1(6,6,5)
       dcf2(11,6,5)=dcf1(11,6,5)
       dcf2(12,6,5)=dcf1(12,6,5)
       dcf2(13,6,5)=dcf1(13,6,5) 
       dcf2(4,4,6)=-sgg*saa/2.d0
       dcf2(5,4,6)=cgg*caa/2.d0 
       dcf2(11,4,6)=dcf2(4,4,6)
       dcf2(12,4,6)=dcf2(5,4,6)
        dcf2(4,4,5)=dcf1(4,4,5)
        dcf2(5,4,5)=dcf1(5,4,5)
        dcf2(6,4,5)=dcf1(6,4,5)
        dcf2(11,4,5)=dcf1(11,4,5)
        dcf2(12,4,5)=dcf1(12,4,5)
        dcf2(13,4,5)=dcf1(13,4,5)
        dcf2(4,5,4)=dcf1(4,5,4)
        dcf2(6,5,4)=dcf1(6,5,4)
        dcf2(11,5,4)=dcf1(11,5,4)
        dcf2(13,5,4)=dcf1(13,5,4)
        dcf2(4,5,6)=caa/2.d0
        dcf2(11,5,6)=dcf2(4,5,6)
        dcf2(4,6,6)=-cgg*saa/2.d0
        dcf2(11,6,6)=-cgg*saa/2.d0
        dcf2(5,6,6)=-sgg*caa/2.d0
       dcf2(13,5,5)=-dcf2(6,5,5) 
        dcf2(12,6,6)=-sgg*caa/2.d0
c
       deti(5,4,4)=sg1
       deti(5,5,4)=-cg1*ta1
       deti(4,5,4)=-sg1/ca1/ca1
       deti(5,6,4)=cg1/ca1
       deti(4,6,4)=sg1*sa1/ca1/ca1
       deti(5,4,6)=cg1
       deti(5,5,6)=sg1*ta1
       deti(4,5,6)=-cg1/ca1/ca1
       deti(5,6,6)=-sg1/ca1
       deti(4,6,6)=cg1*sa1/ca1/ca1
       deti(12,11,11)=sg2
       deti(12,12,11)=-cg2*ta2
       deti(11,12,11)=-sg2/ca2/ca2
       deti(12,13,11)=cg2/ca2
       deti(11,13,11)=sg2*sa2/ca2/ca2
       deti(12,11,13)=cg2
       deti(12,12,13)=sg2*ta2
       deti(11,12,13)=-cg2/ca2/ca2
       deti(12,13,13)=-sg2/ca2
       deti(11,13,13)=cg2*sa2/ca2/ca2
c
c
       da(1,1,1)=sg1*t(2,1)
       da(2,1,1)=ta1*t(2,1)
       da(3,1,1)=cg1*t(2,1)
       da(4,1,1)=t(2,1)*(ru*sg1*ta1-rv+rw*cg1*ta1)
       da(5,1,1)=(-ru*cg1+rw*sg1)*t(2,1)
       da(6,1,1)=(-ru*sg1-rv*ta1-rw*cg1)*t(2,2)
       da(8,1,1)=t(2,1)*(-sg1)
       da(9,1,1)=t(2,1)*(-ta1)
       da(10,1,1)=t(2,1)*(-cg1)
c
       da(1,1,2)=-t(3,1)
       da(3,1,2)=t(1,1)
       da(4,1,2)=-t(2,1)*ru*cg1+t(2,1)*rw*sg1-cp1*sa1
       da(5,1,2)=-t(1,1)*ru-t(3,1)*rw+ca1*sa1*sp1
       da(6,1,2)=-t(1,2)*rw+t(3,2)*ru
       da(8,1,2)=t(3,1)
       da(10,1,2)=-t(1,1)
c
       da(1,1,3)=-t(1,2)
       da(2,1,3)=-t(2,2)
       da(3,1,3)=-t(3,2)
       da(4,1,3)=-t(2,2)*(sg1*ru+ta1*rv+cg1*rw)+t(1,2)*sa1*sg1
     &           -t(2,2)*ca1+t(3,2)*sa1*cg1
       da(5,1,3)=t(3,2)*ru-t(1,2)*rw-t(1,2)*ca1*cg1+t(3,2)*ca1*sg1
       da(6,1,3)=-t(1,1)*ru-t(2,1)*rv-t(3,1)*rw
c
c       da(4,1,3)=-ru*t(2,1)*sg1+t(1,1)*sa1*sg1-rv*t(2,2)*ta1
c     &           -t(2,2)*ca1-rw*t(2,2)*cg1+t(3,2)*sa1*cg1
c       da(5,1,3)=ru*t(3,1)-t(1,1)*ca1*cg1-rw*t(1,2)+t(3,2)*ca1*sg1
c       da(6,1,3)=ru*t(1,2)-rv*t(2,1)-rw*t(3,1)
c
       da(8,1,3)=t(1,2)
       da(9,1,3)=t(2,2)
       da(10,1,3)=t(3,2)
c
       da(1,2,1)=t(2,2)*sg1
       da(2,2,1)=t(2,2)*ta1
       da(3,2,1)=t(2,2)*cg1
       da(4,2,1)=t(2,2)*(ru*sg1*ta1-rv+rw*ta1*cg1)
       da(5,2,1)=t(2,2)*(-ru*cg1+rw*sg1)
       da(6,2,1)=t(2,1)*(ru*sg1+rv*ta1+rw*cg1)
       da(8,2,1)=-t(2,2)*sg1
       da(9,2,1)=-t(2,2)*ta1
       da(10,2,1)=-t(2,2)*cg1
c
       da(1,2,2)=-t(3,2)
       da(3,2,2)=t(1,2)
       da(4,2,2)=t(2,2)*(-ru*cg1+rw*sg1)+t(3,2)*sa1*sg1-t(1,2)*sa1*cg1
       da(5,2,2)=-t(1,2)*ru-t(3,2)*rw-t(3,2)*ca1*cg1-t(1,2)*ca1*sg1
       da(6,2,2)=-t(3,1)*ru+t(1,1)*rw
c
c       da(4,2,2)=sa1*sp1-t(2,2)*cg1+t(1,3)*cp1
c       da(5,2,2)=-ca1*sa1*cp1-t(2,2)*cg1+t(1,3)*cp1
c       da(6,2,2)=-t(3,1)*ru-t(1,1)*rw
c
       da(8,2,2)=t(3,2)
       da(10,2,2)=-t(1,2)
c
       do 50 i=1,3
       da(i,2,3)=t(i,1)
       da(7+i,2,3)=-t(i,1)
50     continue
       da(4,2,3)=t(2,1)*ru*sg1+t(2,1)*rv*ta1+t(2,1)*rw*cg1
     &           -t(1,1)*sa1*sg1+t(2,1)*ca1-t(3,1)*sa1*cg1
       da(5,2,3)=ca1*cp1+t(1,1)*rw-t(3,1)*ru
       da(6,2,3)=-t(1,2)*ru-t(2,2)*rv-t(3,2)*rw
c
       da(1,3,1)=sa1*sg1
       da(2,3,1)=-ca1
       da(3,3,1)=sa1*cg1
       da(4,3,1)=-ca1*sg1*ru-sa1*rv-ca1*cg1*rw-1
c
c       da(4,3,1)=sa1**2*sg1**2-ca1**2+sa1**2*cg1**2-ru*sg1*ca1-rv*sa1
c     &           -rw*cg1*ca1
c
       da(5,3,1)=-ru*t(2,3)*cg1+rw*sg1*t(2,3)
       da(8,3,1)=-t(2,3)*sg1
       da(9,3,1)=ca1
       da(10,3,1)=-t(2,3)*cg1
c
       da(1,3,2)=-ca1*cg1
       da(3,3,2)=t(1,3)
       da(4,3,2)=-sa1*(ru*cg1-rw*sg1)
       da(5,3,2)=-ca1**2-ca1*(ru*sg1+rw*cg1)
       da(8,3,2)=ca1*cg1
       da(10,3,2)=-t(1,3)
c
       da(4,4,1)=-cgg*sg1*t(2,1)+phh*sgg*saa*sg1*t(2,1)
     &  -phh*caa*ta1*t(2,1)-ry*t(2,1)/ca1/ca1+sgg*cg1*t(2,1)
     &  +phh*cgg*saa*cg1*t(2,1)+(rx*sg1+ry*ta1+rz*cg1)*ta1*t(2,1)
       da(11,4,1)=cgg*sg1*t(2,1)+phh*sgg*saa*sg1*t(2,1)
     &  -phh*caa*ta1*t(2,1)-sgg*cg1*t(2,1)+phh*cgg*saa*cg1*t(2,1)
       da(5,4,1)=-allp*sgg*sg1*t(2,1)-phh*cgg*caa*sg1*t(2,1)
     &  -rx*cg1*t(2,1)+ta1*t(2,1)-allp*cgg*cg1*t(2,1)
     &  +phh*sgg*caa*cg1*t(2,1)+rz*sg1*t(2,1)
       da(12,4,1)=-allp*sgg*sg1*t(2,1)-phh*cgg*caa*sg1*t(2,1)
     &  -ta1*t(2,1)-allp*cgg*cg1*t(2,1)+phh*sgg*caa*cg1*t(2,1)
       da(6,4,1)=sgg*caa*sg1*t(2,1)+saa*ta1*t(2,1)+cgg*caa*cg1*t(2,1)
     &  +(-rx*sg1-ry*ta1-rz*cg1)*t(2,2)
       da(13,4,1)=-sgg*caa*sg1*t(2,1)-saa*ta1*t(2,1)-cgg*caa*cg1*t(2,1)
c
       da(4,4,2)=-t(2,1)*rx*cg1+t(2,1)*rz*sg1+cgg*t(3,1)
     &  -phh*sgg*saa*t(3,1)+phh*cgg*saa*t(1,1)+sgg*t(1,1)
       da(11,4,2)=-cgg*t(3,1)-phh*sgg*saa*t(3,1)-sgg*t(1,1)
     &  +phh*cgg*saa*t(1,1)
       da(5,4,2)=-t(1,1)*rx-t(3,1)*rz+allp*sgg*t(3,1)
     &  +phh*cgg*caa*t(3,1)-allp*cgg*t(1,1)+phh*sgg*caa*t(1,1)
       da(12,4,2)=allp*sgg*t(3,1)+phh*cgg*caa*t(3,1)-allp*cgg*t(1,1)
     &  +phh*sgg*caa*t(1,1)
       da(6,4,2)=t(3,2)*rx-t(1,2)*rz-sgg*caa*t(3,1)+cgg*caa*t(1,1)
       da(13,4,2)=sgg*caa*t(3,1)-cgg*caa*t(1,1)
c
       da(4,4,3)=-t(2,2)*rx*sg1-t(2,2)*ry*ta1-t(2,2)*rz*cg1+cgg*t(1,2)
     &  -phh*sgg*saa*t(1,2)+phh*caa*t(2,2)-phh*cgg*saa*t(3,2)
     &  -sgg*t(3,2)
       da(11,4,3)=-cgg*t(1,2)-phh*sgg*saa*t(1,2)+phh*caa*t(2,2)
     &  +sgg*t(3,2)-phh*cgg*saa*t(3,2)
       da(5,4,3)=t(3,2)*rx-t(1,2)*rz+allp*sgg*t(1,2)+phh*cgg*caa*t(1,2)
     &  -t(2,2)+allp*cgg*t(3,2)-phh*sgg*caa*t(3,2)
       da(12,4,3)=allp*sgg*t(1,2)+phh*cgg*caa*t(1,2)+t(2,2)
     &  +allp*cgg*t(3,2)-phh*sgg*caa*t(3,2)
       da(6,4,3)=-t(1,1)*rx-t(2,1)*ry-t(3,1)*rz-sgg*caa*t(1,2)
     &  -saa*t(2,2)-cgg*caa*t(3,2)
       da(13,4,3)=sgg*caa*t(1,2)+saa*t(2,2)+cgg*caa*t(3,2)
c
       da(4,5,1)=-cgg*sg1*t(2,2)+phh*sgg*saa*sg1*t(2,2)
     &  -phh*caa*ta1*t(2,2)-ry*t(2,2)/ca1/ca1-sgg*cg1*t(2,2)
     &  -phh*cgg*saa*cg1*t(2,2)+ta1*t(2,2)*(rx*sg1+ry*ta1+rz*cg1)
       da(11,5,1)=cgg*sg1*t(2,2)+phh*sgg*saa*sg1*t(2,2)
     &  -phh*caa*ta1*t(2,2)-sgg*cg1*t(2,2)+phh*cgg*saa*cg1*t(2,2)
       da(5,5,1)=-allp*sgg*sg1*t(2,2)-phh*cgg*caa*sg1*t(2,2)
     &  -rx*t(2,2)*cg1-allp*cgg*cg1*t(2,2)+ta1*t(2,2)
     &  +phh*sgg*caa*cg1*t(2,2)+rz*t(2,2)*sg1
       da(12,5,1)=-allp*sgg*sg1*t(2,2)-phh*cgg*caa*sg1*t(2,2)
     &  -ta1*t(2,2)-allp*cgg*cg1*t(2,2)+phh*sgg*caa*cg1*t(2,2)
       da(6,5,1)=sgg*caa*sg1*t(2,2)+saa*ta1*t(2,2)+cgg*caa*cg1*t(2,2)
     &  +(rx*sg1+ry*ta1+rz*cg1)*t(2,1)
       da(13,5,1)=-sgg*caa*sg1*t(2,2)-saa*ta1*t(2,2)-cgg*caa*cg1*t(2,2)
c
       da(4,5,2)=-cg1*rx*t(2,2)+t(2,2)*rz*sg1+cgg*t(3,2)
     &  -phh*sgg*saa*t(3,2)+sgg*t(1,2)+phh*cgg*saa*t(1,2)
       da(11,5,2)=-cgg*t(3,2)-phh*sgg*saa*t(3,2)-sgg*t(1,2)
     &  +phh*cgg*saa*t(1,2)
       da(5,5,2)=-t(1,2)*rx-t(3,2)*rz+allp*sgg*t(3,2)
     &  +phh*cgg*caa*t(3,2)-allp*cgg*t(1,2)+phh*sgg*caa*t(1,2)
       da(12,5,2)=+allp*sgg*t(3,2)+phh*cgg*caa*t(3,2)-allp*cgg*t(1,2)
     &  +phh*sgg*caa*t(1,2)
       da(6,5,2)=-t(3,1)*rx+t(1,1)*rz-sgg*caa*t(3,2)+cgg*caa*t(1,2)
       da(13,5,2)=sgg*caa*t(3,2)-cgg*caa*t(1,2)
c
       da(4,5,3)=t(2,1)*rx*sg1+t(2,1)*ry*ta1+t(2,1)*rz*cg1-cgg*t(1,1)
     &  +phh*sgg*saa*t(1,1)-phh*caa*t(2,1)+sgg*t(3,1)+phh*cgg*saa*t(3,1)
       da(11,5,3)=cgg*t(1,1)+phh*sgg*saa*t(1,1)-phh*caa*t(2,1)
     &  -sgg*t(3,1)+phh*cgg*saa*t(3,1)
       da(5,5,3)=-t(3,1)*rx+t(1,1)*rz-allp*sgg*t(1,1)
     &  -phh*cgg*caa*t(1,1)+t(2,1)-allp*cgg*t(3,1)+phh*sgg*caa*t(3,1)
       da(12,5,3)=-allp*sgg*t(1,1)-phh*cgg*caa*t(1,1)-t(2,1)
     &  -allp*cgg*t(3,1)+phh*sgg*caa*t(3,1)
       da(6,5,3)=-t(1,2)*rx-t(2,2)*ry-t(3,2)*rz+sgg*caa*t(1,1)
     &  +saa*t(2,1)+cgg*caa*t(3,1)
       da(13,5,3)=-sgg*caa*t(1,1)-saa*t(2,1)-cgg*caa*t(3,1)
c
c       da(4,6,1)=cgg*sg1*t(2,3)+phh*sgg*saa*sg1*t(2,3)
c     &  +phh*caa*ca1-ry*sa1+sgg*cg1*t(2,3)
c     &  +phh*cgg*saa*cg1*t(2,3)+ca1*(-rx*sg1-rz*cg1)
c       da(11,6,1)=-cgg*sg1*t(2,3)+phh*sgg*saa*sg1*t(2,3)-sgg*cg1*t(2,3)
c     &  +phh*caa*ca1+phh*cgg*saa*cg1*t(2,3)
c       da(5,6,1)=allp*sgg*sg1*t(2,3)-phh*caa*cgg*sg1*t(2,3)
c     &  -rx*t(2,3)*cg1-ca1-allp*cgg*cg1*t(2,3)
c     &  +phh*sgg*caa*cg1*t(2,3)
c       da(12,6,1)=allp*sgg*sg1*t(2,3)-phh*cgg*caa*sg1*t(2,3)
c     &  +ca1-allp*cgg*cg1*t(2,3)+phh*sgg*caa*cg1*t(2,3)
c       da(6,6,1)=sgg*caa*sg1*t(2,3)-saa*ca1
c     &  +cgg*caa*cg1*t(2,3)
c       da(13,6,1)=-sgg*caa*sg1*t(2,3)+saa*ca1
c     &  -cgg*caa*cg1*t(2,3)
c
       da(4,6,1)=-ca1*sg1*rx-sa1*ry-ca1*cg1*rw-cgg*sa1*sg1
     &  +phh*sgg*saa*sa1*sg1+phh*caa*ca1+sgg*sa1*cg1
     &  +phh*cgg*saa*sa1*cg1
       da(11,6,1)=cgg*sa1*sg1+phh*sgg*saa*sa1*sg1+phh*caa*ca1
     &  -sgg*sa1*cg1+phh*cgg*saa*sa1*cg1
       da(5,6,1)=-sa1*cg1*rx+sa1*sg1*rw-allp*sgg*sa1*sg1
     &  -phh*cgg*caa*sa1*sg1-ca1-allp*cgg*sa1*cg1+phh*saa*caa*sa1*cg1
       da(12,6,1)=-allp*sgg*sa1*sg1-phh*cgg*caa*sa1*sg1
     &  +ca1-allp*cgg*sa1*cg1+phh*saa*caa*sa1*cg1
       da(6,6,1)=sgg*caa*sa1*sg1-saa*ca1+cgg*caa*sa1*cg1
       da(13,6,1)=-da(6,6,1)
c
       da(4,6,2)=-rx*sa1*cg1+rz*sa1*sg1+cgg*cg1*ca1
     &  -phh*sgg*saa*cg1*ca1+sgg*t(1,3)+phh*cgg*saa*t(1,3)
       da(11,6,2)=-cgg*cg1*ca1-phh*sgg*saa*cg1*ca1
     &  -sgg*t(1,3)+phh*cgg*saa*t(1,3)
       da(5,6,2)=-rz*ca1*cg1-rx*ca1*sg1
     &  +allp*sgg*cg1*ca1+phh*cgg*caa*cg1*ca1
     &  -allp*cgg*t(1,3)+phh*sgg*caa*t(1,3)
       da(12,6,2)=allp*sgg*cg1*ca1+phh*cgg*caa*cg1*ca1
     &  -allp*cgg*t(1,3)+phh*sgg*caa*t(1,3)
       da(6,6,2)=-sgg*caa*cg1*ca1+cgg*caa*t(1,3)
       da(13,6,2)=sgg*caa*cg1*ca1-cgg*caa*t(1,3)
c
c
        da(4,5,1)=
     &   -((-p1 + p2)*Cos((a1 + a2)/2)*Cos(p1)*Sin(a1))/2 - 
     &   Cos(a1)*Cos(p1)*(-g1 + g2 + (-p1 + p2)*Sin((a1 + a2)/2)) - 
     &   Cos(a1)*Cos(g1)*Cos(p1)*(-((-p1 + p2)*Cos((g1 + g2)/2)*
     &          Sin((a1 + a2)/2))/2 - Sin((g1 + g2)/2)) + 
     &   Cos(g1)*Cos(p1)*Sin(a1)*((-p1 + p2)*Cos((a1 + a2)/2)*
     &      Cos((g1 + g2)/2) + (-a1 + a2)*Sin((g1 + g2)/2)) + 
     &   Cos(p1)*Sin(a1)*Sin(g1)*(-((-a1 + a2)*Cos((g1 + g2)/2)) + 
     &     (-p1 + p2)*Cos((a1 + a2)/2)*Sin((g1 + g2)/2)) - 
     &   Cos(a1)*Cos(p1)*Sin(g1)*(Cos((g1 + g2)/2) - 
     &     (-p1 + p2)*Sin((a1 + a2)/2)*Sin((g1 + g2)/2)/2)
c
         da(4,6,1)=
     &   (-p1 + p2)*Cos(a1)*Cos((a1 + a2)/2)/2 - 
     &      Sin(a1)*(-g1 + g2 + (-p1 + p2)*Sin((a1 + a2)/2)) - 
     &      Cos(g1)*Sin(a1)*(-((-p1 + p2)*Cos((g1 + g2)/2)*
     &             Sin((a1 + a2)/2))/2 - Sin((g1 + g2)/2)) - 
     &      Cos(a1)*Cos(g1)*((-p1 + p2)*Cos((a1 + a2)/2)*
     &          Cos((g1 + g2)/2) + (-a1 + a2)*Sin((g1 + g2)/2)) -
     &      Cos(a1)*Sin(g1)*(-((-a1 + a2)*Cos((g1 + g2)/2)) + 
     &         (-p1 + p2)*Cos((a1 + a2)/2)*Sin((g1 + g2)/2)) - 
     &      Sin(a1)*Sin(g1)*(Cos((g1 + g2)/2) - 
     &         (-p1 + p2)*Sin((a1 + a2)/2)*Sin((g1 + g2)/2)/2) 
c
       da(5,6,1)=
     &  -Cos(a1) - Sin(a1)*Sin(g1)*((-p1 + p2)*Cos((a1 + a2)/2)*
     &           Cos((g1 + g2)/2)/2 + (-a1 + a2)*Sin((g1 + g2)/2)/2) + 
     &      Sin(a1)*Sin(g1)*((-p1 + p2)*Cos((a1 + a2)/2)*
     &          Cos((g1 + g2)/2) + (-a1 + a2)*Sin((g1 + g2)/2)) - 
     &      Cos(g1)*Sin(a1)*((-a1 + a2)*Cos((g1 + g2)/2)/2 - 
     &         (-p1 + p2)*Cos((a1 + a2)/2)*Sin((g1 + g2)/2)/2) - 
     &      Cos(g1)*Sin(a1)*(-((-a1 + a2)*Cos((g1 + g2)/2)) + 
     &         (-p1 + p2)*Cos((a1 + a2)/2)*Sin((g1 + g2)/2))
c
c
       a(1,1)=(-ru*sg1-rv*ta1-rw*cg1)*t(2,1)
       a(1,2)=t(3,1)*ru-t(1,1)*rw
       a(1,3)=t(1,2)*ru+t(2,2)*rv+t(3,2)*rw
       a(2,1)=t(2,2)*(-ru*sg1-rv*ta1-rw*cg1)
       a(2,2)=t(3,2)*ru-t(1,2)*rw
       a(2,3)=-t(1,1)*ru-t(2,1)*rv-t(3,1)*rw
       a(3,1)=(-ru*sg1*sa1+rv*ca1-rw*cg1*sa1)
       a(3,2)=ca1*ru*cg1-t(1,3)*rw
       a(3,3)=0.d0
       a(4,1)=t(2,1)*(-rx*sg1-ry*ta1-rz*cg1)
       a(4,2)=t(3,1)*rx-t(1,1)*rz
       a(4,3)=t(1,2)*rx+t(2,2)*ry+t(3,2)*rz
       a(5,1)=t(2,2)*(-rx*sg1-ry*ta1-rz*cg1)
       a(5,2)=t(3,2)*rx-t(1,2)*rz
       a(5,3)=-t(1,1)*rx-t(2,1)*ry-t(3,1)*rz
       a(6,1)=(-rx*sg1*sa1+ry*ca1-rz*cg1*sa1)
       a(6,2)=ca1*rx*cg1-t(1,3)*rz
       a(6,3)=0.d0
c
       do 100 i=1,3
       do 100 j=1,3
           r7x14(i,j)=-t(j,i)
           r7x14(i,7+j)=t(j,i)
100    continue
       do 110 i=1,3
        r7x14(3+i,4)=a(3+i,1)+t(1,i)*cf1(4,4)+t(2,i)*cf1(5,4)
     & +t(3,i)*cf1(6,4)
        r7x14(i,4)=
     &      a(i,1)+t(1,i)*cf1(1,4)+t(2,i)*cf1(2,4)+t(3,i)*cf1(3,4)
        r7x14(i,5)=a(i,2)+t(1,i)*cf1(1,5)+t(3,i)*cf1(3,5)
        r7x14(3+i,5)=a(3+i,2)+t(1,i)*cf1(4,5)-t(2,i)+t(3,i)*cf1(6,5)
        r7x14(i,6)=a(i,3)
        r7x14(3+i,6)=
     &      a(3+i,3)+t(1,i)*cf1(4,6)+t(2,i)*cf1(5,6)+t(3,i)*cf1(6,6)
        r7x14(3+i,11)=t(1,i)*cf2(4,4)+t(2,i)*cf1(5,4)+t(3,i)*cf2(6,4)
        r7x14(3+i,12)=t(1,i)*cf2(4,5)+t(2,i)+t(3,i)*cf1(6,5)
        r7x14(3+i,13)=t(1,i)*cf2(4,6)-t(2,i)*cf1(5,6)-t(3,i)*cf1(6,6)
110    continue
       r7x14(7,14)=1.d0
c
        do 433 k=1,14
        do 433 i=1,7
        do 433 l=1,7
          tdb1(k,i,l)=0.d0
          tdb2(k,i,l)=0.d0
          dtb1(k,i,l)=0.d0
          dtb2(k,i,l)=0.d0
        do 433 j=1,7
       tdb1(k,i,l)=tdb1(k,i,l)+t1(j,i)*dcf1(k,j,l)
       tdb2(k,i,l)=tdb2(k,i,l)+t1(j,i)*dcf2(k,j,l)
       dtb1(k,i,l)=dtb1(k,i,l)+dt1(k,j,i)*cf1(j,l)
       dtb2(k,i,l)=dtb2(k,i,l)+dt1(k,j,i)*cf2(j,l)
433     continue
       do 444 k=1,14
       do 444 i=1,7
       do 444 j=1,7
        dr7x14(k,i,j)=tdb1(k,i,j)+dtb1(k,i,j)
        dr7x14(k,i,7+j)=tdb2(k,i,j)+dtb2(k,i,j)
444     continue
       do 445 k=1,14
       do 445 i=1,6
       do 445 j=1,3
         dr7x14(k,i,j+3)=da(k,i,j)+dr7x14(k,i,j+3)
445    continue
c
       return
       end
c
c 
c 出力した接線行列の固有値解析をせずに、最後のN-R反復で求まる
c (微少な)変位増分を固有ベクトル(座屈モード)の近似として
c 利用する場合のために、近似モードを出力する。
         subroutine modeout(nofe,dd)
         implicit real*8(a-h,o-z)
          dimension dd(7*256+8)
        open(9,file='modeu.d')
         call modeo7(1,nofe,dd)
        close(9)
        open(9,file='modev.d')
         call modeo7(2,nofe,dd)
        close(9)
        open(9,file='modew.d')
         call modeo7(3,nofe,dd)
        close(9)
        open(9,file='modex.d')
         call modeo7(4,nofe,dd)
        close(9)
        open(9,file='modey.d')
         call modeo7(5,nofe,dd)
        close(9)
        open(9,file='modez.d')
         call modeo7(6,nofe,dd)
        close(9)
        open(9,file='modep.d')
         call modeo7(7,nofe,dd)
        close(9)
          return
           end
c
        subroutine modeo7(i7,nofe,dd)
         implicit real*8(a-h,o-z)
          dimension dd(7*256+8)
         do 10 i=1,nofe-1
          dmode=dd(i*7+i7)
          write(9,*) i+1, dmode
10      continue
          return
           end
c
c 変位dに外乱(座屈モードなど)を読み込む
c
         subroutine modein(nofe,d,gmode)
         implicit real*8(a-h,o-z)
         dimension d(14*256)
         open(9,file='modeu.d')
           call modei7(1,nofe,d,gmode)
         close(9)
         open(9,file='modev.d')
           call modei7(2,nofe,d,gmode)
         close(9)
         open(9,file='modew.d')
           call modei7(3,nofe,d,gmode)
         close(9)
         open(9,file='modex.d')
           call modei7(4,nofe,d,gmode)
         close(9)
         open(9,file='modey.d')
           call modei7(5,nofe,d,gmode)
         close(9)
         open(9,file='modez.d')
           call modei7(6,nofe,d,gmode)
         close(9)
         open(9,file='modep.d')
           call modei7(7,nofe,d,gmode)
         close(9)
           return
           end
c
         subroutine modei7(i7,nofe,d,gmode)
         implicit real*8(a-h,o-z)
         dimension d(14*256)
          read(9,*) idummy,dmode
          d(i*14+i7)=d(i*14+i7)+gmode*dmode
         do 10 i=0,nofe-2
          i1=i+1
          read(9,*) idummy,dmode
          d(i*14+7+i7)=d(i*14+7+i7)+gmode*dmode
          d(i1*14+i7)=d(i1*14+i7)+gmode*dmode
10       continue
          read(9,*) idummy,dmode
          n1=nofe-1
          d(n1*14+7+i7)=d(n1*14+7+i7)+gmode*dmode
          return
          end
c
c