c mozi-koodo: EUC 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 05/6/28版 c プログラムを主に作った人:後藤文彦 c 増分式を Mathematica で計算してくれた人:小林裕さん c いくつかの代数計算やスカイライン法のサブルーチンの c 元プログラムの作者:岩熊先生(http://www.civil.tohoku.ac.jp/~bear/) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c このプログラムは無保証です。 c このプログラムは改造/再配布して構いません。 c 但し、私的利用に留まらない場合は、もとのプログラムが c http://www.str.ce.akita-u.ac.jp/~gotou/programoj/ c から取ってきたものと分かるように明示して下さい。 c 改造/再配布後も上記の条件を保持して下さい。 c バグの報告は、 c http://www.str.ce.akita-u.ac.jp/~gotou/png/meeru.png c にお願いします。 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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=(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 一番最初の一発目の増分ステップの弧長:r read(*,*) r c 2ステップ目以降の計算に使う弧長増分:rn read(*,*) rn c c c ############# ここまでがデータ入力 ############# c c 収束精度:収束時の荷重が、N-R反復における一つ前のステップの c 荷重と有効数字何桁まで合うようにしたいかに応じて cnvr の小ささを c 決める cnvr=5.d-8 c 外乱の大きさに応じて収束精度を甘くしたいときとかは、 cc if(gmode.gt.1.d-10) cnvr=1.d-3 c みたいなのを入れたりするとか。 c c 三次元梁の線形剛性行列(反りを含む1節点7自由度)を s(14,14) c に入れる call smx(el,ea,eix,eiy,gj,eiw,s) c c ### 初期条件 ### c c 変位増分ベクトル:dd(自由度数+1) c 変位増分の回転角成分は空間固定3軸回りの微少回転角 c 変位ベクトル:d(節点自由度数*要素数) c トータルな変位は、要素どうしが繋がる2節点のオイラー角を c 別々に覚えておく。 c ddl は dd を各増分ステップの繰り返しで総和したもの c dp は荷重増分(弧長法なので、増分変位と同様に解かれる未知変数) do 100 i=1,n7 dd(i)=0.d0 d(i)=0.d0 ddl(i)=0.d0 100 continue dp=0.d0 p=0.d0 ro=r dpl=r c まだ一度も収束値が求まっていないことを表す:nccl(後述) nccl=0 c ### 各荷重ステップの一発目の弧長増分過程 ### 1000 continue nofcs=0 c 計算開始からの増分繰り返し数:ninc ninc=ninc+1 if(ninc.gt.4000) then write(*,*)'(増分繰り返しが 4000 を越えました)' close(8) stop endif c c 全変位増分ddと全変位dの中から、1要素ぶんずつをそれぞれ c dd14とd14に抜き出す do 1100 i=0,nofe-1 do 120 j=1,14 dd14(j)=dd(i*7+j) d14(j)=d(i*14+j) 120 continue do 130 j=1,4 lg4(j)=lorg(i*4+j) 130 continue c dd14とd14を使って1要素ずつ接線剛性行列を計算していく c 求まった1要素ぶんの接線剛性行列はskに入る call stmx(nofe,i,alp,gmm,dd14,d14,sk,b14,lg4,s) do 140 j=1,14 d(i*14+j)=d14(j) 140 continue c sk をスカイライン法の1次元ベクトルに書き込む c stuh に上三角、stlh に下三角が入る call skyln(nofe,i,sk,sko,stuh,stlh,imx) 1100 continue c do 150 i=1,n7 i8=imx(n8)+n8-i stuh(i8)=-f(i) stlh(i8)=ddl(i)/ro ddl(i)=0.d0 b(i)=0.d0 150 continue do 160 i=1,n8 stlh(imx(i))=b(i) 160 continue stuh(imx(n8))=dpl/ro stlh(imx(n8))=r ro=r dpl=0.d0 c 境界条件を入れる call bncn(nofe,bc,stuh,stlh,imx) c 代数方程式をガウスの消去法で解いて変位ベクトルを求め、 c その際に対角項の符号から、行列の負の固有値の数 jdt を求める call solve(stuh,stlh,dp,imx,msl,nnegtm,anegtm, & det,idt,jdt,n8) do 170 i=1,n8 dd(i)=stlh(imx(i)) 170 continue c c ### ニュートン−ラフソン 収束過程 ### c(基本的に上の「一発目の弧長増分過程」の計算と同じ) 2000 continue c 各増分ステップでの収束までの繰り返し数:nofcs nofcs=nofcs+1 do 180 i=1,n7 ddl(i)=ddl(i)+dd(i) ddi2=ddi2+dd(i)**2 ddl2=ddl2+ddl(i)**2 180 continue dpl=dpl+dp p=p+dp arc=ddi2/ddl2 ddi2=0 ddl2=0 if(nofcs.gt.6) then write(*,*)'反復数が6を越えても収束しないので弧長を半分にする' c 反復数が6を越えて辛うじて収束するようなことが何回も重なると、 c 答えがおかしくなってたりするので、その回数を nlet で数えて、 c 適当な回数に達したら強制的に終了して答えを出力するように c 荷重が収束した(p1=p)したことにしてやる。 nlet=nlet+1 if(nlet.gt.10) then p1=p goto 4004 endif goto 4001 endif if(arc.lt.tol) goto 3000 do 190 i=1,n7 b(i)=p*f(i) 190 continue do 1200 i=0,nofe-1 do 200 j=1,14 dd14(j)=dd(i*7+j) d14(j)=d(i*14+j) 200 b14(j)=b(i*7+j) do 210 j=1,4 lg4(j)=lorg(i*4+j) 210 continue c 接線剛性行列を1要素ずつ計算してskに入れる call stmx(nofe,i,alp,gmm,dd14,d14,sk,b14,lg4,s) do 220 j=1,14 d(i*14+j)=d14(j) 220 continue c skをスカイライン法の1次元配列に書き込む c stuhが上三角行列、stlhが下三角行列 call skyln(nofe,i,sk,sko,stuh,stlh,imx) do 230 j=1,14 b(i*7+j)=b14(j) 230 continue 1200 continue do 240 i=1,n7 i8=imx(n8)+n8-i stuh(i8)=-f(i) stlh(i8)=ddl(i) 240 continue do 250 i=1,n8 stlh(imx(i))=b(i) 250 continue stuh(imx(n8))=p stlh(imx(n8))=0.d0 c 境界条件を入れる call bncn(nofe,bc,stuh,stlh,imx) c 座屈点の探査のとき、荷重が前の増分ステップと誤差範囲に c になったら終了する。モード計算用の出力が必要な場合(modeo=1)は c fort.7 に収束時の接線剛性行列をフォーマットなしで書き出す。 if( (abs((p1-p)/p).lt.cnvr).and.(gmode.lt.1.d-20) ) then if(modeo.eq.1) then write(7) nofe+1,stuh,stlh,imx,imx9 rewind 7 endif close(8) c c 出力した接線行列の固有値解析をせずに、最後のN-R反復で求まる c (微少な)変位増分を固有ベクトル(座屈モード)の近似として c 利用する場合のために、一応、近似モードもここで出力しておく。 c こういう近似ができるということは、 c 久田俊明・野口裕久 著『非線形有限要素法の基礎と応用』 c(丸善株式会社)p.317とかに書いてある。 c call modeout(nofe,dd) c c 但し、接線剛性方程式の回転自由度は空間固定3軸回りの成分なので、 c 座屈モードの回転成分もその成分に対応して求まる。 c 座屈後へ分岐させる外乱としては、3軸回り成分をオイラー角に与えても c 十分に用を足すとは思うが。 c オイラー角成分がほしいときは、変換行列をかけるとかして。 c stop endif c c c 接線剛性方程式を解く call solve(stuh,stlh,dp,imx,msl,nnegtm,anegtm, & det,idt,jdt,n8) do 260 i=1,n8 dd(i)=stlh(imx(i)) 260 continue goto 2000 c c ### 出力 ### c 3000 continue arc=0.d0 c 収束しないために弧長が小さくなり過ぎたのを、また大きくしたい c ときは、nlet の大きさに応じて弧長を大きくしてやるとか c if(nlet.gt.1) rn=rn*(2.d0)**nlet c ここを通過しているということは、増分ステップが収束したという c ことなので増分ステップ内の繰り返し数 nlet を0にする。 nlet=0 c 座屈点の探査は、座屈点を超えたら、越える前の荷重ステップまで c 戻り、弧長を半分にして進み、また行き過ぎたらまた戻って弧長を c 半分にして次の荷重ステップを計算して座屈点に近づいていく。 c この座屈点を跨いだいったりきたりの計算において、同じ点を重複して計算 c しないように、「アキレスと亀」方式で弧長を毎回半分にしていく。 c c 接線剛性行列の負の固有値数 jdt と固有ベクトルとして得られる c 座屈モードの次数とはだいたい対応していることから、jdt が c 指定した座屈モードの次数(mode)を越えたかどうかで座屈点を c 越えたかどうかを判断する。座屈点を超えたかどうかにも c 様々な場合があり得るので、その状況をnacl,nbcl,ncclを使って判断する。 c c 一発目の弧長を大きく取りすぎて、最初から座屈点を超えてしまった場合 c は、正しい答えがまだ求まっていないので、最初から計算をやり直す c 必要がある。そのため、まだ一つも収束値が求まっていない時は c nccl=0, 一つでも収束値が(座屈点を超える前に)求まったら c nccl=1 とする。 if(jdt.lt.mode) nccl=1 write(*,*) 'N-R反復数=',nofcs write(*,570) det, idt,jdt 570 format(' ','det=',f11.3,'d',i4,' , ','負の固有値数:',i4) goto (571,572,573), nhikak 571 write(*,581) ' p=',p,' PL/EA=',p*ell/ea,' Euler=',eul goto 579 572 write(*,581) ' p=',p,' ML^2/8EI=',p*ell**2/8.d0/eix write(*,580) ' 修正Vlasov(正)=',vlp,' 修正Vlasov(負)=',vln goto 579 573 write(*,581)' p=',p,' PL^3/3EI=',p*ell**3/3./eix,' Tra=',tra goto 579 579 continue c c write(*,580)' p=',p,' Ml/EIy=',p*ell/e/yi 580 format(' ',a,1pd13.5,a,1pd13.5) 581 format(' ',a,1pd13.5,a,1pd13.5,a,1pd13.5) c 節点変位は、u, v, w, α, γ, φ, λの順に出力。 c α, γ, φ はオイラー角、λはねじれ率 c d は無次元化されているので、出力するときは有次元化する c c 指定要素の左節点を出力する場合 if(noflo.ne.0) then do 460 i=1,noflo j=noutl(i)-1 write(*,590) noutl(i),'左',d(j*14+1)*el,d(j*14+2)*el, & d(j*14+3)*el,d(j*14+4),d(j*14+5),d(j*14+6),d(j*14+7)/el 460 continue endif c 指定要素の右節点を出力する場合 if(nofro.ne.0) then do 470 i=1,nofro j=noutr(i)-1 write(*,590) noutr(i),'右',d(j*14+8)*el,d(j*14+9)*el, & d(j*14+10)*el,d(j*14+11),d(j*14+12),d(j*14+13),d(j*14+14)/el 470 continue endif c 指定要素の左節点のuvwをファイル出力 if(noflof.ne.0) then do 480 i=1,noflof j=nuvwl(i)-1 write(8,560) p,d(j*14+1)*el,d(j*14+2)*el,d(j*14+3)*el 480 continue endif c 指定要素の右節点のuvwをファイル出力 if(nofrof.ne.0) then do 490 i=1,nofrof j=nuvwr(i)-1 write(8,560) p,d(j*14+8)*el,d(j*14+9)*el,d(j*14+10)*el 490 continue endif 560 format(1p4d11.3) 590 format(i3,a,1p6d11.3,1pd10.2) if(mode.eq.0) goto 4003 c 接線剛性行列の負の固有値数jdtが指定した座屈モード次数を c 越えたら、座屈したと見なす。 if(jdt.ge.mode) then write(*,*)'座屈したみでだよ................' c 最初に座屈点を超えたら、まずnbcl=1,nacl=0としておく。 nbcl=1 nacl=0 endif 4004 continue c c ### 外乱の入れ方 ### c 一発目の弧長増分ステップが収束した後に入れる。 c 予め座屈点を計算しておいて、一発目の弧長を座屈点の直前の c 分岐しやすい荷重で外乱が入るように調整する。 if (nccl.eq.1) then c c 「アキレスと亀」方式で座屈点に達したら、そこで外乱を読み込む c 場合は、 c if(abs((p1-p)/p).lt.1.d-1) then c とか。 if(abs(gmode).lt.1.d-20) goto 5002 c 但し、外乱は一度しか読み込まない if(iyonda.eq.1) goto 5002 write(*,*)'外乱読み込み中' c 外乱を読み込んだ後に収束精度を甘くしたいときは、 c cnvr=1.d-3 c とか。外乱を1度 読んだらiyonda=1 iyonda=1 c c call modein(nofe,d,gmode) c c rn=rn*1.5d0 endif 5002 continue p1=p c 座屈点を超えたために座屈前に戻って1ステップ以上、座屈点を超えない c 点を解いているとき(nbcl=0, nacl=1)は c 弧長は半分にするけど、座屈してない状態だから前のステップには c 戻らなくていい。rn=rn/2の後、4002まで飛んで1000に戻る。 if(nacl.eq.1) rn=rn/2.d0 c 座屈点を超えたばかりの状態(nbcl=1, nacl=0)のときは、 c 座屈前に戻って弧長を半分にする(それでも、また座屈した時は c 上でnbcl=1になってるので、またここを通過する)。 c 次にここを通ときに座屈前に戻っていることを教えるために c nbcl=0, nacl=1 とする(nbcl=nacl=1の組み合わせはない) if(nbcl.eq.1) then nbcl=0 nacl=1 goto 4001 endif c こういうふうに goto が交差した書き方をしたのは、 c 一旦ここでif文を閉じておいて、上から4001へ飛べるようにするため。 goto 4002 4001 continue read(11) dd,d,ddl,dp,dpl,ro,p rewind 11 rn=rn/2.d0 c 一発目の弧長が大きすぎて座屈してしまった場合(nccl=0)は、 c 弧長を一割小さくして初期状態から計算しなおす。 if(nccl.eq.0) then r=r*0.9d0 ro=r dpl=ro endif 4002 continue c 一度でも収束値が求まっている(nccl=1)なら、収束するたびに、 c 変位をfort.11にフォーマットなしで書き出す。但し、 c 座屈点を超えた場合は、上記のread(11)で、座屈前の変位を読み込んだ c ものが、ここで書き出される。つまり、ここでは常に座屈していない c 状態での最も後に計算した変位が出力される。そして、次に座屈点を c 越えたときに、この変位の点まで戻れるようにしておく。 if(nccl.eq.1) then write(11) dd,d,ddl,dp,dpl,ro,p rewind 11 r=rn endif 4003 continue write(*,*)'----------------------------------------------------' c ここまでで一つの弧長増分ステップの終わり。次の弧長増分ステップ c に戻る。 goto 1000 end c c ものすごく長いけど、ここまでがメインプログラム。 c 以下がサブルーチン c c 等曲げを受ける円弧アーチの横ねじれ座屈に対するウラソフの解析解 subroutine vlasov(n,blen,eix,eiy,gj,eiw,th,pi,vln,vlp) implicit real*8(a-h,o-z) if(th.lt.1.d-6) th=1.d-6 r=blen/th f=(n*pi/blen)**2*eiw+gj a=eiy/eix-1.d0+(1.d0/eix-eiy/eix/eix)*f b=(eiy+(1.d0-2.d0*eiy/eix)*f)/r c=((n*pi/blen)**2-1.d0/r/r)*eiy*f d=b**2-4.d0*a*c vln=(-b+sqrt(d))/2.d0/a vlp=(-b-sqrt(d))/2.d0/a return end c c 梁の線形剛性行列(無次元化してある) c 2節点14自由度 subroutine smx(el,ea,eix,eiy,gj,eiw,s) implicit real*8(a-h,o-z) dimension s(14,14) s(3,3)= ea*el*el/eiy s(10,3)= -ea*el*el/eiy s(10,10)= ea*el*el/eiy s(2,2)= 12.d0*eix/eiy s(4,2)= -6.d0*eix/eiy s(4,4)= 4.d0*eix/eiy s(11,2)= -6.d0*eix/eiy s(11,4)= 2.d0*eix/eiy s(1,1)= 12.d0 s(5,1)= 6.d0 s(5,5)= 4.d0 s(12,1)= 6.d0 s(12,5)= 2.d0 s(6,6)= 12.d0*eiw/eiy/el/el+gj/eiy*6.d0/5.d0 s(13,6)=-12.d0*eiw/eiy/el/el-gj/eiy*6.d0/5.d0 s(7,6)= -6.d0*eiw/eiy/el/el-gj/eiy/10.d0 s(14,6)= -6.d0*eiw/eiy/el/el-gj/eiy/10.d0 s(7,7)= 4.d0*eiw/eiy/el/el+gj/eiy*2.d0/15.d0 s(14,7)= 2.d0*eiw/eiy/el/el-gj/eiy/30.d0 s(14,14)= 4.d0*eiw/eiy/el/el+gj/eiy*2.d0/15.d0 s(9,2)=-s(2,2) s(9,9)=s(2,2) s(9,4)=-s(4,2) s(11,9)=s(9,4) s(11,11)=s(4,4) s(8,1)=-s(1,1) s(8,5)=-s(5,1) s(8,8)=s(1,1) s(12,8)=-s(5,1) s(12,12)=s(5,5) s(13,7)=-s(7,6) s(13,13)=s(6,6) s(14,13)=s(13,7) c do 100 i=1,13 do 100 j=i+1,14 s(i,j)=s(j,i) 100 continue c return end c c 接線剛性行列(変位に関する増分式)の計算 c SUBROUTINE STMX(nofe,is,alp,gmm,dd,d,sl,b,lg,ss) implicit real*8 (a-h,o-z) DIMENSION ss(14,14),sk(14,14),sl(14,14),d(14),dl(14),dd(14), & t(14,14),ti(14,14),t14(14,14),eti(14,14),alp(256),gmm(256), & cf1(7,7),cf2(7,7),t1(7,7),ta(7,7),tg(7,7),tp(7,7),b(15),lg(4), & deti(14,14,14),r7x14(7,14),dr7x14(14,7,14),dss(14,14,14), & dt0(7,14,14),dmy1(14),dmy2(14), & tngnt1(14,14),tngnt3(14,14),tngnt4(14,14),sngnt4(14,14) call cl4rms(14,t14, 14,sl, 14,sk, 14,ti) call cl4rms(7,ta, 7,tg, 7,tp, 7,cf1) call cl4rms(7,cf2, 14,eti, 7,t1, 14,t) call clrv(14,dl) call cl4rms(14,tngnt1,14,tngnt3,14,tngnt4,14,sngnt4) call cl5rms(deti,r7x14,dr7x14,dss,dt0) call clrmx(14,sngnt4) c c 曲がり梁の要素の初期方向を与える座標変換:ti call trns(ti,alp(is+1),gmm(is+1),0.d0) c 3軸回りの回転成分をオイラー角成分に変換する行列:eti call euler(eti,d(4),d(5),d(11),d(12)) c 端部の境界条件が局所系で与えられていた時は、その変位に c 掛けられるtiを外す(単位行列に変える) c 左端の並進変位 if(lg(1).eq.1) then do 1 i=1,3 do 1 j=1,3 ti(i,j)=0.d0 ti(i,i)=1.d0 1 continue endif c 左端の回転変位 if(lg(2).eq.1) then do 2 i=4,6 do 2 j=4,6 ti(i,j)=0.d0 ti(i,i)=1.d0 2 continue endif c 右端の並進変位 if(lg(3).eq.1) then do 3 i=8,10 do 3 j=8,10 ti(i,j)=0.d0 ti(i,i)=1.d0 3 continue endif c 右端の回転変位 if(lg(4).eq.1) then do 4 i=11,13 do 4 j=11,13 ti(i,j)=0.d0 ti(i,i)=1.d0 4 continue endif c 変位増分ddをti系のオイラー角成分に変換する c dmy1=ti^T * dd    (^T は転置) call mxtv(14,ti,dd,dmy1) c dd=eti * dmy1 = eti * ti^T * dd call mxv(14,eti,dmy1,dd) c オイラー角成分にしてから足し合わせる do 80 i=1,14 d(i)=d(i)+dd(i) 80 continue c オイラー角が更新されたところで、座標変換行列を再定義:t call trns(t,d(4),d(5),d(6)) c 3軸成分をオイラー角成分に変換する行列も再定義:eti call euler(eti,d(4),d(5),d(11),d(12)) c 変位ベクトル u1=d(1) v1=d(2) w1=d(3) a1=d(4) g1=d(5) p1=d(6) q1=d(7) u2=d(8) v2=d(9) w2=d(10) a2=d(11) g2=d(12) p2=d(13) q2=d(14) ca1=cos(a1) cg1=cos(g1) cp1=cos(p1) ca2=cos(a2) cg2=cos(g2) cp2=cos(p2) sa1=sin(a1) sg1=sin(g1) sp1=sin(p1) sa2=sin(a2) sg2=sin(g2) sp2=sin(p2) c c 相対変位 dl(7)=q1 dl(8)=u2-u1-ca1*sg1 dl(9)=v2-v1-sa1 dl(10)=w2-w1-(ca1*cg1-1.d0) gg=(g1+g2)/2.d0 aa=(a1+a2)/2.d0 pp=(p1+p2)/2.d0 cpp=cos(pp) caa=cos(aa) cgg=cos(gg) spp=sin(pp) saa=sin(aa) sgg=sin(gg) dl(11)=-cgg*(a2-a1)+sgg*caa*(p2-p1) dl(12)=(g2-g1)+saa*(p2-p1) dl(13)=sgg*(a2-a1)+cgg*caa*(p2-p1) dl(14)=q2 c c t の増分 ta(1,1)=-sp1*ca1*sg1 ta(1,2)=-cp1*ca1*sg1 ta(1,3)=-sa1*sg1 ta(2,1)=-sp1*sa1 ta(2,2)=-cp1*sa1 ta(2,3)=ca1 ta(3,1)=-sp1*ca1*cg1 ta(3,2)=-cp1*ca1*cg1 ta(3,3)=-sa1*cg1 tg(1,1)=-cp1*sg1-sp1*sa1*cg1 tg(1,2)=sp1*sg1-cp1*sa1*cg1 tg(1,3)=ca1*cg1 tg(3,1)=-cp1*cg1+sp1*sa1*sg1 tg(3,2)=sp1*cg1+cp1*sa1*sg1 tg(3,3)=-ca1*sg1 tp(1,1)=-sp1*cg1-cp1*sa1*sg1 tp(1,2)=-cp1*cg1+sp1*sa1*sg1 tp(2,1)=cp1*ca1 tp(2,2)=-sp1*ca1 tp(3,1)=sp1*sg1-cp1*sa1*cg1 tp(3,2)=cp1*sg1+sp1*sa1*cg1 do 162 i=1,3 do 162 j=1,3 ta(3+i,3+j)=ta(i,j) tg(3+i,3+j)=tg(i,j) tp(3+i,3+j)=tp(i,j) 162 continue c do 100 i=1,6 do 100 j=1,6 dt0(4,i,j)=ta(i,j) dt0(4,i+7,j+7)=ta(i,j) dt0(5,i,j)=tg(i,j) dt0(5,i+7,j+7)=tg(i,j) dt0(6,i,j)=tp(i,j) dt0(6,i+7,j+7)=tp(i,j) 100 continue c c 相対変位の増分 cf1(1,1)=-1.d0 cf1(1,4)=sa1*sg1 cf1(1,5)=-ca1*cg1 cf1(2,2)=-1.d0 cf1(2,4)=-ca1 cf1(3,3)=-1.d0 cf1(3,4)=sa1*cg1 cf1(3,5)=ca1*sg1 cf1(4,4)=cgg-(p2-p1)/2.d0*sgg*saa cf1(4,5)=(a2-a1)/2.d0*sgg+(p2-p1)/2.d0*cgg*caa cf1(4,6)=-sgg*caa cf1(5,4)=(p2-p1)/2.d0*caa cf1(5,5)=-1.d0 cf1(5,6)=-saa cf1(6,4)=-sgg-(p2-p1)/2.d0*cgg*saa cf1(6,5)=(a2-a1)/2.d0*cgg-(p2-p1)/2.d0*sgg*caa cf1(6,6)=-cgg*caa cf2(1,1)=1.d0 cf2(2,2)=1.d0 cf2(3,3)=1.d0 cf2(4,4)=-cgg-(p2-p1)/2.d0*sgg*saa cf2(4,5)=cf1(4,5) cf2(4,6)=sgg*caa cf2(5,4)=cf1(5,4) cf2(5,5)=1.d0 cf2(5,6)=saa cf2(6,4)=cf1(6,4)+2.d0*sgg cf2(6,5)=cf1(6,5) cf2(6,6)=cgg*caa cf2(7,7)=1.d0 c c その他、煩雑な増分式は、farejo の中で作る call farejo(d,dl,t,cf1,cf2,dt0,r7x14,dr7x14,deti) c call mxtv8(t,dl,dmy1) call mxv8(ss,dmy1,dmy2) do 140 i=1,14 do 140 j=1,14 tngnt1(i,j)=0.d0 do 140 k=1,7 c dr7x17 no naka no sori no zoubun ha 0 tngnt1(i,j)=tngnt1(i,j)+dr7x14(j,k,i)*dmy2(k+7) 140 continue c call mxm14(ss,r7x14,sngnt4) call mtmx14(r7x14,sngnt4,tngnt3) c do 260 i=1,14 do 260 j=1,14 sl(i,j)=tngnt1(i,j)+tngnt3(i,j) 260 continue c call mxmx(14,sl,eti,sk) call mxtmx(14,eti,sk,sl) c call mtv14(r7x14,dmy2,dmy1) c do 300 i=1,14 do 300 j=1,14 tngnt4(i,j)=0.d0 do 300 k=1,14 tngnt4(i,j)=tngnt4(i,j)+deti(j,k,i)*dmy1(k) 300 continue c call mxmx(14,tngnt4,eti,sngnt4) c do 290 i=1,14 do 290 j=1,14 sk(i,j)=sngnt4(i,j)+sl(i,j) 290 continue c call mxmxt(14,sk,ti,sngnt4) call mxmx(14,ti,sngnt4,sl) c c 不釣り合い力の計算 c call mxtv(14,t,dl,dmy2) call mxv(14,ss,dmy2,dmy1) call mtv14(r7x14,dmy1,dmy2) call mxtv(14,eti,dmy2,dmy1) call mxv(14,ti,dmy1,dmy2) do 550 i=1,14 b(i)=b(i)-dmy2(i) 550 continue return c888 format(14f5.1) end c c オイラー角で座標変換行列を定義する subroutine trns(t,al,gm,ph) IMPLICIT REAL*8 (A-H,O-Z) dimension t(14,14) c write(*,*)'I am in sub trns' t(1,1)=cos(ph)*cos(gm)-sin(ph)*sin(al)*sin(gm) t(1,2)=-sin(ph)*cos(gm)-cos(ph)*sin(al)*sin(gm) t(1,3)=cos(al)*sin(gm) t(2,1)=sin(ph)*cos(al) t(2,2)=cos(ph)*cos(al) t(2,3)=sin(al) t(3,1)=-cos(ph)*sin(gm)-sin(ph)*sin(al)*cos(gm) t(3,2)=sin(ph)*sin(gm)-cos(ph)*sin(al)*cos(gm) t(3,3)=cos(al)*cos(gm) do 10 i=1,3 do 10 j=1,3 t(i+3,j+3)=t(i,j) t(i+7,j+7)=t(i,j) t(i+10,j+10)=t(i,j) 10 continue t(7,7)=1.d0 t(14,14)=1.d0 return end c c 3軸回りの回転角成分をオイラー角成分に変換する subroutine euler(t,a1,g1,a2,g2) IMPLICIT REAL*8 (A-H,O-Z) dimension t(14,14) do 10 i=1,14 t(i,i)=1.d0 10 continue t(4,4)=-cos(g1) t(4,6)=sin(g1) t(5,4)=-sin(g1)*tan(a1) t(5,6)=-cos(g1)*tan(a1) t(6,4)=sin(g1)/cos(a1) t(6,6)=cos(g1)/cos(a1) t(11,11)=-cos(g2) t(11,13)=sin(g2) t(12,11)=-sin(g2)*tan(a2) t(12,13)=-cos(g2)*tan(a2) t(13,11)=sin(g2)/cos(a2) t(13,13)=cos(g2)/cos(a2) return end c c c 以下の mx...とかmt...という一連のサブルーチンは行列の c かけ算を行うためのもの。 c subroutine mtmx14(a,b,c) implicit real*8(a-h,o-z) dimension a(7,14),b(14,14),c(14,14) do 10 i=1,7 do 10 k=1,7 c(i,k)=0.d0 do 10 j=1,7 c(i,k)=c(i,k)+a(j,i)*b(j+7,k) c(i,k+7)=c(i,k+7)+a(j,i)*b(j+7,k+7) c(i+7,k)=c(i+7,k)+a(j,i+7)*b(j+7,k) c(i+7,k+7)=c(i+7,k+7)+a(j,i+7)*b(j+7,k+7) 10 continue do 20 j=1,14 c(7,j)=c(7,j)+b(7,j) 20 continue return end c c subroutine mxm14(a,b,c) implicit real*8(a-h,o-z) dimension a(14,14),b(7,14),c(14,14) do 10 i=1,7 do 10 k=1,7 c(i,k)=0.d0 do 10 j=1,7 c(i,k)=c(i,k)+a(i,j+7)*b(j,k) c(i,k+7)=c(i,k+7)+a(i,j+7)*b(j,k+7) c(i+7,k)=c(i+7,k)+a(i+7,j+7)*b(j,k) c(i+7,k+7)=c(i+7,k+7)+a(i+7,j+7)*b(j,k+7) 10 continue do 20 i=1,14 c(i,7)=c(i,7)+a(i,7) 20 continue return end c subroutine mxmx(n,a,b,c) implicit real*8(a-h,o-z) dimension a(n,n),b(n,n),c(n,n) do 10 i=1,n do 10 k=1,n c(i,k)=0.d0 do 10 j=1,n c(i,k)=c(i,k)+a(i,j)*b(j,k) 10 continue return end c subroutine mxtmx(n,a,b,c) implicit real*8(a-h,o-z) dimension a(n,n),b(n,n),c(n,n) do 10 i=1,n do 10 k=1,n c(i,k)=0.d0 do 10 j=1,n c(i,k)=c(i,k)+a(j,i)*b(j,k) 10 continue return end c subroutine mxmxt(n,a,b,c) implicit real*8(a-h,o-z) dimension a(n,n),b(n,n),c(n,n) do 10 i=1,n do 10 k=1,n c(i,k)=0.d0 do 10 j=1,n c(i,k)=c(i,k)+a(i,j)*b(k,j) 10 continue return end c subroutine mxv(n,a,b,c) implicit real*8(a-h,o-z) dimension a(n,n),b(n),c(n) do 10 i=1,n c(i)=0.d0 do 10 j=1,n c(i)=c(i)+a(i,j)*b(j) 10 continue return end c subroutine mxtv(n,a,b,c) implicit real*8(a-h,o-z) dimension a(n,n),b(n),c(n) do 10 i=1,n c(i)=0.d0 do 10 j=1,n c(i)=c(i)+a(j,i)*b(j) 10 continue return end c subroutine mxv8(a,b,c) implicit real*8(a-h,o-z) dimension a(14,14),b(14),c(14) do 10 i=1,14 c(i)=0.d0 do 10 j=1,7 c(i)=c(i)+a(i,j+7)*b(j+7) 10 continue do 20 i=1,14 c(i)=c(i)+a(i,7)*b(7) 20 continue return end c subroutine mxtv8(a,b,c) implicit real*8(a-h,o-z) dimension a(14,14),b(14),c(14) do 10 i=1,14 c(i)=0.d0 do 10 j=1,7 c(i)=c(i)+a(j+7,i)*b(j+7) 10 continue do 20 i=1,14 c(i)=c(i)+a(7,i)*b(7) 20 continue return end c subroutine mtv14(a,b,c) implicit real*8(a-h,o-z) dimension a(7,14),b(14),c(14) do 10 i=1,14 c(i)=0.d0 do 10 j=1,7 c(i)=c(i)+a(j,i)*b(j+7) 10 continue c(7)=c(7)+b(7) return end c subroutine maxa(nofe,n8,imx) dimension imx(n8+1) c imx(1)=1 do 10 j=2,15 imx(j)=imx(j-1)+j-1 10 continue do 20 nj=1,nofe-1 do 20 j=9,15 imx(7*nj+j)=imx(7*nj+j-1)+j-1 20 continue imx(n8)=imx(n8-1)+14 imx(n8+1)=imx(n8)+n8 return end c c 各要素ごとの接線剛性行列を、スカイライン法の上三角行列と c 下三角行列にわけながら全要素について重ね合わせる subroutine skyln(nofe,i,s,so,stuh,stlh,imx) implicit real*8(a-h,o-z) dimension s(14,14),so(7,7), & stuh(84*nofe+36),stlh(84*nofe+36),imx(7*nofe+9) c do 10 j=1,7 do 10 k=j,7 ij=7*i+j ik=7*i+k im=imx(ik)+ik-ij stuh(im)=s(j,k)+so(j,k) stlh(im)=s(k,j)+so(k,j) so(j,k)=s(j+7,k+7) so(k,j)=s(k+7,j+7) 10 continue do 15 j=1,7 do 15 k=1,7 ij=7*i+j ik=7*i+k im7=imx(ik+7)+ik+7-ij stuh(im7)=s(j,k+7) stlh(im7)=s(k+7,j) 15 continue if(i.eq.nofe-1) then do 20 j=1,7 do 20 k=j,7 ij=7*i+j ik=7*i+k i7=imx(ik+7)+ik-ij stuh(i7)=s(j+7,k+7) stlh(i7)=s(k+7,j+7) so(j,k)=0.d0 so(k,j)=0.d0 20 continue endif return end c c 境界条件を入れる subroutine bncn(nofe,bc,stuh,stlh,imx) implicit real*8(a-h,o-z) dimension bc(7*nofe+8), & stuh(84*nofe+36),stlh(84*nofe+36),imx(7*nofe+9) c do 10 i=0,nofe-1 do 10 j=1,7 do 10 k=j,14 ij=7*i+j ik=7*i+k im=imx(ik)+ik-ij stuh(im)=bc(ij)*bc(ik)*stuh(im) stlh(im)=bc(ij)*bc(ik)*stlh(im) 10 continue do 15 j=7*nofe+1,7*nofe+7 do 15 k=j,7*nofe+7 im=imx(k)+k-j stuh(im)=bc(j)*bc(k)*stuh(im) stlh(im)=bc(j)*bc(k)*stlh(im) 15 continue n8=7*nofe+8 do 20 i=1,n8 if(bc(i).gt.1.d-6) goto 20 stuh(imx(i))=1.d0 stlh(imx(i))=0.d0 im=imx(n8)+n8-i stuh(im)=bc(i)*stuh(im) stlh(im)=bc(i)*stlh(im) 20 continue return end c c c 接線剛性方程式を解いて変位増分を求めるとともに、接線剛性行列の c 負の固有値数を求める。 c 岩熊先生が、東京大学の大計のライブラリーかどこかの「ガウスの消去法」の c サブルーチンを2次元骨組をスカイライン法で解くために書き換えたものを c 拝借した。 c subroutine solve(stuh,stlh,dp,imx,msl,nnegtm,anegtm, & det,idt,jdt,n01) c c 方程式は掃き出し法によりここで解かれる c implicit real*8 (a-h,o-z) dimension stuh(12*n01-60),stlh(12*n01-60) dimension imx(n01+1),msl(n01),nnegtm(n01),anegtm(n01) c c 行列式 = det D idt c jdt : 負の固有値の数 c do 100 i=1,n01 msl(i)=1+i+imx(i)-imx(i+1) 100 continue det=stuh(1) idt=0 jdt=0 if( det.gt.0.d0 ) go to 80 jdt=1 nnegtm(jdt)=1 anegtm(jdt)=stuh(1) 80 continue c c 前進消去 c do 10 j=2,n01 ist=msl(j)+1 ied=j-1 if( ist.gt.ied ) go to 20 do 30 i=ist,ied kst=msl(j) if( msl(i).gt.kst ) kst=msl(i) ked=i-1 ij=imx(j)+j-i do 30 k=kst,ked kj=imx(j)+j-k ki=imx(i)+i-k stuh(ij)=stuh(ij)-stuh(kj)*stlh(ki) stlh(ij)=stlh(ij)-stlh(kj)*stuh(ki) 30 continue 20 continue do 40 i=msl(j),ied ij=imx(j)+j-i stuh(ij)=stuh(ij)/stuh(imx(i)) stlh(ij)=stlh(ij)/stuh(imx(i)) 40 continue do 50 k=msl(j),ied kj=imx(j)+j-k stuh(imx(j))=stuh(imx(j))-stuh(kj)*stlh(kj)*stuh(imx(k)) stlh(imx(j))=stlh(imx(j))-stlh(kj)*stlh(imx(k)) 50 continue c if( j.eq.n01 ) go to 10 if( stuh(imx(j)).gt.0.d0 ) go to 90 jdt=jdt+1 nnegtm(jdt)=j anegtm(jdt)=stuh(imx(j)) 90 continue det=det*stuh(imx(j)) jjjdt=idint(dlog10(dabs(det))) idt=idt+jjjdt det=det/10.d0**(jjjdt) c 10 continue c c 後進消去 c do 60 j=1,n01 stlh(imx(j))=stlh(imx(j))/stuh(imx(j)) 60 continue do 70 j=2,n01 i=n01+2-j ked=i-1 do 70 k=msl(i),ked ii=imx(i)+i-k stlh(imx(k))=stlh(imx(k))-stuh(ii)*stlh(imx(i)) 70 continue dp=stlh(imx(n01)) c return end c c 以下のcl..というサブルーチン群は配列の中に 0 を書き込んで c 初期化するためのもの c subroutine clrv(n,a) implicit real*8(a-h,o-z) dimension a(n) do 10 i=1,n a(i)=0.d0 10 continue return end c subroutine cliv(n,l) implicit real*8(a-h,o-z) dimension l(n) do 10 i=1,n l(i)=0 10 continue return end c subroutine cl4rvs(n1,a1,n2,a2,n3,a3,n4,a4) implicit real*8(a-h,o-z) dimension a1(n1),a2(n2),a3(n3),a4(n4) do 10 i=1,n1 a1(i)=0.d0 10 continue do 20 i=1,n2 a2(i)=0.d0 20 continue do 30 i=1,n3 a3(i)=0.d0 30 continue do 40 i=1,n4 a4(i)=0.d0 40 continue return end c subroutine cl4ivs(n1,l1,n2,l2,n3,l3,n4,l4) implicit real*8(a-h,o-z) dimension l1(n1),l2(n2),l3(n3),l4(n4) do 10 i=1,n1 l1(i)=0 10 continue do 20 i=1,n2 l2(i)=0 20 continue do 30 i=1,n3 l3(i)=0 30 continue do 40 i=1,n4 l4(i)=0 40 continue return end c subroutine clrmx(n,a) implicit real*8(a-h,o-z) dimension a(n,n) do 10 i=1,n do 10 j=1,n a(i,j)=0.d0 10 continue return end c subroutine cl4rms(n1,a1,n2,a2,n3,a3,n4,a4) implicit real*8(a-h,o-z) dimension a1(n1,n1),a2(n2,n2),a3(n3,n3),a4(n4,n4) do 10 i=1,n1 do 10 j=1,n1 a1(i,j)=0.d0 10 continue do 20 i=1,n2 do 20 j=1,n2 a2(i,j)=0.d0 20 continue do 30 i=1,n3 do 30 j=1,n3 a3(i,j)=0.d0 30 continue do 40 i=1,n4 do 40 j=1,n4 a4(i,j)=0.d0 40 continue return end c subroutine cl5rms(deti,r7x14,dr7x14,dss,dt0) implicit real*8(a-h,o-z) dimension deti(14,14,14),r7x14(7,14),dr7x14(14,7,14), & dss(14,14,14),dt0(7,14,14) c do 10 i=1,14 do 10 j=1,14 do 10 k=1,14 deti(i,j,k)=0.d0 dss(i,j,k)=0.d0 10 continue c do 20 i=1,7 do 20 j=1,14 r7x14(i,j)=0.d0 do 20 k=1,14 dr7x14(j,i,k)=0.d0 dt0(i,j,k)=0.d0 20 continue return end c c c c 増分行列の幾つかがここで作られる c 尚、このサブルーチン内の増分式は、 c 小林裕さんが Mathematica で計算してくれました。 c subroutine farejo(d,dl,t,cf1,cf2,dt0,r7x14,dr7x14,deti) IMPLICIT REAL*8 (A-H,O-Z) dimension d(14),dl(14),t(14,14),cf1(7,7),cf2(7,7), & dt0(7,14,14),r7x14(7,14),dr7x14(14,7,14),deti(14,14,14), & a(6,3),da(14,6,3),dcf1(14,7,7),dcf2(14,7,7) dimension tdb1(14,7,7),tdb2(14,7,7),dtb1(14,7,7), & dtb2(14,7,7),dt1(14,14,14),t1(7,7) do 1007 i=1,6 do 1007 j=1,3 a(i,j)=0.d0 do 1007 k=1,14 da(k,i,j)=0.d0 1007 continue do 1013 k=1,14 do 1013 i=1,14 do 1013 j=1,14 deti(k,i,j)=0.d0 1013 continue do 1018 i=1,7 do 1018 j=1,14 r7x14(i,j)=0.d0 do 1018 k=1,14 dr7x14(k,i,j)=0.d0 1018 continue do 1024 k=1,14 do 1024 i=1,7 do 1024 j=1,7 dcf1(k,i,j)=0.d0 dcf2(k,i,j)=0.d0 1024 continue do 460 k=1,7 do 460 i=1,14 do 460 j=1,14 dt1(k,i,j)=dt0(k,i,j) dt1(7+k,i,j)=0.d0 460 continue do 470 i=1,7 do 470 j=1,7 t1(i,j)=t(i,j) 470 continue u1=d(1) v1=d(2) w1=d(3) a1=d(4) g1=d(5) p1=d(6) q1=d(7) u2=d(8) v2=d(9) w2=d(10) a2=d(11) g2=d(12) p2=d(13) q2=d(14) c ca1=cos(a1) cg1=cos(g1) cp1=cos(p1) ca2=cos(a2) cg2=cos(g2) cp2=cos(p2) sa1=sin(a1) sg1=sin(g1) sp1=sin(p1) sa2=sin(a2) sg2=sin(g2) sp2=sin(p2) ta1=tan(a1) ta2=tan(a2) tg1=tan(g1) tg2=tan(g2) tp1=tan(p1) tp2=tan(p2) c gmm=(g1+g2)/2.d0 all=(a1+a2)/2.d0 phh=(p2-p1)/2.d0 allp=(a2-a1)/2.d0 c cgg=cos(gmm) caa=cos(all) cph=cos(phh) callp=cos(allp) sgg=sin(gmm) saa=sin(all) sph=sin(phh) sallp=sin(allp) c ru=dl(8) rv=dl(9) rw=dl(10) rx=dl(11) ry=dl(12) rz=dl(13) c dcf1(4,1,4)=ca1*sg1 dcf1(5,1,4)=sa1*cg1 dcf1(4,1,5)=sa1*cg1 dcf1(5,1,5)=ca1*sg1 dcf1(4,2,4)=sa1 dcf1(4,3,4)=ca1*cg1 dcf1(5,3,4)=-sa1*sg1 dcf1(4,3,5)=-sa1*sg1 dcf1(5,3,5)=ca1*cg1 dcf1(4,4,4)=-phh*sgg*caa/2.d0 dcf1(11,4,4)=dcf1(4,4,4) dcf1(5,4,4)=(-sgg-phh*cgg*saa)/2.d0 dcf1(12,4,4)=dcf1(5,4,4) dcf1(6,4,4)=sgg*saa/2.d0 dcf1(13,4,4)=-dcf1(6,4,4) dcf1(4,4,5)=(-sgg-phh*cgg*saa)/2.d0 dcf1(5,4,5)=(allp*cgg-phh*sgg*caa)/2.d0 dcf1(6,4,5)=-cgg*caa/2.d0 dcf1(11,4,5)=(sgg-phh*cgg*saa)/2.d0 dcf1(12,4,5)=dcf1(5,4,5) dcf1(13,4,5)=-dcf1(6,4,5) dcf1(4,4,6)=sgg*saa/2.d0 dcf1(5,4,6)=-cgg*caa/2.d0 dcf1(11,4,6)=dcf1(4,4,6) dcf1(12,4,6)=dcf1(5,4,6) dcf2(4,4,4)=-phh*sgg*caa/2.d0 dcf2(5,4,4)=(sgg-phh*cgg*saa)/2.d0 dcf2(6,4,4)=sgg*saa/2.d0 dcf2(11,4,4)=dcf2(4,4,4) dcf2(12,4,4)=dcf2(5,4,4) dcf2(13,4,4)=-dcf2(6,4,4) dcf1(4,5,4)=-saa*phh/2.d0 dcf1(11,5,4)=dcf1(4,5,4) dcf1(4,5,6)=-caa/2.d0 dcf1(11,5,6)=dcf1(4,5,6) dcf1(6,5,4)=-caa/2.d0 dcf1(13,5,4)=-dcf1(6,5,4) dcf1(4,6,4)=-phh*cgg*caa/2.d0 dcf1(5,6,4)=(-cgg+phh*sgg*saa)/2.d0 dcf1(6,6,4)=cgg*saa/2.d0 dcf1(11,6,4)=dcf1(4,6,4) dcf1(12,6,4)=dcf1(5,6,4) dcf1(13,6,4)=-dcf1(6,6,4) dcf1(4,6,5)=(-cgg+phh*sgg*saa)/2.d0 dcf1(5,6,5)=(-allp*sgg-phh*cgg*caa)/2.d0 dcf1(6,6,5)=sgg*caa/2.d0 dcf1(11,6,5)=(cgg+phh*sgg*saa)/2.d0 dcf1(12,6,5)=dcf1(5,6,5) dcf1(13,6,5)=-dcf1(6,6,5) dcf1(4,6,6)=cgg*saa/2.d0 dcf1(5,6,6)=sgg*caa/2.d0 dcf1(11,6,6)=dcf1(4,6,6) dcf1(12,6,6)=dcf1(5,6,6) dcf2(4,6,4)=-phh*cgg*caa/2.d0 dcf2(5,6,4)=(cgg+phh*sgg*saa)/2.d0 dcf2(6,6,4)=cgg*saa/2.d0 dcf2(11,6,4)=dcf2(4,6,4) dcf2(12,6,4)=dcf2(5,6,4) dcf2(13,6,4)=-dcf2(6,6,4) dcf2(4,6,5)=dcf1(4,6,5) dcf2(5,6,5)=dcf1(5,6,5) dcf2(6,6,5)=dcf1(6,6,5) dcf2(11,6,5)=dcf1(11,6,5) dcf2(12,6,5)=dcf1(12,6,5) dcf2(13,6,5)=dcf1(13,6,5) dcf2(4,4,6)=-sgg*saa/2.d0 dcf2(5,4,6)=cgg*caa/2.d0 dcf2(11,4,6)=dcf2(4,4,6) dcf2(12,4,6)=dcf2(5,4,6) dcf2(4,4,5)=dcf1(4,4,5) dcf2(5,4,5)=dcf1(5,4,5) dcf2(6,4,5)=dcf1(6,4,5) dcf2(11,4,5)=dcf1(11,4,5) dcf2(12,4,5)=dcf1(12,4,5) dcf2(13,4,5)=dcf1(13,4,5) dcf2(4,5,4)=dcf1(4,5,4) dcf2(6,5,4)=dcf1(6,5,4) dcf2(11,5,4)=dcf1(11,5,4) dcf2(13,5,4)=dcf1(13,5,4) dcf2(4,5,6)=caa/2.d0 dcf2(11,5,6)=dcf2(4,5,6) dcf2(4,6,6)=-cgg*saa/2.d0 dcf2(11,6,6)=-cgg*saa/2.d0 dcf2(5,6,6)=-sgg*caa/2.d0 dcf2(13,5,5)=-dcf2(6,5,5) dcf2(12,6,6)=-sgg*caa/2.d0 c deti(5,4,4)=sg1 deti(5,5,4)=-cg1*ta1 deti(4,5,4)=-sg1/ca1/ca1 deti(5,6,4)=cg1/ca1 deti(4,6,4)=sg1*sa1/ca1/ca1 deti(5,4,6)=cg1 deti(5,5,6)=sg1*ta1 deti(4,5,6)=-cg1/ca1/ca1 deti(5,6,6)=-sg1/ca1 deti(4,6,6)=cg1*sa1/ca1/ca1 deti(12,11,11)=sg2 deti(12,12,11)=-cg2*ta2 deti(11,12,11)=-sg2/ca2/ca2 deti(12,13,11)=cg2/ca2 deti(11,13,11)=sg2*sa2/ca2/ca2 deti(12,11,13)=cg2 deti(12,12,13)=sg2*ta2 deti(11,12,13)=-cg2/ca2/ca2 deti(12,13,13)=-sg2/ca2 deti(11,13,13)=cg2*sa2/ca2/ca2 c c da(1,1,1)=sg1*t(2,1) da(2,1,1)=ta1*t(2,1) da(3,1,1)=cg1*t(2,1) da(4,1,1)=t(2,1)*(ru*sg1*ta1-rv+rw*cg1*ta1) da(5,1,1)=(-ru*cg1+rw*sg1)*t(2,1) da(6,1,1)=(-ru*sg1-rv*ta1-rw*cg1)*t(2,2) da(8,1,1)=t(2,1)*(-sg1) da(9,1,1)=t(2,1)*(-ta1) da(10,1,1)=t(2,1)*(-cg1) c da(1,1,2)=-t(3,1) da(3,1,2)=t(1,1) da(4,1,2)=-t(2,1)*ru*cg1+t(2,1)*rw*sg1-cp1*sa1 da(5,1,2)=-t(1,1)*ru-t(3,1)*rw+ca1*sa1*sp1 da(6,1,2)=-t(1,2)*rw+t(3,2)*ru da(8,1,2)=t(3,1) da(10,1,2)=-t(1,1) c da(1,1,3)=-t(1,2) da(2,1,3)=-t(2,2) da(3,1,3)=-t(3,2) da(4,1,3)=-t(2,2)*(sg1*ru+ta1*rv+cg1*rw)+t(1,2)*sa1*sg1 & -t(2,2)*ca1+t(3,2)*sa1*cg1 da(5,1,3)=t(3,2)*ru-t(1,2)*rw-t(1,2)*ca1*cg1+t(3,2)*ca1*sg1 da(6,1,3)=-t(1,1)*ru-t(2,1)*rv-t(3,1)*rw c c da(4,1,3)=-ru*t(2,1)*sg1+t(1,1)*sa1*sg1-rv*t(2,2)*ta1 c & -t(2,2)*ca1-rw*t(2,2)*cg1+t(3,2)*sa1*cg1 c da(5,1,3)=ru*t(3,1)-t(1,1)*ca1*cg1-rw*t(1,2)+t(3,2)*ca1*sg1 c da(6,1,3)=ru*t(1,2)-rv*t(2,1)-rw*t(3,1) c da(8,1,3)=t(1,2) da(9,1,3)=t(2,2) da(10,1,3)=t(3,2) c da(1,2,1)=t(2,2)*sg1 da(2,2,1)=t(2,2)*ta1 da(3,2,1)=t(2,2)*cg1 da(4,2,1)=t(2,2)*(ru*sg1*ta1-rv+rw*ta1*cg1) da(5,2,1)=t(2,2)*(-ru*cg1+rw*sg1) da(6,2,1)=t(2,1)*(ru*sg1+rv*ta1+rw*cg1) da(8,2,1)=-t(2,2)*sg1 da(9,2,1)=-t(2,2)*ta1 da(10,2,1)=-t(2,2)*cg1 c da(1,2,2)=-t(3,2) da(3,2,2)=t(1,2) da(4,2,2)=t(2,2)*(-ru*cg1+rw*sg1)+t(3,2)*sa1*sg1-t(1,2)*sa1*cg1 da(5,2,2)=-t(1,2)*ru-t(3,2)*rw-t(3,2)*ca1*cg1-t(1,2)*ca1*sg1 da(6,2,2)=-t(3,1)*ru+t(1,1)*rw c c da(4,2,2)=sa1*sp1-t(2,2)*cg1+t(1,3)*cp1 c da(5,2,2)=-ca1*sa1*cp1-t(2,2)*cg1+t(1,3)*cp1 c da(6,2,2)=-t(3,1)*ru-t(1,1)*rw c da(8,2,2)=t(3,2) da(10,2,2)=-t(1,2) c do 50 i=1,3 da(i,2,3)=t(i,1) da(7+i,2,3)=-t(i,1) 50 continue da(4,2,3)=t(2,1)*ru*sg1+t(2,1)*rv*ta1+t(2,1)*rw*cg1 & -t(1,1)*sa1*sg1+t(2,1)*ca1-t(3,1)*sa1*cg1 da(5,2,3)=ca1*cp1+t(1,1)*rw-t(3,1)*ru da(6,2,3)=-t(1,2)*ru-t(2,2)*rv-t(3,2)*rw c da(1,3,1)=sa1*sg1 da(2,3,1)=-ca1 da(3,3,1)=sa1*cg1 da(4,3,1)=-ca1*sg1*ru-sa1*rv-ca1*cg1*rw-1 c c da(4,3,1)=sa1**2*sg1**2-ca1**2+sa1**2*cg1**2-ru*sg1*ca1-rv*sa1 c & -rw*cg1*ca1 c da(5,3,1)=-ru*t(2,3)*cg1+rw*sg1*t(2,3) da(8,3,1)=-t(2,3)*sg1 da(9,3,1)=ca1 da(10,3,1)=-t(2,3)*cg1 c da(1,3,2)=-ca1*cg1 da(3,3,2)=t(1,3) da(4,3,2)=-sa1*(ru*cg1-rw*sg1) da(5,3,2)=-ca1**2-ca1*(ru*sg1+rw*cg1) da(8,3,2)=ca1*cg1 da(10,3,2)=-t(1,3) c da(4,4,1)=-cgg*sg1*t(2,1)+phh*sgg*saa*sg1*t(2,1) & -phh*caa*ta1*t(2,1)-ry*t(2,1)/ca1/ca1+sgg*cg1*t(2,1) & +phh*cgg*saa*cg1*t(2,1)+(rx*sg1+ry*ta1+rz*cg1)*ta1*t(2,1) da(11,4,1)=cgg*sg1*t(2,1)+phh*sgg*saa*sg1*t(2,1) & -phh*caa*ta1*t(2,1)-sgg*cg1*t(2,1)+phh*cgg*saa*cg1*t(2,1) da(5,4,1)=-allp*sgg*sg1*t(2,1)-phh*cgg*caa*sg1*t(2,1) & -rx*cg1*t(2,1)+ta1*t(2,1)-allp*cgg*cg1*t(2,1) & +phh*sgg*caa*cg1*t(2,1)+rz*sg1*t(2,1) da(12,4,1)=-allp*sgg*sg1*t(2,1)-phh*cgg*caa*sg1*t(2,1) & -ta1*t(2,1)-allp*cgg*cg1*t(2,1)+phh*sgg*caa*cg1*t(2,1) da(6,4,1)=sgg*caa*sg1*t(2,1)+saa*ta1*t(2,1)+cgg*caa*cg1*t(2,1) & +(-rx*sg1-ry*ta1-rz*cg1)*t(2,2) da(13,4,1)=-sgg*caa*sg1*t(2,1)-saa*ta1*t(2,1)-cgg*caa*cg1*t(2,1) c da(4,4,2)=-t(2,1)*rx*cg1+t(2,1)*rz*sg1+cgg*t(3,1) & -phh*sgg*saa*t(3,1)+phh*cgg*saa*t(1,1)+sgg*t(1,1) da(11,4,2)=-cgg*t(3,1)-phh*sgg*saa*t(3,1)-sgg*t(1,1) & +phh*cgg*saa*t(1,1) da(5,4,2)=-t(1,1)*rx-t(3,1)*rz+allp*sgg*t(3,1) & +phh*cgg*caa*t(3,1)-allp*cgg*t(1,1)+phh*sgg*caa*t(1,1) da(12,4,2)=allp*sgg*t(3,1)+phh*cgg*caa*t(3,1)-allp*cgg*t(1,1) & +phh*sgg*caa*t(1,1) da(6,4,2)=t(3,2)*rx-t(1,2)*rz-sgg*caa*t(3,1)+cgg*caa*t(1,1) da(13,4,2)=sgg*caa*t(3,1)-cgg*caa*t(1,1) c da(4,4,3)=-t(2,2)*rx*sg1-t(2,2)*ry*ta1-t(2,2)*rz*cg1+cgg*t(1,2) & -phh*sgg*saa*t(1,2)+phh*caa*t(2,2)-phh*cgg*saa*t(3,2) & -sgg*t(3,2) da(11,4,3)=-cgg*t(1,2)-phh*sgg*saa*t(1,2)+phh*caa*t(2,2) & +sgg*t(3,2)-phh*cgg*saa*t(3,2) da(5,4,3)=t(3,2)*rx-t(1,2)*rz+allp*sgg*t(1,2)+phh*cgg*caa*t(1,2) & -t(2,2)+allp*cgg*t(3,2)-phh*sgg*caa*t(3,2) da(12,4,3)=allp*sgg*t(1,2)+phh*cgg*caa*t(1,2)+t(2,2) & +allp*cgg*t(3,2)-phh*sgg*caa*t(3,2) da(6,4,3)=-t(1,1)*rx-t(2,1)*ry-t(3,1)*rz-sgg*caa*t(1,2) & -saa*t(2,2)-cgg*caa*t(3,2) da(13,4,3)=sgg*caa*t(1,2)+saa*t(2,2)+cgg*caa*t(3,2) c da(4,5,1)=-cgg*sg1*t(2,2)+phh*sgg*saa*sg1*t(2,2) & -phh*caa*ta1*t(2,2)-ry*t(2,2)/ca1/ca1-sgg*cg1*t(2,2) & -phh*cgg*saa*cg1*t(2,2)+ta1*t(2,2)*(rx*sg1+ry*ta1+rz*cg1) da(11,5,1)=cgg*sg1*t(2,2)+phh*sgg*saa*sg1*t(2,2) & -phh*caa*ta1*t(2,2)-sgg*cg1*t(2,2)+phh*cgg*saa*cg1*t(2,2) da(5,5,1)=-allp*sgg*sg1*t(2,2)-phh*cgg*caa*sg1*t(2,2) & -rx*t(2,2)*cg1-allp*cgg*cg1*t(2,2)+ta1*t(2,2) & +phh*sgg*caa*cg1*t(2,2)+rz*t(2,2)*sg1 da(12,5,1)=-allp*sgg*sg1*t(2,2)-phh*cgg*caa*sg1*t(2,2) & -ta1*t(2,2)-allp*cgg*cg1*t(2,2)+phh*sgg*caa*cg1*t(2,2) da(6,5,1)=sgg*caa*sg1*t(2,2)+saa*ta1*t(2,2)+cgg*caa*cg1*t(2,2) & +(rx*sg1+ry*ta1+rz*cg1)*t(2,1) da(13,5,1)=-sgg*caa*sg1*t(2,2)-saa*ta1*t(2,2)-cgg*caa*cg1*t(2,2) c da(4,5,2)=-cg1*rx*t(2,2)+t(2,2)*rz*sg1+cgg*t(3,2) & -phh*sgg*saa*t(3,2)+sgg*t(1,2)+phh*cgg*saa*t(1,2) da(11,5,2)=-cgg*t(3,2)-phh*sgg*saa*t(3,2)-sgg*t(1,2) & +phh*cgg*saa*t(1,2) da(5,5,2)=-t(1,2)*rx-t(3,2)*rz+allp*sgg*t(3,2) & +phh*cgg*caa*t(3,2)-allp*cgg*t(1,2)+phh*sgg*caa*t(1,2) da(12,5,2)=+allp*sgg*t(3,2)+phh*cgg*caa*t(3,2)-allp*cgg*t(1,2) & +phh*sgg*caa*t(1,2) da(6,5,2)=-t(3,1)*rx+t(1,1)*rz-sgg*caa*t(3,2)+cgg*caa*t(1,2) da(13,5,2)=sgg*caa*t(3,2)-cgg*caa*t(1,2) c da(4,5,3)=t(2,1)*rx*sg1+t(2,1)*ry*ta1+t(2,1)*rz*cg1-cgg*t(1,1) & +phh*sgg*saa*t(1,1)-phh*caa*t(2,1)+sgg*t(3,1)+phh*cgg*saa*t(3,1) da(11,5,3)=cgg*t(1,1)+phh*sgg*saa*t(1,1)-phh*caa*t(2,1) & -sgg*t(3,1)+phh*cgg*saa*t(3,1) da(5,5,3)=-t(3,1)*rx+t(1,1)*rz-allp*sgg*t(1,1) & -phh*cgg*caa*t(1,1)+t(2,1)-allp*cgg*t(3,1)+phh*sgg*caa*t(3,1) da(12,5,3)=-allp*sgg*t(1,1)-phh*cgg*caa*t(1,1)-t(2,1) & -allp*cgg*t(3,1)+phh*sgg*caa*t(3,1) da(6,5,3)=-t(1,2)*rx-t(2,2)*ry-t(3,2)*rz+sgg*caa*t(1,1) & +saa*t(2,1)+cgg*caa*t(3,1) da(13,5,3)=-sgg*caa*t(1,1)-saa*t(2,1)-cgg*caa*t(3,1) c c da(4,6,1)=cgg*sg1*t(2,3)+phh*sgg*saa*sg1*t(2,3) c & +phh*caa*ca1-ry*sa1+sgg*cg1*t(2,3) c & +phh*cgg*saa*cg1*t(2,3)+ca1*(-rx*sg1-rz*cg1) c da(11,6,1)=-cgg*sg1*t(2,3)+phh*sgg*saa*sg1*t(2,3)-sgg*cg1*t(2,3) c & +phh*caa*ca1+phh*cgg*saa*cg1*t(2,3) c da(5,6,1)=allp*sgg*sg1*t(2,3)-phh*caa*cgg*sg1*t(2,3) c & -rx*t(2,3)*cg1-ca1-allp*cgg*cg1*t(2,3) c & +phh*sgg*caa*cg1*t(2,3) c da(12,6,1)=allp*sgg*sg1*t(2,3)-phh*cgg*caa*sg1*t(2,3) c & +ca1-allp*cgg*cg1*t(2,3)+phh*sgg*caa*cg1*t(2,3) c da(6,6,1)=sgg*caa*sg1*t(2,3)-saa*ca1 c & +cgg*caa*cg1*t(2,3) c da(13,6,1)=-sgg*caa*sg1*t(2,3)+saa*ca1 c & -cgg*caa*cg1*t(2,3) c da(4,6,1)=-ca1*sg1*rx-sa1*ry-ca1*cg1*rw-cgg*sa1*sg1 & +phh*sgg*saa*sa1*sg1+phh*caa*ca1+sgg*sa1*cg1 & +phh*cgg*saa*sa1*cg1 da(11,6,1)=cgg*sa1*sg1+phh*sgg*saa*sa1*sg1+phh*caa*ca1 & -sgg*sa1*cg1+phh*cgg*saa*sa1*cg1 da(5,6,1)=-sa1*cg1*rx+sa1*sg1*rw-allp*sgg*sa1*sg1 & -phh*cgg*caa*sa1*sg1-ca1-allp*cgg*sa1*cg1+phh*saa*caa*sa1*cg1 da(12,6,1)=-allp*sgg*sa1*sg1-phh*cgg*caa*sa1*sg1 & +ca1-allp*cgg*sa1*cg1+phh*saa*caa*sa1*cg1 da(6,6,1)=sgg*caa*sa1*sg1-saa*ca1+cgg*caa*sa1*cg1 da(13,6,1)=-da(6,6,1) c da(4,6,2)=-rx*sa1*cg1+rz*sa1*sg1+cgg*cg1*ca1 & -phh*sgg*saa*cg1*ca1+sgg*t(1,3)+phh*cgg*saa*t(1,3) da(11,6,2)=-cgg*cg1*ca1-phh*sgg*saa*cg1*ca1 & -sgg*t(1,3)+phh*cgg*saa*t(1,3) da(5,6,2)=-rz*ca1*cg1-rx*ca1*sg1 & +allp*sgg*cg1*ca1+phh*cgg*caa*cg1*ca1 & -allp*cgg*t(1,3)+phh*sgg*caa*t(1,3) da(12,6,2)=allp*sgg*cg1*ca1+phh*cgg*caa*cg1*ca1 & -allp*cgg*t(1,3)+phh*sgg*caa*t(1,3) da(6,6,2)=-sgg*caa*cg1*ca1+cgg*caa*t(1,3) da(13,6,2)=sgg*caa*cg1*ca1-cgg*caa*t(1,3) c c da(4,5,1)= & -((-p1 + p2)*Cos((a1 + a2)/2)*Cos(p1)*Sin(a1))/2 - & Cos(a1)*Cos(p1)*(-g1 + g2 + (-p1 + p2)*Sin((a1 + a2)/2)) - & Cos(a1)*Cos(g1)*Cos(p1)*(-((-p1 + p2)*Cos((g1 + g2)/2)* & Sin((a1 + a2)/2))/2 - Sin((g1 + g2)/2)) + & Cos(g1)*Cos(p1)*Sin(a1)*((-p1 + p2)*Cos((a1 + a2)/2)* & Cos((g1 + g2)/2) + (-a1 + a2)*Sin((g1 + g2)/2)) + & Cos(p1)*Sin(a1)*Sin(g1)*(-((-a1 + a2)*Cos((g1 + g2)/2)) + & (-p1 + p2)*Cos((a1 + a2)/2)*Sin((g1 + g2)/2)) - & Cos(a1)*Cos(p1)*Sin(g1)*(Cos((g1 + g2)/2) - & (-p1 + p2)*Sin((a1 + a2)/2)*Sin((g1 + g2)/2)/2) c da(4,6,1)= & (-p1 + p2)*Cos(a1)*Cos((a1 + a2)/2)/2 - & Sin(a1)*(-g1 + g2 + (-p1 + p2)*Sin((a1 + a2)/2)) - & Cos(g1)*Sin(a1)*(-((-p1 + p2)*Cos((g1 + g2)/2)* & Sin((a1 + a2)/2))/2 - Sin((g1 + g2)/2)) - & Cos(a1)*Cos(g1)*((-p1 + p2)*Cos((a1 + a2)/2)* & Cos((g1 + g2)/2) + (-a1 + a2)*Sin((g1 + g2)/2)) - & Cos(a1)*Sin(g1)*(-((-a1 + a2)*Cos((g1 + g2)/2)) + & (-p1 + p2)*Cos((a1 + a2)/2)*Sin((g1 + g2)/2)) - & Sin(a1)*Sin(g1)*(Cos((g1 + g2)/2) - & (-p1 + p2)*Sin((a1 + a2)/2)*Sin((g1 + g2)/2)/2) c da(5,6,1)= & -Cos(a1) - Sin(a1)*Sin(g1)*((-p1 + p2)*Cos((a1 + a2)/2)* & Cos((g1 + g2)/2)/2 + (-a1 + a2)*Sin((g1 + g2)/2)/2) + & Sin(a1)*Sin(g1)*((-p1 + p2)*Cos((a1 + a2)/2)* & Cos((g1 + g2)/2) + (-a1 + a2)*Sin((g1 + g2)/2)) - & Cos(g1)*Sin(a1)*((-a1 + a2)*Cos((g1 + g2)/2)/2 - & (-p1 + p2)*Cos((a1 + a2)/2)*Sin((g1 + g2)/2)/2) - & Cos(g1)*Sin(a1)*(-((-a1 + a2)*Cos((g1 + g2)/2)) + & (-p1 + p2)*Cos((a1 + a2)/2)*Sin((g1 + g2)/2)) c c a(1,1)=(-ru*sg1-rv*ta1-rw*cg1)*t(2,1) a(1,2)=t(3,1)*ru-t(1,1)*rw a(1,3)=t(1,2)*ru+t(2,2)*rv+t(3,2)*rw a(2,1)=t(2,2)*(-ru*sg1-rv*ta1-rw*cg1) a(2,2)=t(3,2)*ru-t(1,2)*rw a(2,3)=-t(1,1)*ru-t(2,1)*rv-t(3,1)*rw a(3,1)=(-ru*sg1*sa1+rv*ca1-rw*cg1*sa1) a(3,2)=ca1*ru*cg1-t(1,3)*rw a(3,3)=0.d0 a(4,1)=t(2,1)*(-rx*sg1-ry*ta1-rz*cg1) a(4,2)=t(3,1)*rx-t(1,1)*rz a(4,3)=t(1,2)*rx+t(2,2)*ry+t(3,2)*rz a(5,1)=t(2,2)*(-rx*sg1-ry*ta1-rz*cg1) a(5,2)=t(3,2)*rx-t(1,2)*rz a(5,3)=-t(1,1)*rx-t(2,1)*ry-t(3,1)*rz a(6,1)=(-rx*sg1*sa1+ry*ca1-rz*cg1*sa1) a(6,2)=ca1*rx*cg1-t(1,3)*rz a(6,3)=0.d0 c do 100 i=1,3 do 100 j=1,3 r7x14(i,j)=-t(j,i) r7x14(i,7+j)=t(j,i) 100 continue do 110 i=1,3 r7x14(3+i,4)=a(3+i,1)+t(1,i)*cf1(4,4)+t(2,i)*cf1(5,4) & +t(3,i)*cf1(6,4) r7x14(i,4)= & a(i,1)+t(1,i)*cf1(1,4)+t(2,i)*cf1(2,4)+t(3,i)*cf1(3,4) r7x14(i,5)=a(i,2)+t(1,i)*cf1(1,5)+t(3,i)*cf1(3,5) r7x14(3+i,5)=a(3+i,2)+t(1,i)*cf1(4,5)-t(2,i)+t(3,i)*cf1(6,5) r7x14(i,6)=a(i,3) r7x14(3+i,6)= & a(3+i,3)+t(1,i)*cf1(4,6)+t(2,i)*cf1(5,6)+t(3,i)*cf1(6,6) r7x14(3+i,11)=t(1,i)*cf2(4,4)+t(2,i)*cf1(5,4)+t(3,i)*cf2(6,4) r7x14(3+i,12)=t(1,i)*cf2(4,5)+t(2,i)+t(3,i)*cf1(6,5) r7x14(3+i,13)=t(1,i)*cf2(4,6)-t(2,i)*cf1(5,6)-t(3,i)*cf1(6,6) 110 continue r7x14(7,14)=1.d0 c do 433 k=1,14 do 433 i=1,7 do 433 l=1,7 tdb1(k,i,l)=0.d0 tdb2(k,i,l)=0.d0 dtb1(k,i,l)=0.d0 dtb2(k,i,l)=0.d0 do 433 j=1,7 tdb1(k,i,l)=tdb1(k,i,l)+t1(j,i)*dcf1(k,j,l) tdb2(k,i,l)=tdb2(k,i,l)+t1(j,i)*dcf2(k,j,l) dtb1(k,i,l)=dtb1(k,i,l)+dt1(k,j,i)*cf1(j,l) dtb2(k,i,l)=dtb2(k,i,l)+dt1(k,j,i)*cf2(j,l) 433 continue do 444 k=1,14 do 444 i=1,7 do 444 j=1,7 dr7x14(k,i,j)=tdb1(k,i,j)+dtb1(k,i,j) dr7x14(k,i,7+j)=tdb2(k,i,j)+dtb2(k,i,j) 444 continue do 445 k=1,14 do 445 i=1,6 do 445 j=1,3 dr7x14(k,i,j+3)=da(k,i,j)+dr7x14(k,i,j+3) 445 continue c return end c c c 出力した接線行列の固有値解析をせずに、最後のN-R反復で求まる c (微少な)変位増分を固有ベクトル(座屈モード)の近似として c 利用する場合のために、近似モードを出力する。 subroutine modeout(nofe,dd) implicit real*8(a-h,o-z) dimension dd(7*256+8) open(9,file='modeu.d') call modeo7(1,nofe,dd) close(9) open(9,file='modev.d') call modeo7(2,nofe,dd) close(9) open(9,file='modew.d') call modeo7(3,nofe,dd) close(9) open(9,file='modex.d') call modeo7(4,nofe,dd) close(9) open(9,file='modey.d') call modeo7(5,nofe,dd) close(9) open(9,file='modez.d') call modeo7(6,nofe,dd) close(9) open(9,file='modep.d') call modeo7(7,nofe,dd) close(9) return end c subroutine modeo7(i7,nofe,dd) implicit real*8(a-h,o-z) dimension dd(7*256+8) do 10 i=1,nofe-1 dmode=dd(i*7+i7) write(9,*) i+1, dmode 10 continue return end c c 変位dに外乱(座屈モードなど)を読み込む c subroutine modein(nofe,d,gmode) implicit real*8(a-h,o-z) dimension d(14*256) open(9,file='modeu.d') call modei7(1,nofe,d,gmode) close(9) open(9,file='modev.d') call modei7(2,nofe,d,gmode) close(9) open(9,file='modew.d') call modei7(3,nofe,d,gmode) close(9) open(9,file='modex.d') call modei7(4,nofe,d,gmode) close(9) open(9,file='modey.d') call modei7(5,nofe,d,gmode) close(9) open(9,file='modez.d') call modei7(6,nofe,d,gmode) close(9) open(9,file='modep.d') call modei7(7,nofe,d,gmode) close(9) return end c subroutine modei7(i7,nofe,d,gmode) implicit real*8(a-h,o-z) dimension d(14*256) read(9,*) idummy,dmode d(i*14+i7)=d(i*14+i7)+gmode*dmode do 10 i=0,nofe-2 i1=i+1 read(9,*) idummy,dmode d(i*14+7+i7)=d(i*14+7+i7)+gmode*dmode d(i1*14+i7)=d(i1*14+i7)+gmode*dmode 10 continue read(9,*) idummy,dmode n1=nofe-1 d(n1*14+7+i7)=d(n1*14+7+i7)+gmode*dmode return end c c