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 ２ステップ目以降の計算に使う弧長増分：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 トータルな変位は、要素どうしが繋がる２節点のオイラー角を
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 接線剛性行列を１要素ずつ計算して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