c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c+ + c+ 4-node plate bending element program + c+ by T. Nishihara + c+ + c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ program kani implicit real*8 (a-h,o-z) c parameter( nele = 5, node = 12, ndof = 6 * node ) parameter( nbcc = 2, nout = 1 ) c dimension bc(ndof), forc(ndof),dd(ndof+1), d(ndof,24) dimension s(24,24), sk(24,24), sinc(ndof+1,ndof+1) dimension ddl(ndof), xb(ndof+1), d24(24), dd24(24) dimension ba24(24), ba(ndof+1), alp3(3), cor(ndof,2) dimension alp(ndof), snr2(ndof,ndof) dimension jcon(4), nodn(ndof,4), mds(nout), mmds(nout) dimension nbc(nbcc) c character*20 file1 c open(7,file='input.dat') open(8,file='outkani.dat') c c******************************************************************** c データの読み込み c******************************************************************** c call clrv (ndof,alp) c read(7,*) th,e,um c do 11 i=1,node read(7,*) cor(i,1), cor(i,2) 11 continue c do 13 i=1,nele read(7,*) nodn(i,1), nodn(i,2), nodn(i,3), nodn(i,4) 13 continue c do 50 i=1,ndof bc(i) = 1.d0 50 continue c read(7,*) nxx, nyy, nzz, nrx, nry read(7,*) (nbc(i), i=1, nxx) do 130 i=1,nxx bc ( nbc(i) * 6 - 5 ) = 0.d0 130 continue read(7,*) (nbc(i), i=1, nyy) do 131 i=1,nyy bc ( nbc(i) * 6 - 4 ) = 0.d0 131 continue read(7,*) (nbc(i), i=1, nzz) do 132 i=1,nzz bc ( nbc(i) * 6 - 3 ) = 0.d0 132 continue read(7,*) (nbc(i), i=1, nrx) do 133 i=1,nrx bc ( nbc(i) * 6 - 2 ) = 0.d0 133 continue read(7,*) (nbc(i), i=1, nry) do 134 i=1,nry bc ( nbc(i) * 6 - 1 ) = 0.d0 134 continue do 135 i=1,node bc ( i * 6 ) = 0.d0 135 continue c call clrv (ndof,forc) read(7,*) nf do 60 i=1,nf read(7,*) nodd, xfor, yfor, zfor, xrot, yrot, zrot forc( nodd * 6 - 5 ) = xfor forc( nodd * 6 - 4 ) = yfor forc( nodd * 6 - 3 ) = zfor forc( nodd * 6 - 2 ) = xrot forc( nodd * 6 - 1 ) = yrot forc( nodd * 6 ) = zrot 60 continue c read(7,*) np do 560 i=1,np read(7,561) mds(i), mmds(i), file1 open(20+i,file=file1) 560 continue c 後藤追記(04/1/19) c Linux 上の g77 では大丈夫なようだけど、Cygwin の g77 とかでは、 c 以下のフォーマットを入れないとちゃんと読めないようだ。 c だから入力データもここのフォーマットと合わせて書く必要がある 561 format(2i3, a6) c read(7,*) ma if ( ma.eq.0 ) go to 695 do 570 ie=1,nele read(7,*) alp(ie*3-2),alp(ie*3-1),alp(ie*3) 570 continue c 695 continue c write(8,*)'====================================================' write(8,*)' output of program' write(8,*)'====================================================' write(8,*) write(8,*) write(8,*)'----------------------------------------------------' write(8,*)' parameter' write(8,*)'----------------------------------------------------' write(8,2005) nele, node write(8,2001) e, um write(8,2002) th c 2001 format('Young modulas = ',1pe15.7,5x,'Poisson ration = ', & 1pe15.7) 2002 format('Thickness = ',1pe15.7) 2005 format('Number of element:',i5,5x,'Number of node:',i5) c write(8,*) write(8,*)'----------------------------------------------------' write(8,*)' co-ordinate of nodes' write(8,*)'----------------------------------------------------' write(8,2014) do 12 i=1,node write(8,2015) i,cor(i,1), cor(i,2) 12 continue 2014 format(2x,'node number',12x,'x',20x,'y') 2015 format('node = ',i5,5x,1pe15.7,5x,1pe15.7) c write(8,*) write(8,*)'----------------------------------------------------' write(8,*)' element' write(8,*)'----------------------------------------------------' do 40 i=1,nele write(8,2007) i,nodn(i,1),nodn(i,2),nodn(i,3),nodn(i,4) 40 continue 2007 format('element = ',i5,5x,i5,i5,i5,i5) c write(8,*) write(8,*)'----------------------------------------------------' write(8,*)' boundary condition' write(8,*)'----------------------------------------------------' write(8,2009) do 2008 i=1,node write(8,2010) i,bc(6*i-5),bc(6*i-4),bc(6*i-3), & bc(6*i-2),bc(6*i-1),bc(6*i) 2008 continue 2009 format(3x,'node number',8x,'x',9x,'y',9x,'z',7x, & 'rot-x',5x,'rot-y',5x,'rot-z') 2010 format('node = ',i5,1x, & 1pe10.1,1pe10.1,1pe10.1,1pe10.1,1pe10.1,1pe10.1) c write(8,*) write(8,*)'----------------------------------------------------' write(8,*)' nodal force' write(8,*)'----------------------------------------------------' write(8,2009) do 2011 i=1,node write(8,2010) i,forc(6*i-5),forc(6*i-4),forc(6*i-3), & forc(6*i-2),forc(6*i-1),forc(6*i) 2011 continue c if ( ma.eq.0 ) go to 696 write(8,*) write(8,*)'----------------------------------------------------' write(8,*)' initial angle of element' write(8,*)'----------------------------------------------------' write(8,2017) do 571 ie=1,nele write(8,2016) ie,alp(ie*3-2),alp(ie*3-1),alp(ie*3) 571 continue 2017 format(2x,'element number',10x,'alpha',15x,'beta',15x,'gamma') 2016 format('element = ',i5,5x,1pe15.7,5x,1pe15.7,5x,1pe15.7) 696 continue c write(*,*)'start loding step' c c******************************************************************** c 第1増分ステップの第1収斂ステップ c******************************************************************** read(7,*) r write(8,*) write(8,*) write(8,2012) r write(8,*) read(7,*) rn write(8,2013) rn 2012 format('first arc = ',1pe15.7) 2013 format('next arc = ',1pe15.7) c do 225 i=1,ndof do 225 j=1,24 d(i,j) = 0.d0 225 continue call clrv (ndof+1,dd) call clrv (ndof,ddl) call clrv (ndof+1,ba) dp = 0.d0 p = 0.d0 dpl = 0.d0 r0 = 0.d0 write(3) dd,d,ddl,dp,dpl,r0,p rewind 3 r0 = r dpl = r nacl = 0 nbcl = 0 nccl = 0 konv = 0 c c******************************************************************** c 弧長増分過程 c******************************************************************** c 1000 continue nofcs = 0 ninc = ninc + 1 c if(ninc.gt.1000) then c write(12,*) c write(12,*)'収束しないよ!' c close(12) c close(13) c stop c endif c call clrmx (ndof+1,sinc) c do 2000 ie=1,nele do 440 i = 1,4 jcon(i) = nodn(ie, i) 440 continue ib1 = (jcon(1) - 1) * 6 ib2 = (jcon(2) - 1) * 6 ib3 = (jcon(3) - 1) * 6 ib4 = (jcon(4) - 1) * 6 c a = cor( jcon(2), 1 ) - cor( jcon(1), 1 ) b = cor( jcon(4), 2 ) - cor( jcon(1), 2 ) call smx(s,th,e,um,a,b) c do 800 i=1,3 alp3(i) = alp((ie-1)*3+i) 800 continue c do 150 i=1,6 dd24(i) = dd(ib1+i) dd24(i+6) = dd(ib2+i) dd24(i+12) = dd(ib3+i) dd24(i+18) = dd(ib4+i) d24(i) = d(ie,i) d24(i+6) = d(ie,i+6) d24(i+12) = d(ie,i+12) d24(i+18) = d(ie,i+18) 150 continue c call stmx (s,sk,a,b,d24,dd24,ba24,alp3) c do 160 i=1,6 d(ie,i) = d24(i) d(ie,i+6) = d24(i+6) d(ie,i+12) = d24(i+12) d(ie,i+18) = d24(i+18) 160 continue c do 170 i=1,6 do 170 j=1,6 sinc(ib1+i,ib1+j) = sinc(ib1+i,ib1+j) + sk(i,j) sinc(ib1+i,ib2+j) = sinc(ib1+i,ib2+j) + sk(i,j+6) sinc(ib1+i,ib3+j) = sinc(ib1+i,ib3+j) + sk(i,j+12) sinc(ib1+i,ib4+j) = sinc(ib1+i,ib4+j) + sk(i,j+18) sinc(ib2+i,ib1+j) = sinc(ib2+i,ib1+j) + sk(i+6,j) sinc(ib2+i,ib2+j) = sinc(ib2+i,ib2+j) + sk(i+6,j+6) sinc(ib2+i,ib3+j) = sinc(ib2+i,ib3+j) + sk(i+6,j+12) sinc(ib2+i,ib4+j) = sinc(ib2+i,ib4+j) + sk(i+6,j+18) sinc(ib3+i,ib1+j) = sinc(ib3+i,ib1+j) + sk(i+12,j) sinc(ib3+i,ib2+j) = sinc(ib3+i,ib2+j) + sk(i+12,j+6) sinc(ib3+i,ib3+j) = sinc(ib3+i,ib3+j) + sk(i+12,j+12) sinc(ib3+i,ib4+j) = sinc(ib3+i,ib4+j) + sk(i+12,j+18) sinc(ib4+i,ib1+j) = sinc(ib4+i,ib1+j) + sk(i+18,j) sinc(ib4+i,ib2+j) = sinc(ib4+i,ib2+j) + sk(i+18,j+6) sinc(ib4+i,ib3+j) = sinc(ib4+i,ib3+j) + sk(i+18,j+12) sinc(ib4+i,ib4+j) = sinc(ib4+i,ib4+j) + sk(i+18,j+18) 170 continue 2000 continue c do 180 i=1,ndof sinc(ndof+1,i) = ddl(i) / r0 c ddl(i) = 0.d0 sinc(i,ndof+1) = -forc(i) 180 continue sinc(ndof+1,ndof+1) = dpl / r0 c call clrv (ndof+1,xb) xb(ndof+1) = r r0 = r dpl = 0.d0 c do 190 i=1,ndof do 190 j=1,ndof sinc(i,j) = bc(i) * sinc(i,j) sinc(j,i) = bc(i) * sinc(j,i) 190 continue do 200 i=1,ndof if(sinc(i,i) .ne. 0.d0) go to 200 sinc(i,i) = 1.d0 200 continue c call clrv (ndof+1,dd) call solve (sinc,dd,xb,dp,det,idt,ndof,ndof+1,jdt) c dp = dd(ndof+1) c c******************************************************************** c Newton-Raphson収斂過程 c******************************************************************** 1100 nofcs = nofcs + 1 c do 210 i=1,ndof ddl(i) = ddl(i) + dd(i) ddi2 = ddi2 + dd(i)*dd(i) ddl2 = ddl2 + ddl(i)*ddl(i) 210 continue c dpl = dpl + dp p = p + dp c if(ddl2.eq.0) go to 1200 arc = ddi2 / ddl2 ddi2 = 0.d0 ddl2 = 0.d0 if ( nofcs.gt.6 ) then c write(*,*)'6ステップで収斂しないので、弧長を半分にする。' nlet = nlet + 1 konv = konv + 1 if( nlet.gt.10 ) then p1 = p go to 4020 end if go to 4000 end if write(*,3020) arc 3020 format(5x,'arc = ',1pe15.7) c if(arc.lt.1.d-10) go to 1300 c 1200 call clrmx (ndof+1,sinc) call clrv (ndof+1,ba) c do 220 i=1,ndof ba(i) = p * forc(i) 220 continue c do 2100 ie=1,nele do 240 i = 1,4 jcon(i) = nodn(ie, i) 240 continue ib1 = (jcon(1) - 1) * 6 ib2 = (jcon(2) - 1) * 6 ib3 = (jcon(3) - 1) * 6 ib4 = (jcon(4) - 1) * 6 cc a = cor( jcon(2), 1 ) - cor( jcon(1), 1 ) b = cor( jcon(4), 2 ) - cor( jcon(1), 2 ) call smx(s,th,e,um,a,b) c do 810 i=1,3 alp3(i) = alp((ie-1)*3+i) 810 continue c do 250 i=1,6 dd24(i) = dd(ib1+i) dd24(i+6) = dd(ib2+i) dd24(i+12) = dd(ib3+i) dd24(i+18) = dd(ib4+i) d24(i) = d(ie,i) d24(i+6) = d(ie,i+6) d24(i+12) = d(ie,i+12) d24(i+18) = d(ie,i+18) ba24(i) = ba(ib1+i) ba24(i+6) = ba(ib2+i) ba24(i+12) = ba(ib3+i) ba24(i+18) = ba(ib4+i) 250 continue c call stmx (s,sk,a,b,d24,dd24,ba24,alp3) c do 260 i=1,6 d(ie,i) = d24(i) d(ie,i+6) = d24(i+6) d(ie,i+12) = d24(i+12) d(ie,i+18) = d24(i+18) ba(ib1+i) = ba24(i) ba(ib2+i) = ba24(i+6) ba(ib3+i) = ba24(i+12) ba(ib4+i) = ba24(i+18) 260 continue c do 270 i=1,6 do 270 j=1,6 sinc(ib1+i,ib1+j) = sinc(ib1+i,ib1+j) + sk(i,j) sinc(ib1+i,ib2+j) = sinc(ib1+i,ib2+j) + sk(i,j+6) sinc(ib1+i,ib3+j) = sinc(ib1+i,ib3+j) + sk(i,j+12) sinc(ib1+i,ib4+j) = sinc(ib1+i,ib4+j) + sk(i,j+18) sinc(ib2+i,ib1+j) = sinc(ib2+i,ib1+j) + sk(i+6,j) sinc(ib2+i,ib2+j) = sinc(ib2+i,ib2+j) + sk(i+6,j+6) sinc(ib2+i,ib3+j) = sinc(ib2+i,ib3+j) + sk(i+6,j+12) sinc(ib2+i,ib4+j) = sinc(ib2+i,ib4+j) + sk(i+6,j+18) sinc(ib3+i,ib1+j) = sinc(ib3+i,ib1+j) + sk(i+12,j) sinc(ib3+i,ib2+j) = sinc(ib3+i,ib2+j) + sk(i+12,j+6) sinc(ib3+i,ib3+j) = sinc(ib3+i,ib3+j) + sk(i+12,j+12) sinc(ib3+i,ib4+j) = sinc(ib3+i,ib4+j) + sk(i+12,j+18) sinc(ib4+i,ib1+j) = sinc(ib4+i,ib1+j) + sk(i+18,j) sinc(ib4+i,ib2+j) = sinc(ib4+i,ib2+j) + sk(i+18,j+6) sinc(ib4+i,ib3+j) = sinc(ib4+i,ib3+j) + sk(i+18,j+12) sinc(ib4+i,ib4+j) = sinc(ib4+i,ib4+j) + sk(i+18,j+18) 270 continue 2100 continue c do 280 i=1,ndof sinc(ndof+1,i) = ddl(i) sinc(i,ndof+1) = -forc(i) 280 continue sinc(ndof+1,ndof+1) = p ba(ndof+1) = 0.d0 do 290 i=1,ndof do 290 j=1,ndof sinc(i,j) = bc(i) * sinc(i,j) sinc(j,i) = bc(i) * sinc(j,i) ba(i) = bc(i) * ba(i) 290 continue do 300 i=1,ndof if(sinc(i,i) .ne. 0.d0) go to 300 sinc(i,i) = 1.d0 300 continue do 700 i=1,ndof do 700 j=1,ndof snr2(i,j) = sinc(i,j) 700 continue c call solve (sinc,dd,ba,dp,det,idt,ndof,ndof+1,jdt) dp = dd(ndof+1) go to 1100 c c******************************************************************** c 計算結果表示 c******************************************************************** 1300 continue if( abs((p1-p)/p).lt.1.d-7 ) then write(8,*) write(8,5100) write(8,5111) p 5100 format('---------------------- bucled ------------------------') 5111 format('converged p = ',1pe15.7) stop end if c arc = 0.d0 nlet = 0 if( konv.gt.1 ) rn = rn * 1.5d0 if( jdt.eq.0 ) nccl = 1 write(*,*)'----------------------------------------------------' write(*,3010) jdt, p 3010 format(3x,'jdt = ',i2,5x,'p = ',1pe15.7) c c do 950 i=1,np write(20+i,5001) d(mds(i), mmds(i)), p 950 continue c 5000 format(i3,3x,1pd15.7,3x,1pd15.7,3x,1pd15.7) 5001 format(1pe15.7,5x,1pe15.7) c if( jdt.ge.1 ) then write(*,*)'buckled.....' nbcl = 1 nacl = 0 end if c 4020 continue p1 = p if( nacl.eq.1 ) rn = rn / 2.d0 if( nbcl.eq.1 ) then nbcl = 0 nacl = 1 go to 4000 end if go to 4010 c 4000 continue read(3) dd,d,ddl,dp,dpl,r0,p rewind 3 rn = rn / 2.d0 c if( nccl.eq.0 ) then r = r * 0.9d0 r0 = r dpl = r0 end if c 4010 continue if( nccl.eq.1 ) then write(3) dd,d,ddl,dp,dpl,r0,p rewind 3 r = rn end if c c if ( p.gt.8000) then c stop c end if c goto 1000 end c c---------------------- main program end ---------------------------- c c==================================================================== subroutine stmx (s,sk,a,b,d,dd,ba,alp) implicit real*8(a-h,o-z) c dimension d(24), dl(24), ta(24,24), tb(24,24), tg(24,24), & dt0(24,24,24), sk(24,24), cf1(6,6),te(24,24), & cf2(6,6), cf3(6,6), cf4(6,6), cf5(6,6), cf6(6,6), & tngnt1(24,24), tngnt2(24,24), tngnt3(24,24), & tl(24,24), s(24,24), alp(3), & dmy1(24), dmy2(24), dmy3(24,24), dmy4(24,24), & dmy5(24,24), dmy6(24,24), cf(24,24), sk1(24,24), & ba(24), dd(24), sk2(24,24), ti(24,24), sk3(24,24), & dmy7(24), dmy8(24) c c******************************************************************** c call clrv (24,dmy1) call clrmx (24,tl) call clrmx (24,te) call clrmx (24,ti) call trns (ti,alp(1),alp(2),alp(3)) call euler (te,d(4),d(5),d(10),d(11),d(16),d(17),d(22),d(23)) c call mxtv (24,ti,dd,dmy1) call mxv (24,te,dmy1,dd) c do 130 i=1,24 d(i) = d(i) + dd(i) 130 continue c c******************************************************************** c u1 = d(1) v1 = d(2) w1 = d(3) a1 = d(4) b1 = d(5) g1 = d(6) u2 = d(7) v2 = d(8) w2 = d(9) a2 = d(10) b2 = d(11) g2 = d(12) u3 = d(13) v3 = d(14) w3 = d(15) a3 = d(16) b3 = d(17) g3 = d(18) u4 = d(19) v4 = d(20) w4 = d(21) a4 = d(22) b4 = d(23) g4 = d(24) c ca1 = cos(a1) cb1 = cos(b1) cg1 = cos(g1) ca2 = cos(a2) cb2 = cos(b2) cg2 = cos(g2) ca3 = cos(a3) cb3 = cos(b3) cg3 = cos(g3) ca4 = cos(a4) cb4 = cos(b4) cg4 = cos(g4) c sa1 = sin(a1) sb1 = sin(b1) sg1 = sin(g1) sa2 = sin(a2) sb2 = sin(b2) sg2 = sin(g2) sa3 = sin(a3) sb3 = sin(b3) sg3 = sin(g3) sa4 = sin(a4) sb4 = sin(b4) sg4 = sin(g4) c a12 = (a1+a2) / 2.d0 a13 = (a1+a3) / 2.d0 a14 = (a1+a4) / 2.d0 b12 = (b1+b2) / 2.d0 b13 = (b1+b3) / 2.d0 b14 = (b1+b4) / 2.d0 g12 = (g1+g2) / 2.d0 g13 = (g1+g3) / 2.d0 g14 = (g1+g4) / 2.d0 c caa2 = cos(a12) caa3 = cos(a13) caa4 = cos(a14) cbb2 = cos(b12) cbb3 = cos(b13) cbb4 = cos(b14) cgg2 = cos(g12) cgg3 = cos(g13) cgg4 = cos(g14) saa2 = sin(a12) saa3 = sin(a13) saa4 = sin(a14) sbb2 = sin(b12) sbb3 = sin(b13) sbb4 = sin(b14) sgg2 = sin(g12) sgg3 = sin(g13) sgg4 = sin(g14) c c******************************************************************** c call clrv (24,dl) c dl(7) = u2-u1+a-a*cg1*cb1 dl(8) = v2-v1-a*sg1*ca1-a*cg1*sb1*sa1 dl(9) = w2-w1-a*sg1*sa1+a*cg1*sb1*ca1 dl(10) = (a2-a1)+(g2-g1)*sbb2 dl(11) = -(g2-g1)*cbb2*saa2+(b2-b1)*caa2 dl(12) = (g2-g1)*cbb2*caa2+(b2-b1)*saa2 dl(13) = u3-u1+a-a*cg1*cb1+b*sg1*cb1 dl(14) = v3-v1-a*sg1*ca1-a*cg1*sb1*sa1+b-b*cg1*ca1+b*sg1*sb1*sa1 dl(15) = w3-w1-a*sg1*sa1+a*cg1*sb1*ca1-b*cg1*sa1-b*sg1*sb1*ca1 dl(16) = (a3-a1)+(g3-g1)*sbb3 dl(17) = -(g3-g1)*cbb3*saa3+(b3-b1)*caa3 dl(18) = (g3-g1)*cbb3*caa3+(b3-b1)*saa3 dl(19) = u4-u1+b*sg1*cb1 dl(20) = v4-v1+b-b*cg1*ca1+b*sg1*sb1*sa1 dl(21) = w4-w1-b*cg1*sa1-b*sg1*sb1*ca1 dl(22) = (a4-a1)+(g4-g1)*sbb4 dl(23) = -(g4-g1)*cbb4*saa4+(b4-b1)*caa4 dl(24) = (g4-g1)*cbb4*caa4+(b4-b1)*saa4 c c******************************************************************** c call clrmx (24,ta) call clrmx (24,tb) call clrmx (24,tg) c ta(1,1) = 0.d0 ta(1,2) = 0.d0 ta(1,3) = 0.d0 ta(2,1) = -sg1*sa1+cg1*sb1*ca1 ta(2,2) = -cg1*sa1-ca1*sg1*sb1 ta(2,3) = -cb1*ca1 ta(3,1) = sg1*ca1+cg1*sb1*sa1 ta(3,2) = ca1*cg1-sg1*sb1*sa1 ta(3,3) = -cb1*sa1 tb(1,1) = -cg1*sb1 tb(1,2) = sb1*sg1 tb(1,3) = cb1 tb(2,1) = cg1*cb1*sa1 tb(2,2) = -sa1*sg1*cb1 tb(2,3) = sb1*sa1 tb(3,1) = -cg1*cb1*ca1 tb(3,2) = sg1*cb1*ca1 tb(3,3) = -sb1*ca1 tg(1,1) = -sg1*cb1 tg(1,2) = -cg1*cb1 tg(1,3) = 0.d0 tg(2,1) = cg1*ca1-sg1*sb1*sa1 tg(2,2) = -sg1*ca1-cg1*sb1*sa1 tg(2,3) = 0.d0 tg(3,1) = cg1*sa1+sg1*sb1*ca1 tg(3,2) = -sg1*sa1+cg1*sb1*ca1 tg(3,3) = 0.d0 c do 10 i=1,3 do 10 j=1,3 ta(i+3,j+3) = ta(i,j) ta(i+6,j+6) = ta(i,j) ta(i+9,j+9) = ta(i,j) ta(i+12,j+12) = ta(i,j) ta(i+15,j+15) = ta(i,j) ta(i+18,j+18) = ta(i,j) ta(i+21,j+21) = ta(i,j) tb(i+3,j+3) = tb(i,j) tb(i+6,j+6) = tb(i,j) tb(i+9,j+9) = tb(i,j) tb(i+12,j+12) = tb(i,j) tb(i+15,j+15) = tb(i,j) tb(i+18,j+18) = tb(i,j) tb(i+21,j+21) = tb(i,j) tg(i+3,j+3) = tg(i,j) tg(i+6,j+6) = tg(i,j) tg(i+9,j+9) = tg(i,j) tg(i+12,j+12) = tg(i,j) tg(i+15,j+15) = tg(i,j) tg(i+18,j+18) = tg(i,j) tg(i+21,j+21) = tg(i,j) 10 continue c call clmx3 (24,dt0) c do 20 i=1,24 do 20 j=1,24 dt0(4,i,j) = ta(i,j) dt0(5,i,j) = tb(i,j) dt0(6,i,j) = tg(i,j) 20 continue c c******************************************************************** c coefficient of relative displacement vector c call clrmx (6,cf1) call clrmx (6,cf2) call clrmx (6,cf3) call clrmx (6,cf4) call clrmx (6,cf5) call clrmx (6,cf6) c cf1(1,1) = -1.d0 cf1(1,5) = a*cg1*sb1 cf1(1,6) = a*sg1*cb1 cf1(2,2) = -1.d0 cf1(2,4) = a*sg1*sa1-a*cg1*sb1*ca1 cf1(2,5) = -a*cg1*cb1*sa1 cf1(2,6) = -a*cg1*ca1+a*sg1*sb1*sa1 cf1(3,3) = -1.d0 cf1(3,4) = -a*sg1*ca1-a*cg1*sb1*sa1 cf1(3,5) = a*cg1*cb1*ca1 cf1(3,6) = -a*cg1*sa1-a*sg1*sb1*ca1 cf1(4,4) = -1.d0 cf1(4,5) = (g2-g1)/2.d0*cbb2 cf1(4,6) = -sbb2 cf1(5,4) = -(g2-g1)/2.d0*cbb2*caa2-(b2-b1)/2.d0*saa2 cf1(5,5) = (g2-g1)/2.d0*sbb2*saa2-caa2 cf1(5,6) = cbb2*saa2 cf1(6,4) = -(g2-g1)/2.d0*cbb2*saa2+(b2-b1)/2.d0*caa2 cf1(6,5) = -saa2-(g2-g1)/2.d0*sbb2*caa2 cf1(6,6) = -cbb2*caa2 c cf2(1,1) = 1.d0 cf2(2,2) = 1.d0 cf2(3,3) = 1.d0 cf2(4,4) = 1.d0 cf2(4,5) = cf1(4,5) cf2(4,6) = -cf1(4,6) cf2(5,4) = cf1(5,4) cf2(5,5) = (g2-g1)/2.d0*sbb2*saa2+caa2 cf2(5,6) = -cf1(5,6) cf2(6,4) = cf1(6,4) cf2(6,5) = -(g2-g1)/2.d0*sbb2*caa2+saa2 cf2(6,6) = -cf1(6,6) c cf3(1,1) = -1.d0 cf3(1,5) = a*cg1*sb1-b*sg1*sb1 cf3(1,6) = a*sg1*cb1+b*cg1*cb1 cf3(2,2) = -1.d0 cf3(2,4) = a*sg1*sa1-a*cg1*sb1*ca1+b*cg1*sa1+b*sg1*sb1*ca1 cf3(2,5) = -a*cg1*cb1*sa1+b*sg1*cb1*sa1 cf3(2,6) = -a*cg1*ca1+a*sg1*sb1*sa1+b*sg1*ca1+b*cg1*sb1*sa1 cf3(3,3) = -1.d0 cf3(3,4) = -a*sg1*ca1-a*cg1*sb1*sa1-b*cg1*ca1+b*sg1*sb1*sa1 cf3(3,5) = a*cg1*cb1*ca1-b*sg1*cb1*ca1 cf3(3,6) = -a*cg1*sa1-a*sg1*sb1*ca1+b*sg1*sa1-b*cg1*sb1*ca1 cf3(4,4) = -1.d0 cf3(4,5) = (g3-g1)/2.d0*cbb3 cf3(4,6) = -sbb3 cf3(5,4) = -(g3-g1)/2.d0*cbb3*caa3-(b3-b1)/2.d0*saa3 cf3(5,5) = (g3-g1)/2.d0*sbb3*saa3-caa3 cf3(5,6) = cbb3*saa3 cf3(6,4) = -(g3-g1)/2.d0*cbb3*saa3+(b3-b1)/2.d0*caa3 cf3(6,5) = -saa3-(g3-g1)/2.d0*sbb3*caa3 cf3(6,6) = -cbb3*caa3 c cf4(1,1) = 1.d0 cf4(2,2) = 1.d0 cf4(3,3) = 1.d0 cf4(4,4) = 1.d0 cf4(4,5) = cf3(4,5) cf4(4,6) = -cf3(4,6) cf4(5,4) = cf3(5,4) cf4(5,5) = (g3-g1)/2.d0*sbb3*saa3+caa3 cf4(5,6) = -cf3(5,6) cf4(6,4) = cf3(6,4) cf4(6,5) = -(g3-g1)/2.d0*sbb3*caa3+saa3 cf4(6,6) = -cf3(6,6) c cf5(1,1) = -1.d0 cf5(1,5) = -b*sg1*sb1 cf5(1,6) = b*cg1*cb1 cf5(2,2) = -1.d0 cf5(2,4) = b*cg1*sa1+b*sg1*sb1*ca1 cf5(2,5) = b*sg1*cb1*sa1 cf5(2,6) = b*sg1*ca1+b*cg1*sb1*sa1 cf5(3,3) = -1.d0 cf5(3,4) = -b*cg1*ca1+b*sg1*sb1*sa1 cf5(3,5) = -b*sg1*cb1*ca1 cf5(3,6) = b*sg1*sa1-b*cg1*sb1*ca1 cf5(4,4) = -1.d0 cf5(4,5) = (g4-g1)/2.d0*cbb4 cf5(4,6) = -sbb4 cf5(5,4) = -(g4-g1)/2.d0*cbb4*caa4-(b4-b1)/2.d0*saa4 cf5(5,5) = (g4-g1)/2.d0*sbb4*saa4-caa4 cf5(5,6) = cbb4*saa4 cf5(6,4) = -(g4-g1)/2.d0*cbb4*saa4+(b4-b1)/2.d0*caa4 cf5(6,5) = -saa4-(g4-g1)/2.d0*sbb4*caa4 cf5(6,6) = -cbb4*caa4 c cf6(1,1) = 1.d0 cf6(2,2) = 1.d0 cf6(3,3) = 1.d0 cf6(4,4) = 1.d0 cf6(4,5) = cf5(4,5) cf6(4,6) = -cf5(4,6) cf6(5,4) = cf5(5,4) cf6(5,5) = (g4-g1)/2.d0*sbb4*saa4+caa4 cf6(5,6) = -cf5(5,6) cf6(6,4) = cf5(6,4) cf6(6,5) = -(g4-g1)/2.d0*sbb4*caa4+saa4 cf6(6,6) = -cf5(6,6) c call clrmx (24,cf) do 1000 i=1,6 do 1000 j=1,6 cf(i+6,j) = cf1(i,j) cf(i+6,j+6) = cf2(i,j) cf(i+12,j) = cf3(i,j) cf(i+12,j+12) = cf4(i,j) cf(i+18,j) = cf5(i,j) cf(i+18,j+18) = cf6(i,j) 1000 continue c call trns(tl,d(4),d(5),d(6)) c call clrmx (24,tngnt1) call clrmx (24,tngnt2) call clrmx (24,tngnt3) call clrv (24,dmy1) call clrv (24,dmy2) call clrmx (24,dmy3) call clrmx (24,dmy4) call clrmx (24,dmy5) call clrmx (24,dmy6) c c******************************************************************** c call mxtv (24,tl,dl,dmy1) call mxv (24,s,dmy1,dmy2) do 30 k=1,24 do 30 i=1,24 do 30 j=1,24 tngnt1(i,k) = tngnt1(i,k) + dt0(k,i,j) * dmy2(j) 30 continue c c******************************************************************** c do 40 k=1,24 do 40 i=1,24 do 40 j=1,24 dmy3(i,k) = dmy3(i,k) + dt0(k,j,i) * dl(j) 40 continue c call mxmx (24,s,dmy3,dmy4) call mxmx (24,tl,dmy4,tngnt2) c c******************************************************************** c call mxtmx (24,tl,cf,dmy5) call mxmx (24,s,dmy5,dmy6) call mxmx(24,tl,dmy6,tngnt3) c c******************************************************************** c call clrmx (24,sk1) call clrmx (24,sk2) call clrmx (24,sk3) call clrmx (24,sk) c do 80 i=1,24 do 80 j=1,24 sk1(i,j) = tngnt1(i,j) + tngnt2(i,j) + tngnt3(i,j) 80 continue c call euler (te,d(4),d(5),d(10),d(11),d(16),d(17),d(22),d(23)) call mxmx (24,ti,sk1,sk2) call mxmx (24,sk2,te,sk3) call mxmxt (24,sk3,ti,sk) c c******************************************************************** c call clrv (24,dmy1) call clrv (24,dmy2) call clrv (24,dmy7) call clrv (24,dmy8) call mxtv (24,tl,dl,dmy1) call mxv (24,s,dmy1,dmy2) call mxv (24,tl,dmy2,dmy7) call mxv (24,ti,dmy7,dmy8) do 120 i=1,24 ba(i) = ba(i) - dmy8(i) 120 continue c return end c c==================================================================== subroutine mxtmx(n,a,b,c) implicit real*8(a-h,o-z) dimension a(n,n),b(n,n),c(n,n) do 10 i=1,n do 10 k=1,n c(i,k)=0.d0 do 10 j=1,n c(i,k)=c(i,k)+a(j,i)*b(j,k) 10 continue return end c c==================================================================== subroutine mxmx(n,a,b,c) implicit real*8(a-h,o-z) dimension a(n,n),b(n,n),c(n,n) do 10 i=1,n do 10 k=1,n c(i,k)=0.d0 do 10 j=1,n c(i,k)=c(i,k)+a(i,j)*b(j,k) 10 continue return end c c==================================================================== subroutine mxtv(n,a,b,c) implicit real*8(a-h,o-z) dimension a(n,n),b(n),c(n) do 10 i=1,n c(i)=0.d0 do 10 j=1,n c(i)=c(i)+a(j,i)*b(j) 10 continue return end c c==================================================================== subroutine mxv(n,a,b,c) implicit real*8(a-h,o-z) dimension a(n,n),b(n),c(n) do 10 i=1,n c(i)=0.d0 do 10 j=1,n c(i)=c(i)+a(i,j)*b(j) 10 continue return end c c==================================================================== subroutine mxmxt(n,a,b,c) implicit real*8(a-h,o-z) dimension a(n,n),b(n,n),c(n,n) do 10 i=1,n do 10 k=1,n c(i,k)=0.d0 do 10 j=1,n c(i,k)=c(i,k)+a(i,j)*b(k,j) 10 continue return end c c==================================================================== subroutine smx(s,th,e,um,a,b) implicit real*8(a-h,o-z) dimension s(24,24) c call clrmx (24,s) c p1 = ( b / ( 3.d0 * a ) + ( 1.d0 - um ) * a / ( 6.d0 * b ) ) & * e * th / ( 1.d0 - um * um ) p2 = ( um / 4.d0 + ( 1 - um ) / 8.d0 ) & * e * th / ( 1.d0 - um * um ) p3 = ( a / ( 3.d0 * b ) + ( 1.d0 - um ) * b / ( 6.d0 * a ) ) & * e * th / ( 1.d0 - um * um ) p4 = ( - b / ( 3.d0 * a ) + ( 1.d0 - um ) * a / ( 12.d0 * b ) ) & * e * th / ( 1 - um * um ) p5 = ( um / 4.d0 - ( 1.d0 - um ) / 8.d0 ) & * e * th / ( 1.d0 - um * um ) p6 = ( a / ( 6.d0 * b ) - ( 1.d0 - um ) * b / ( 6.d0 & * a ) ) * e * th / ( 1.d0 - um * um ) p7 = ( b / ( 6.d0 * a ) - ( 1.d0 - um ) * a / ( 6.d0 * b ) ) & * e * th / ( 1.d0 - um * um ) p8 = ( - a / ( 3.d0 * b ) + ( 1.d0 - um ) * b / ( 12.d0 * a ) ) & * e * th / ( 1.d0 - um * um ) c s(1,1) = p1 s(1,7) = p4 s(1,13) = -p1 / 2.d0 s(1,19) = p7 s(1,2) = p2 s(1,8) = p5 s(1,14) = -p2 s(1,20) = -p5 s(7,1) = p4 s(7,7) = p1 s(7,13) = p7 s(7,19) = -p1 / 2.d0 s(7,2) = -p5 s(7,8) = -p2 s(7,14) = p5 s(7,20) = p2 s(13,1) = -p1 / 2.d0 s(13,7) = p7 s(13,13) = p1 s(13,19) = p4 s(13,2) = -p2 s(13,8) = -p5 s(13,14) = p2 s(13,20) = p5 s(19,1) = p7 s(19,7) = -p1 / 2.d0 s(19,13) = p4 s(19,19) = p1 s(19,2) = p5 s(19,8) = p2 s(19,14) = -p5 s(19,20) = -p2 s(2,1) = p2 s(2,7) = -p5 s(2,13) = -p2 s(2,19) = p5 s(2,2) = p3 s(2,8) = p6 s(2,14) = -p3 / 2.d0 s(2,20) = p8 s(8,1) = p5 s(8,7) = -p2 s(8,13) = -p5 s(8,19) = p2 s(8,2) = p6 s(8,8) = p3 s(8,14) = p8 s(8,20) = -p3 / 2.d0 s(14,1) = -p2 s(14,7) = p5 s(14,13) = p2 s(14,19) = -p5 s(14,2) = -p3 / 2.d0 s(14,8) = p8 s(14,14) = p3 s(14,20) = p6 s(20,1) = -p5 s(20,7) = p2 s(20,13) = p5 s(20,19) = -p2 s(20,2) = p8 s(20,8) = -p3 / 2.d0 s(20,14) = p6 s(20,20) = p3 c p15 = e * th * th * th / ( 360.d0 * ( 1.d0 - um ** 2 ) * a * b ) s(3,3) = 120.d0 * b * b / a / a + 120.d0 * a * a / b / b & + 84.d0 - 24.d0 * um s(9,3) = - 120.d0 * b * b / a / a + 60.d0 * a * a / b / b & - 84.d0 + 24.d0 * um s(15,3) = - 60.d0 * b * b / a / a - 60.d0 * a * a / b / b & + 84.d0 - 24.d0 * um s(21,3) = 60.d0 * b * b / a / a - 120.d0 * a * a / b / b & - 84.d0 + 24.d0 * um s(3,3) = s(3,3) * p15 s(9,3) = s(9,3) * p15 s(15,3) = s(15,3) * p15 s(21,3) = s(21,3) * p15 s(9,9) = s(3,3) s(15,15) = s(3,3) s(21,21) = s(3,3) s(3,9) = s(9,3) s(21,15) = s(9,3) s(15,21) = s(9,3) s(21,9) = s(15,3) s(3,15) = s(15,3) s(9,21) = s(15,3) s(15,9) = s(21,3) s(9,15) = s(21,3) s(3,21) = s(21,3) c s(4,3) = 60.d0 * a * a / b + 6.d0 * b + 24.d0 * b * um s(10,3) = 30.d0 * a * a / b - 6.d0 * b - 24.d0 * b * um s(16,3) = 30.d0 * a * a / b - 6.d0 * b + 6.d0 * b * um s(22,3) = 60.d0 * a * a / b + 6.d0 * b - 6.d0 * b * um s(4,3) = s(4,3) * p15 s(10,3) = s(10,3) * p15 s(16,3) = s(16,3) * p15 s(22,3) = s(22,3) * p15 s(10,9) = s(4,3) s(16,15) = - s(4,3) s(22,21) = - s(4,3) s(4,9) = s(10,3) s(22,15) = - s(10,3) s(16,21) = - s(10,3) s(22,9) = s(16,3) s(4,15) = - s(16,3) s(10,21) = - s(16,3) s(16,9) = s(22,3) s(10,15) = - s(22,3) s(4,21) = - s(22,3) s(3,4) = s(4,3) s(3,10) = s(10,3) s(3,16) = s(16,3) s(3,22) = s(22,3) s(9,10) = s(10,9) s(15,16) = s(16,15) s(21,22) = s(22,21) s(9,4) = s(4,9) s(15,22) = s(22,15) s(21,16) = s(16,21) s(9,22) = s(22,9) s(15,4) = s(4,15) s(21,10) = s(10,21) s(9,16) = s(16,9) s(15,10) = s(10,15) s(21,4) = s(4,21) c s(5,3) = - 60.d0 * b * b / a - 6.d0 * a - 24.d0 * um * a s(11,3) = - 60.d0 * b * b / a - 6.d0 * a + 6.d0 * um * a s(17,3) = - 30.d0 * b * b / a + 6.d0 * a - 6.d0 * um * a s(23,3) = - 30.d0 * b * b / a + 6.d0 * a + 24.d0 * um * a s(5,3) = s(5,3) * p15 s(11,3) = s(11,3) * p15 s(17,3) = s(17,3) * p15 s(23,3) = s(23,3) * p15 s(11,9) = - s(5,3) s(17,15) = - s(5,3) s(23,21) = s(5,3) s(5,9) = - s(11,3) s(23,15) = - s(11,3) s(17,21) = s(11,3) s(23,9) = - s(17,3) s(5,15) = - s(17,3) s(11,21) = s(17,3) s(17,9) = - s(23,3) s(11,15) = - s(23,3) s(5,21) = s(23,3) s(3,5) = s(5,3) s(3,11) = s(11,3) s(3,17) = s(17,3) s(3,23) = s(23,3) s(9,11) = s(11,9) s(15,17) = s(17,15) s(21,23) = s(23,21) s(9,5) = s(5,9) s(15,23) = s(23,15) s(21,17) = s(17,21) s(9,23) = s(23,9) s(15,5) = s(5,15) s(21,11) = s(11,21) s(9,17) = s(17,9) s(15,11) = s(11,15) s(21,5) = s(5,21) c s(4,4) = 8.d0 * b * b + 40.d0 * a * a - 8.d0 * b * b * um s(10,4) = - 8.d0 * b * b + 20.d0 * a * a + 8.d0 * b * b * um s(16,4) = 2.d0 * b * b + 10.d0 * a * a - 2.d0 * b * b * um s(22,4) = - 2.d0 * b * b + 20.d0 * a * a + 2.d0 * b * b * um s(4,4) = s(4,4) * p15 s(10,4) = s(10,4) * p15 s(16,4) = s(16,4) * p15 s(22,4) = s(22,4) * p15 s(10,10) = s(4,4) s(16,16) = s(4,4) s(22,22) = s(4,4) s(4,10) = s(10,4) s(22,16) = s(10,4) s(16,22) = s(10,4) s(22,10) = s(16,4) s(4,16) = s(16,4) s(10,22) = s(16,4) s(16,10) = s(22,4) s(10,16) = s(22,4) s(4,22) = s(22,4) c s(5,4) = - 30.d0 * um * a * b s(11,4) = 0.d0 s(17,4) = 0.d0 s(23,4) = 0.d0 s(5,4) = s(5,4) * p15 s(11,4) = s(11,4) * p15 s(17,4) = s(17,4) * p15 s(23,4) = s(23,4) * p15 s(11,10) = - s(5,4) s(17,16) = s(5,4) s(23,22) = - s(5,4) s(5,10) = - s(11,4) s(23,16) = s(11,4) s(17,22) = - s(11,4) s(23,10) = - s(17,4) s(5,16) = s(17,4) s(11,22) = - s(17,4) s(17,10) = - s(23,4) s(11,16) = s(23,4) s(5,22) = - s(23,4) s(4,5) = s(5,4) s(4,11) = s(11,4) s(4,17) = s(17,4) s(4,23) = s(23,4) s(10,11) = s(11,10) s(16,17) = s(17,16) s(22,23) = s(23,22) s(10,5) = s(5,10) s(16,23) = s(23,16) s(22,17) = s(17,22) s(10,23) = s(23,10) s(16,5) = s(5,16) s(22,11) = s(11,22) s(10,17) = s(17,10) s(16,11) = s(11,16) s(22,5) = s(5,22) c s(5,5) = 40.d0 * b * b + 8.d0 * a * a - 8.d0 * a * a * um s(11,5) = 20.d0 * b * b - 2.d0 * a * a + 2.d0 * a * a * um s(17,5) = 10.d0 * b * b + 2.d0 * a * a - 2.d0 * a * a * um s(23,5) = 20.d0 * b * b - 8.d0 * a * a + 8.d0 * a * a * um s(5,5) = s(5,5) * p15 s(11,5) = s(11,5) * p15 s(17,5) = s(17,5) * p15 s(23,5) = s(23,5) * p15 s(11,11) = s(5,5) s(17,17) = s(5,5) s(23,23) = s(5,5) s(5,11) = s(11,5) s(23,17) = s(11,5) s(17,23) = s(11,5) s(23,11) = s(17,5) s(5,17) = s(17,5) s(11,23) = s(17,5) s(17,11) = s(23,5) s(11,17) = s(17,11) s(5,23) = s(17,11) c p9 = e * th * a * b * 0.03d0 s(6, 6) = p9 s(6, 12) = -p9 / 3.d0 s(6, 18) = -p9 / 3.d0 s(6, 24) = -p9 / 3.d0 s(12, 6) = -p9 / 3.d0 s(12, 12) = p9 s(12, 18) = -p9 / 3.d0 s(12, 24) = -p9 / 3.d0 s(18, 6) = -p9 / 3.d0 s(18, 12) = -p9 / 3.d0 s(18, 18) = p9 s(18, 24) = -p9 / 3.d0 s(24, 6) = -p9 / 3.d0 s(24, 12) = -p9 / 3.d0 s(24, 18) = -p9 / 3.d0 s(24, 24) = p9 c return end c==================================================================== subroutine trns (tl,al,bt,gm) implicit real*8(a-h,o-z) dimension t(3,3), tl(24,24) c do 110 i=1,24 do 110 j=1,24 tl(i,j) = 0.d0 110 continue c check!! t(1,1) = cos(gm)*cos(bt) t(1,2) = -sin(gm)*cos(bt) t(1,3) = sin(bt) t(2,1) = sin(gm)*cos(al)+cos(gm)*sin(bt)*sin(al) t(2,2) = cos(gm)*cos(al)-sin(gm)*sin(bt)*sin(al) t(2,3) = -cos(bt)*sin(al) t(3,1) = sin(gm)*sin(al)-cos(gm)*sin(bt)*cos(al) t(3,2) = cos(gm)*sin(al)+sin(gm)*sin(bt)*cos(al) t(3,3) = cos(bt)*cos(al) c do 10 i=1,3 do 10 j=1,3 tl(i,j) = t(i,j) tl(i+3,j+3) = t(i,j) tl(i+6,j+6) = t(i,j) tl(i+9,j+9) = t(i,j) tl(i+12,j+12) = t(i,j) tl(i+15,j+15) = t(i,j) tl(i+18,j+18) = t(i,j) tl(i+21,j+21) = t(i,j) 10 continue c return end c c==================================================================== subroutine euler(te,a1,b1,a2,b2,a3,b3,a4,b4) implicit real*8 (a-h,o-z) dimension te(24,24) c do 10 i=1,24 do 10 j=1,24 te(i,j) = 0.d0 10 continue c do 20 i=1,24 te(i,i)=1.d0 20 continue c te(4,5) = sin(a1)*tan(b1) te(4,6) = -cos(a1)*tan(b1) te(5,5) = cos(a1) te(5,6) = sin(a1) te(6,5) = -sin(a1)/cos(b1) te(6,6) = cos(a1)/cos(b1) te(10,11) = sin(a2)*tan(b2) te(10,12) = -cos(a2)*tan(b2) te(11,11) = cos(a2) te(11,12) = sin(a2) te(12,11) = -sin(a2)/cos(b2) te(12,12) = cos(a2)/cos(b2) te(16,17) = sin(a3)*tan(b3) te(16,18) = -cos(a3)*tan(b3) te(17,17) = cos(a3) te(17,18) = -sin(a3) te(18,17) = -sin(a3)/cos(b3) te(18,18) = cos(a3)/cos(b3) te(22,23) = sin(a4)*tan(b4) te(22,24) = -cos(a4)*tan(b4) te(23,23) = cos(a4) te(23,24) = sin(a4) te(24,23) = -sin(a4)/cos(b4) te(24,24) = cos(a4)/cos(b4) return end c c==================================================================== subroutine clrv(n,a) implicit real*8(a-h,o-z) dimension a(n) do 10 i=1,n a(i)=0.d0 10 continue return end c c==================================================================== subroutine clmx3 (n,a) implicit real*8(a-h,o-z) dimension a(n,n,n) do 10 k=1,n do 10 i=1,n do 10 j=1,n a(k,i,j) = 0.d0 10 continue return end c c==================================================================== subroutine clrmx (n,a) implicit real*8(a-h,o-z) c dimension a(n,n) do 10 i=1,n do 10 j=1,n a(i,j)=0.d0 10 continue c return end c==================================================================== subroutine solve (sinc,dd,df,dp,det,idt,n0,n01,jdt) c c Gaussian Elimination Method c implicit real*8 (a-h,o-z) dimension sinc(n01,n01),dd(n01),df(n01) c det=1.d0 idt=0 jdt=0 c Forward Elimination do 100 i=1,n0 inext=i+1 if(sinc(i,i).le.0.d0) jdt=jdt+1 det=det*sinc(i,i) j=idint(dlog10(dabs(det))) idt=idt+j det=det/10.d0**(j) c sinc(i,i)=1.d0/sinc(i,i) df(i)=df(i)*sinc(i,i) do 110 j=inext,n01 sinc(i,j)=sinc(i,j)*sinc(i,i) df(j)=df(j)-sinc(j,i)*df(i) do 110 k=inext,n01 sinc(k,j)=sinc(k,j)-sinc(k,i)*sinc(i,j) 110 continue 100 continue c Backward Substitution dd(n01)=df(n01)/sinc(n01,n01) dp=dd(n01) do 120 i=1,n0 k=n01-i dd(k)=df(k) knext=k+1 do 120 j=knext,n01 dd(k)=dd(k)-sinc(k,j)*dd(j) 120 continue c return end