c mozi-koodo: EUC
c 集成材梁の有限変位・弾塑性解析プログラム。
c 圧縮側の軸方向応力が降伏応力を越えたら塑性化。
c 引張側の軸方向応力が破断応力を越えたら破壊。
c 金属材料は想定してないので、von Mises の降伏条件は使ってない。
c 断面は長方形のみ。
c 載荷履歴に依存する弾塑性問題では、不釣り合い力をトータルな値
c として扱えないので（後述）、荷重増分を十分に小さく与えること。
c その意味で（接線剛性方程式の定式化はトータルラグラジアンだが）
c 疑似的なアップデートラグラジアンになっている。
c harisen.f から改造。
c 
c 定式化は、
c 『後藤文彦，小林　裕，斉木　功，岩熊哲夫：
c   空間固定三軸回りの回転自由度を用いた空間梁解析
c   応用力学論文集, 1, 1998, pp.319-327.』
c の式(48)で座標変換をオイラー角で表した接線剛性方程式を用いた。
c 尚、反りに関して、ねじり率の自由度は、
c  『後藤文彦，小林裕，岩熊哲夫：
c   オイラー角を用いた簡潔な有限変位解析手法
c   構造工学論文集, 43A, 1997, pp.333-338.
c の定式化のように扱った。
c
c 04/8/26版
c プログラムを主に作った人：後藤文彦
c 増分式を Mathematica で計算してくれた人：小林裕さん
c いくつかの代数計算やスカイライン法のサブルーチンの
c 元プログラムの作者：岩熊先生(http://www.civil.tohoku.ac.jp/~bear/)
c 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c このプログラムは無保証です。                           c
c このプログラムは改造／再配布して構いません。           c
c 但し、私的利用に留まらない場合は、もとのプログラムが   c
c http://www.str.ce.akita-u.ac.jp/~gotou/programoj/      c
c から取ってきたものと分かるように明示して下さい。       c
c 改造／再配布後も上記の条件を保持して下さい。           c
c バグの報告は、                                         c
c http://www.str.ce.akita-u.ac.jp/~gotou/jpg/meeru.jpg   c
c にお願いします。                                       c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         program  harisori
        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),
     &  dld(14*256),vec14(14),sgm(256,200,200),q(7*256+7),
     &  szenbu(256,14,14),dld14(14),ts14(14,14)
c
c 頭にclと付くサブルーチンは、配列類の初期化
         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 clrv(7*256+7,q)
         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 変数の初期化
         ninc=0
         ddi2=0
         ddl2=0
         nlet=0
         p1=0.d0
         nbcl=0
         nacl=0
         iyonda=0
c
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)
c
c 収束精度：収束時の荷重が、N-R反復における一つ前のステップの
c 荷重と有効数字何桁まで合うようにしたいかに応じて tol の小ささを
c 決める
        tol=1.d-10
c 外乱の大きさに応じて収束精度を甘くしたいときとかは、
cc       if(gmode.gt.1.d-10) tol=tol*1.d3
c みたいなのを入れたりするとか。cnvr は座屈点探査の精度。
c
c 初期状態における全ての要素の方向をオイラー角で与える
c （このオイラー角は右手系ではないので注意）
        do 8 i=1,nofe
         read(*,510) alp(i),gmm(i),phi(i)
8       continue
510     format(1p3d17.9) 
c
c  断面定数
c ヤング率：e
          read(*,*) e
c 降伏後の軸（繊維）方向接線弾性係数：et
          read(*,*) et
c 軸圧縮降伏応力：sgmy
          read(*,*) sgmy
c 圧縮：負なので sgmy の符号を負にしておく
          sgmy=-abs(sgmy)
c 軸引張破断応力 sgmu
          read(*,*) sgmu
c せん断弾性係数：g
          read(*,*) g
c 長方形断面の幅 haba
          read(*,*) haba
c 長方形断面の高さ taka
          read(*,*) taka
c 長方形断面の幅の断面要素分割数 nhaba
          read(*,*) nhaba
c 長方形断面の高さの断面要素分割数 ntaka
          read(*,*) ntaka
c 断面積：a
          read(*,*) a
          ea=e*a
c x軸回りの断面二次モーメント:xi
          read(*,*) xi
          eix=e*xi
c y軸回りの断面二次モーメント:yi
          read(*,*) yi
          eiy=e*yi
c ねじり定数：cj
          read(*,*) cj
          gj=g*cj
c そりねじり定数:wi
          read(*,*) wi
          eiw=e*wi
c el:一要素の長さ
        read(*,*) el
c 梁の軸長：ell
        ell=el*nofe
c y軸回りの曲げによるxz面のせん断変形に対するせん断補正係数:xzk
        read(*,*) xzk
c x軸回りの曲げによるyz面のせん断変形に対するせん断補正係数:yzk
        read(*,*) yzk
c 構造工学シリーズ7『構造工学における計算力学の基礎と応用』(土木学会)
c のp.94,95 のφxy,φxz に対応するパラメーター：pxz, pyz
        pxz=eiy/(g*xzk*a*el**2)
        pyz=eix/(g*yzk*a*el**2)
c enko.fでyz面内に横たわるアーチにした場合の中心角：th
        th=-alp(nofe)*2.d0
c π(rad)：pi
        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=(mode*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 軸応力を（gnuplot作図用に）出力したい要素の要素番号：nsgm
          read(*,*) nsgm
c
c 一番最初の一発目の増分ステップの弧長：r
        read(*,*) r
c ２ステップ目以降の計算に使う弧長増分：rn
        read(*,*) rn
c
c
c ############# ここまでがデータ入力 #############
c
c 座屈点探査の精度：「アキレスと亀」（の系：分割のパラドックス
c 「目的地に達するにはその半分の点に達しなければならず、その半分の
c 点に達するには更にその半分の点に達しなければならず……」）
c 方式で座屈点を探査する際に、前の荷重ステップのN-R収束荷重と、
c 今の荷重ステップのN-R収束荷重とが、
c 有効数字何桁まで合うようにしたいかに応じて cnvr の小ささを
c 決める
        cnvr=5.d-8
c 各N-R収束過程における収束精度は tol で与える。
c
c N-R収束荷重pが前のステップの収束荷重p1とある精度内で同じ(p1=p)に
c なった回数：ip1p
        ip1p=0
c
c ###  初期条件 ### 
c
c 変位増分ベクトル：dd(自由度数+1）
c 変位増分の回転角成分は空間固定3軸回りの微少回転角
c 変位ベクトル：d(節点自由度数*要素数）
c トータルな変位は、要素どうしが繋がる２節点のオイラー角を
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
c 各断面要素の軸応力 sgm(要素番号,x座標,y座標)の初期化
         do i=1,nofe
          do l=1,200
          do m=1,200
           sgm(i,l,m)=0.d0
          end do
          end do
         end do
c 弧長の初期値
         ro=r
         dpl=r
c まだ一度も収束値が求まっていないことを表す:nccl（後述）
         nccl=0
c  ### 各荷重ステップの一発目の弧長増分過程 ###
1000   continue        
           nofcs=0
c 計算開始からの増分繰り返し数：ninc
        ninc=ninc+1
c         if(ninc.gt.6000) then
c           write(*,*)'(増分繰り返しが 6000 を越えました)'
c           close(8)
c           stop
c         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)
c
c dldは各増分ステップで（1発目の弧長+N-R過程で）収束した荷重増分
c 要素ごとに前のステップのdld14をdldから読み出す。
           dld14(j)=dld(i*14+j)
c dld14を読み出したらdldはゼロに戻して良い。
           dld(i*14+j)=0.d0
c
120     continue
         do 130 j=1,4
          lg4(j)=lorg(i*4+j)
130      continue
c
c 前のステップのトータルな d14 で要素ごとに座標変換 ts14 を作る
        call trns(ts14,d14(4),d14(5),d14(6))
c （直線梁に限定するので初期形状の座標変換は作らない）
c  ts14 の転置を掛けて居所系の dld14 を求める
         call mxtv(14,ts14,dld14,vec14)
c （直線梁に限定するので初期形状の座標変換はなし）
c
c dd14とd14を使って1要素ずつ接線剛性行列を計算していく
c 求まった1要素ぶんの接線剛性行列はskに入る
c
c ----------------------------------------
c 要素番号：nyou
            nyou=i+1
c 前のステップのdld14を使って剛性行列を作る
c まずは断面要素の降伏判定をして断面定数を断面積分で求める
c 弾性なら、c11=EA,c22=EIy,c33=EIx,c44=EIw,c55=GJ
c
      call danmen(nyou,nofe,nhaba,ntaka,el,haba,taka,
     &            e,et,g,sgmy,sgmu,sgm,vec14,pi,c11,c22,c33,c44)
c
c 求まった断面定数で線形剛性行列（反りを含む1節点7自由度）を
c s(14,14) s に入れる
c
        call smx(el,eiy,pxz,pyz,c11,c22,c33,c44,gj,s)
c
c  要素ごとの剛性行列 s を全要素について szenbu に覚えさせる
         do ii=1,14
         do ij=1,14
            szenbu(nyou,ii,ij)=s(ii,ij)
         end do
         end do
c
c        stmxに行く前にdld14はゼロにする必要がある        
         call clrv(14,dld14)
c
c 接線剛性行列を計算する。
        call stmx(nofe,i,alp,gmm,dd14,d14,sk,b14,lg4,s,dld14)
c
        do 140 j=1,14
          d(i*14+j)=d14(j)
140      continue
c
c stmx から返ってきたdld14は局所系だから全体系に変換する。
c ts14 を掛けて全体系の vec14 を求める
        call mxv(14,ts14,dld14,vec14)
cc vec14 を固定座標まで回転させる　ヤメ（理由は忘れだけど：
cc        call mxv(14,ti14,dld14,vec14)
cc 初期座標方向の局所系で dld14 を定義すべきと判断したのかな）
c
c このステップの一発めの dld14 を dld の中に要素全部に
c ついて覚えさせる．
        do j=1,14
         dld(i*14+j)=vec14(j)
        end do
c
c sk をスカイライン法の1次元ベクトルに書き込む
c stuh に上三角、stlh に下三角が入る
        call skyln(nofe,i,sk,sko,stuh,stlh,imx)
1100   continue
c
cccccccccccccccccccccccccccccccccccccccccccccccc
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 適当な回数に達したら強制的に終了
        nlet=nlet+1
          if(nlet.gt.10) then
           write(*,*) '収束ステップが10を越えたので強制終了'
           stop
          endif
        goto 4001
       endif
        if(arc.lt.tol) goto 3000
c 不釣合力 b の初期化（最終的な不釣合力を出すのは1200の下）
        do 190 i=1,n7
          b(i)=0.d0
cc トータルラグラジアンの場合は、トータルな不釣合力を用いる
cc        b(i)=p*f(i)
190     continue
c
c N-R 収束過程の要素ごとのループのはじまり
c
        do 1200 i=0,nofe-1
        do 200 j=1,14
        dd14(j)=dd(i*7+j)
        d14(j)=d(i*14+j)
c 前の段階でのdld14をdld から読み出してくる
        dld14(j)=dld(i*14+j)
        b14(j)=b(i*7+j)
200     continue
c
        do 210 j=1,4
         lg4(j)=lorg(i*4+j)
210     continue
c
c stmx へ局所系のdldを送るため dld14 を局所系まで回転させる
cc dld14 は初期形状座標で使うので（上述）、ここでdld14を
cc 初期形状座標まで回転させる必要はない。
cc        call mxtv(14,ti14,dld14,vec14)　ヤメ
c dld14 に ts14 の転置を掛けて居所系の dld14 を求める
        call mxtv(14,ts14,dld14,vec14)
c
c szenbu から要素ごとの剛性行列 s を読み出す
c 要素番号：nyou
        nyou=i+1
        do ii=1,14
        do ij=1,14
           s(ii,ij)=szenbu(nyou,ii,ij)
        end do
        end do
c
c 接線剛性行列を１要素ずつ計算してskに入れる
        call stmx(nofe,i,alp,gmm,dd14,d14,sk,b14,lg4,s,dld14)
c
c         print *, 'sk'
c         print '(1p7d11.3)', ( (sk(i,j),j=8,14),i=8,14 )
c         print '(1p7d11.3)',  (b14(j),j=1,7)
c         print '(1p7d11.3)',  (b14(j),j=8,14)
c         stop
c
c
        do 220 j=1,14
        d(i*14+j)=d14(j)
220     continue
c
c stmx から返ってきたdld14は局所系だから全体系に変換する。
c ts14 を掛けて全体系の vec14 を求める
        call mxv(14,ts14,dld14,vec14)
cc vec14 を固定座標まで回転させる　ヤメ（理由は忘れだけど：
cc        call mxv(14,ti14,dld14,vec14)
cc 初期座標方向の局所系で dld14 を定義すべきと判断したのかな）
c
c このステップの一発めの dld14 を dld の中に要素全部に
c ついて覚えさせる．
        do j=1,14
         dld(i*14+j)=vec14(j)
        end do
c
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
ccccccccc 要素ごとのループの終わり cccccccccccccccc
c
c 毎ステップで求まる増分変位Δdで求めた KΔd を累積した
c b=b+KΔd を目標荷重 pf から引いて不釣合力を求める。
c
        do i=1,n7
         q(i)=q(i)+b(i)
         b(i)=p*f(i)-q(i)
        end do
c
        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
c       弾塑性解析で増分荷重を小さくしすぎると座屈点に達する前の
c       何でもないとこでp1=pとなって終了してしまうことがあるので、
c       p1=pが2回ぐらい連続しないと終了しないようにする。
        ip1p=ip1p+1
        print *,'前ステップと収束荷重が',ip1p,'回同じになりました'
         if(ip1p.gt.1) then
           rn=rn*2.d0
           call sgmgnu(nsgm,nhaba,ntaka,sgm)
              if(modeo.eq.1) then
               write(7) nofe+1,stuh,stlh,imx,imx9
               rewind 7
              endif
          close(7)
          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
        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
c
c         print *, 'dd'
c         print '(1p7d11.3)',  (dd(j),j=8,14)
c
        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/EIy=',p*ell/eiy
c      write(*,581) ' p=',p,'   ML^2/8EI=',p*ell**2/8.d0/eix
       write(*,580) ' 修Vla正=',vlp*ell/eiy,'   修Vla負=',vln*ell/eiy
       goto 579
573    continue
        rmd=ell/sqrt(eiy/ea)
        write(*,581)' p=',p,' PL^2/EIy=',p*ell**2/eiy,' λ=',rmd
       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(*,*) ' zakutu'
         write(*,*) ' 座屈したみでだよ................'
           call sgmgnu(nsgm,nhaba,ntaka,sgm)
         stop
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 断面要素の降伏判定をしながら断面定数を断面積分する
      subroutine danmen(nyou,nofe,nhaba,ntaka,el,haba,taka,
     &            e,et,g,sgmy,sgmu,sgm,dld14,pi,c11,c22,c33,c44)
c
      implicit real*8 (a-h,o-z)
      dimension x(200),y(200),depsz1(200,200),depsz2(200,200),
     & dsigz(200,200),eet(200,200),sgm(256,200,200),dld14(14)
c
c まず、この要素の各断面要素の et を求める
      do l=1,nhaba
      do m=1,ntaka
c
c nyou番目の要素の断面要素（座標：l,m）の軸応力：sigma
         sigma=sgm(nyou,l,m)
c
c 集成材の場合は、圧縮のみ弾塑性、引張は弾性域で破断なので
c 降伏判定は圧縮側つまりsigma<0のときのみ行う。
c （引張側は破断判定のみ行う）圧縮：負なのでsgmy<0で与えてある
       if(sigma.lt.sgmy) then
        eet(l,m)=et
         if((nyou.eq.2).and.(l.eq.1).and.(m.eq.1)) then
          write(*,*) ' sosei'
          write(*,*) '要素,断面座標,et,sgm=',nyou,l,m,eet(l,m),sigma
         endif
c 降伏判定された sigma は sgmy に置き換えて、次の dsgm で
c 除荷されたら、eet=e になるようにしておく（除荷は e で戻るように）。
        sgm(nyou,l,m)=sgmy
       else 
         if(sigma.gt.sgmu) then
          write(*,*) ' hadan'
          write(*,*) ' 要素番号：',nyou,' 断面要素座標：',l,m,
     &               'で破断が生じました'
          write(*,*) 'sigma=',sigma
           call sgmgnu(nyou,nhaba,ntaka,sgm)
          stop
         end if
        eet(l,m)=e
      end if
      end do
      end do
c
          an=haba/float(nhaba)
          bn=taka/float(ntaka)
          f=an*bn
c 断面図心を原点とした場合の断面要素の座標値を求める
          nhaba2=nhaba/2
          ntaka2=ntaka/2
        do i=1,nhaba2
         x(nhaba2-i+1)= -(an*i -an/2.d0)
         x(nhaba2+i)=     an*i -an/2.d0
        end do
        do j=1,ntaka2
         y(ntaka2-j+1)=   bn*j -bn/2.d0
         y(ntaka2+j)=   -(bn*j -bn/2.d0)
        end do
c
c 節点変位を以下の変数に書き写す
      u1=  dld14(1) *el
      v1=  dld14(2) *el
      w1=  dld14(3) *el
      thx1=dld14(4)
      thy1=dld14(5)
      thz1=dld14(6)
      q1=  dld14(7) /el
      u2=  dld14(8) *el
      v2=  dld14(9) *el
      w2=  dld14(10)*el
      thx2=dld14(11)
      thy2=dld14(12)
      thz2=dld14(13)
      q2=  dld14(14)/el
c
c 断面定数を断面積分する前に初期化
      c11=0.d0
      c22=0.d0
      c33=0.d0
      c44=0.d0
c     c55=0.d0
c
c
c 要素両端の各断面要素の軸方向歪成分を計算
c
       do l=1,nhaba
       do m=1,ntaka
           xl=x(l)
           ym=y(m)
           call sori(7,haba,taka,xl,ym,pi,w)
c
      depsz1(l,m)=(w2-w1)/el
     &         -(6.d0/(el**2)*(u2-u1)+2.d0/el*(2.d0*thy1+thy2))
     &          *(xl-ym*thz1)
     &         -(6.d0/(el**2)*(v2-v1)+2.d0/el*(2.d0*thx1+thx2))
     &          *(ym+xl*thz1)
     &         -w*(6.d0/(el**2)*(thz2-thz1)+2.d0/el*(2.d0*q1+q2))
     &         +1.d0/2.d0*((thy1**2)+(thx1**2)
     &         +((xl**2)+(ym**2))*(q1**2))
      depsz2(l,m)=(w2-w1)/el
     &         -(6.d0/(el**2)*(u1-u2)-2.d0/el*(thy1+2.d0*thy2))
     &          *(xl-ym*thz2)
     &         -(6.d0/(el**2)*(v1-v2)-2.d0/el*(thx1+2.d0*thx2))
     &          *(ym+xl*thz2)
     &         -w*(6.d0/(el**2)*(thz1-thz2)-2.d0/el*(q1+2.d0*q2))
     &         +1.d0/2.d0*((thy2**2)+(thx2**2)
     &         +((xl**2)+(ym**2))*(q2**2))
c
c 各断面要素（ファイバー要素）の軸方向応力成分を求める
       dsigz(l,m)=eet(l,m)*( depsz1(l,m) + depsz2(l,m) )/2.d0
       sgm(nyou,l,m)=sgm(nyou,l,m)+dsigz(l,m)
c
c
c 断面定数の断面積分
c 初期不整のある問題や座屈後など、断面の降伏が非対称に生じる問題では、
c c12,c13などの断面一次モーメントも計算する必要がある。
c
           c11=c11+eet(l,m)*f
           c22=c22+eet(l,m)*xl**2*f
           c33=c33+eet(l,m)*ym**2*f
           c44=c44+eet(l,m)* w**2 *f
       end do
       end do
c
         haba2=haba/2.d0
         taka2=taka/2.d0
c         call timj(7,haba2,taka2,pi,tj)
c         c55=g*tj * taka*haba**3/3.d0
c
c       write(*,*) 'eet='
c       write(*,*) (eet(i,i),i=1,5)
c       write(*,*) 'a=',c11/eet(1,1)
c       write(*,*) 'xi=',c33/eet(1,1)
c       write(*,*) 'yi=',c22/eet(1,1)
c       write(*,*) 'wi=',c44/eet(1,1)
c       write(*,*) 'ej=', c55/g
c       stop
c
      return
      end
c
c 初期形状が直線梁でないものや大きな回転を伴う問題を解くために
c 断面要素の軸応力を3*3の応力成分の33成分に入れて座標変換するには、
c 以下のような応力の座標変換が必要となる。これを全断面要素について
c やるには、sgm(256,200,200,3,3)のようなサイズの配列が必要になるし、
c 計算時間も増える。その割には、直線梁の座屈点の計算程度には、殆ど
c 影響しないと考えて今回は使っていない。曲がり梁でも座屈点までなら
c 軸応力成分のみを座標変換なしで扱ってもまあ大丈夫かも。
c
c 全体系の応力を局所系に座標変換 Sl=T^T Sg T
c       subroutine ttsgt(t,sg,sl)
c       implicit real*8(a-h,o-z)
c       dimension t(14,14),sg(3,3),sl(3,3),c(3,3)
c       do 10 i=1,3
c       do 10 k=1,3
c        c(i,k)=0.d0
c       do 10 j=1,3
c        c(i,k)=c(i,k)+sg(i,j)*t(j,k)
c10     continue
c       do 20 i=1,3
c       do 20 k=1,3
c        sl(i,k)=0.d0
c       do 20 j=1,3
c        sl(i,k)=sl(i,k)+t(j,i)*c(j,k)
c20     continue
c       return
c       end
c
c 局所系の応力を全体系に座標変換 Sg=T Sl T^T
c      subroutine tsltt(t,sl,sg)
c       implicit real*8(a-h,o-z)
c       dimension t(14,14),sl(3,3),sg(3,3),c(3,3)
c       do 10 i=1,3
c       do 10 k=1,3
c        c(i,k)=0.d0
c       do 10 j=1,3
c        c(i,k)=c(i,k)+sl(i,j)*t(k,j)
c10     continue
c       do 20 i=1,3
c       do 20 k=1,3
c        sg(i,k)=0.d0
c       do 20 j=1,3
c        sg(i,k)=sg(i,k)+t(i,j)*c(j,k)
c20     continue
c       return
c       end
c
c
c そりねじり定数
c 正規化そり関数の級数について和を取る
         subroutine sori(kyuu,a,b,x,y,pi,w)
          implicit real*8 (a-h, o-z)
          w=0.d0
          do 10 n=1,kyuu,2
           w=w+(-1.d0)**( (n+1)/2 )/ n**3 *
     &       sinh(n*pi*y/2.d0/a) / cosh(n*pi*b/2.d0/a)
     &       *sin(n*pi*x/2.d0/a)
10        continue
           w=32.d0*a**2/pi**3* w + x*y
          return
          end
c
cc ねじり定数も断面積分する場合
cc ねじり定数
cc Timoshenko が計算した正解
cc 手で積分して、級数も極限を取れる部分は取ってある
cc         subroutine timj(kyuu,a,b,pi,tj)
cc          implicit real*8 (a-h, o-z)
cc           tj=0.d0
cc          do 10 n=1,kyuu,2
cc           tj=tj+1.d0/n**5*tanh(n*pi*b/2.d0/a)
cc10        continue
cc           tj=1.d0-192.d0/pi**5*a/b*tj
cc          return
cc          end
c
c
c 梁の線形剛性行列（無次元化してある）
c 2節点14自由度
        subroutine smx(el,eiy,pxz,pyz,c11,c22,c33,c44,c55,s)
        implicit real*8(a-h,o-z)
        dimension s(14,14)
        s(3,3)=   c11*el*el/eiy
        s(10,3)= -c11*el*el/eiy
        s(10,10)= c11*el*el/eiy
        s(2,2)=            12.d0*c33/eiy /(1.d0+12.d0*pyz)
        s(4,2)=            -6.d0*c33/eiy /(1.d0+12.d0*pyz)
        s(4,4)= (4.d0+12.d0*pyz)*c33/eiy /(1.d0+12.d0*pyz)
        s(11,2)=           -6.d0*c33/eiy /(1.d0+12.d0*pyz)
        s(11,4)=(2.d0-12.d0*pyz)*c33/eiy /(1.d0+12.d0*pyz)
        s(1,1)= 12.d0           *c22/eiy /(1.d0+12.d0*pxz)
        s(5,1)=  6.d0           *c22/eiy /(1.d0+12.d0*pxz)
        s(5,5)= (4.d0+12.d0*pxz)*c22/eiy /(1.d0+12.d0*pxz)
        s(12,1)= 6.d0           *c22/eiy /(1.d0+12.d0*pxz)
        s(12,5)=(2.d0-12.d0*pxz)*c22/eiy /(1.d0+12.d0*pxz)
        s(6,6)=  12.d0*c44/eiy/el/el+c55/eiy*6.d0/5.d0
        s(13,6)=-12.d0*c44/eiy/el/el-c55/eiy*6.d0/5.d0
        s(7,6)=  -6.d0*c44/eiy/el/el-c55/eiy/10.d0
        s(14,6)= -6.d0*c44/eiy/el/el-c55/eiy/10.d0
        s(7,7)=   4.d0*c44/eiy/el/el+c55/eiy*2.d0/15.d0
        s(14,7)=  2.d0*c44/eiy/el/el-c55/eiy/30.d0
        s(14,14)= 4.d0*c44/eiy/el/el+c55/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,dld14)
        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),dd3(14),dld14(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 dd3=ti^T * dd　　　　（^T は転置）
       call mxtv(14,ti,dd,dd3)
c この局所系 dd3 をオイラー化する前に、各ステップの
c 増分歪を計算するための増分変位（ステップごとの1発目の弧長+NR収束）
c dld14 を計算する。
       do i=1,14
        dld14(i)=dld14(i)+dd3(i)
       end do
c dd3 をオイラー角成分に変換
c dd=eti * dd3 = eti * ti^T * dd
       call mxv(14,eti,dd3,dd)
c オイラー角成分にしてから足し合わせてトータルな d を求める
       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
c stmx 内の上のddやdd3は既にtiの転置がかかった局所座標だけど、
c 最終的に解くddを全体系にするには、このskの右からもtiを掛ける
       call mxmxt(14,sk,ti,sngnt4)
       call mxmx(14,ti,sngnt4,sl)
c
c 不釣合力 b の計算（なぜトータルに求めてはいけないか）
c
c 載荷履歴に依存する弾塑性を解くには、トータルなdに対して
c 不釣り合い力を求めることはできない（f=K(d)d を材料非線形に対しても
c dの関数として与えることができてdの増分が取れるならできるかも
c 知れない？けど）ので、ここではΣK(Δd)Δd をbに累積しておいて、
c 行番号1200 の下辺りで、目標荷重pfからbを引いて不釣合力を求める。
c K(d)d≒ΣK(Δd)Δd が成り立つには、Δd が十分に小さい必要があり、
c そのためには、増分荷重Δfを、そのためには増分荷重Δrをじゅうぶんに
c 小さく与える必要がある。つまり、この部分でアップデートラグラジアン的
c な解き方となる。
c
c もし、不釣合力をトータルなまま計算したとすると、例えば
c f
c ↑ ／
c  /   →d  のようなバイリニアなモデルを与えても、
c
c f
c ↑ __／
c  /   →d  みたいな答えになってしまう。これは傾きがK1からK2に
c 変わったときに、pf-K2(d)で不釣合力を計算しているために、答えが、
c 原点を通る f=k2(d) の解に修正されてしまうからである。
c こうならないようにするには、
c pf-(K1Δd+K1Δd+K1Δd+...)+(K2Δd+K2Δd+...) のように履歴が
c 累積される形で不釣合力を求めなければならない。
c
c
c 上記のslには右から既にtiの転置が掛かってしまっているので、
c dd3 にtiを掛けて全体系にする
       call mxv(14,ti,dd3,dmy1)
       call mxv(14,sl,dmy1,dmy2)
c
       do i=1,14
        b(i)=b(i)+dmy2(i)
       end do
c
c         print *, 'ti'
c         print '(1p7d11.3)', ( (ti(i,j),j=8,14),i=8,14 )
c         print *, 'sl'
c         print '(1p7d11.3)', ( (sl(i,j),j=8,14),i=8,14 )
c         print *, 'dl'
c         print '(1p7d11.3)',  (dl(j),j=8,14)
c         print *, 'b'
c         print '(1p7d11.3)',  (b(j),j=1,7)
c         print '(1p7d11.3)',  (b(j),j=8,14)
c
c
cc 弾性問題をトータルラグラジアンとして解く場合
cc       call mxtv(14,t,dl,dmy2)
cc      call mxv(14,ss,dmy2,dmy1)
cc       call mtv14(r7x14,dmy1,dmy2)
cc       call mxtv(14,eti,dmy2,dmy1)
cc       call mxv(14,ti,dmy1,dmy2)
cc       do 550 i=1,14
cc トータルに求める場合は+dmy2ではなく-dmy2となる
cc 不釣り合い力の初期化のところも参照
cc       b(i)=b(i)-dmy2(i)
cc550    continue
c         print *, 'b-total'
c         print '(1p7d11.3)',  (b(j),j=1,7)
c         print '(1p7d11.3)',  (b(j),j=8,14)
c         iprint=iprint+1
c         if(iprint.eq.2) stop
       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
c 収束時のddの微少量として近似した座屈モードは節点辺り7自由度だけど
c トータルなdはオイラー角の関係で節点辺り14自由度で覚えさせてるから
c そこを重複して書かせる。
         subroutine modei7(i7,nofe,d,gmode)
         implicit real*8(a-h,o-z)
         dimension d(14*256)
          read(9,*) idummy,dmode
          i=0
          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 nsgm番目の要素の軸応力sgmをgnuplot作図用にファイル出力
       subroutine sgmgnu(nsgm,nhaba,ntaka,sgm)
       implicit real*8(a-h,o-z)
       dimension sgm(256,200,200)
        open(10,file='gnusgm.out')
c 縦横が200*200とかの分割数をsplotで描くと線で塗りつぶされるので
c 適当に間引くステップ幅を nstep で指定
          nstep=9
          do i=1,nhaba,nstep
           do j=1,ntaka,nstep
            write(10,560) i,j,sgm(nsgm,i,j)
           end do
           write(10,*)
          end do
          close(10)
         return
560      format( 2i4,1pd11.3)
         end
c
