wilson时程积分法计算结构的地震响应时程的程序
$debugIMPLICIT REAL*8(A-H,O-Z)
DIMENSION AMZL(7,7),AKK1(7,7),AKGD(7,7),AMZL1(7,7)
DIMENSION XX(7),XXG0(2720),XXG(2720)
DIMENSION FF(7,1),XX2PP(2720,7),XX1PP(2720,7),XXPP(2720,7)
DIMENSION xpie1(7,1),xpie2(7,1),xpie(7,1)
DIMENSION guodu1(7,1),guodu2(7,1),guodu(7,1)
DIMENSION IS(7),JS(7)
DOUBLE PRECISION AKK1
DOUBLE PRECISION AMZL,AKGD
DATA AMZL/4818000.,0,0,0,0,0,0,
* 0,1475000.,0,0,0,0,0,
* 0,0,1302000.,0,0,0,0,
* 0,0,0,1533000.,0,0,0,
* 0,0,0,0,1840000.,0,0,
* 0,0,0,0,0,2205000.,0,
* 0,0,0,0,0,0,3738000./
DATA AKGD/655000000.,-366000000.,0,0,0,0,0,
* -366000000.,734000000.,-368000000.,0,0,0,0,
* 0,-368000000.,778000000.,-410000000.,0,0,0,
* 0,0,-410000000.,901000000.,-491000000.,0,0,
* 0,0,0,-491000000.,1086000000.,-595000000.,0,
* 0,0,0,0,-595000000.,1344000000.,-749000000.,
* 0,0,0,0,0,-749000000.,749000000./
C DATA II/1,1,1,1,1,1,1/
THITA=1.4
DATAT=0.02
alfa=0.1956
bata=0.0095
DO 75 I=1,7
DO 75 J=1,7
AMZL1(I,J)=AMZL(I,J)*(6./(THITA**2*DATAT**2)+
* 3.*alfa/(THITA*DATAT))
75 CONTINUE
DO 85 I=1,7
DO 85 J=1,7
AKK1(I,J)=AKGD(I,J)*(1.+3.*bata/(THITA*DATAT))+AMZL1(I,J)
85 CONTINUE
CALL INV(AKK1,7,L,IS,JS)
NMB=100
OPEN(1,FILE='TAFTNE21.txt',STATUS='OLD')
do 18 j=1,NMB
18 read(1,*)XXG0(j)
close(1)
XXG(1)=XXG0(1)
do 28 j=2,NMB
28 XXG(J)=XXG0(j)-XXG0(j-1)
do 38 j=1,7
38 FF(J,1)=-AMZL(J,J)*XXG(1)
CALL MUL(AKK1,FF,7,7,1,XX)
do 48 j=1,7
XX2PP(1,J)=XX(J)*6./(THITA**3.*DATAT**2.)
XX1PP(1,J)=XX2PP(1,J)*DATAT/2.
XXPP(1,J)=XX2PP(1,J)*DATAT**2/6.
48 CONTINUE
do 200 I=2,2720
do 49 j=1,7
xpie2(j,1)=XX2PP(I-1,J)
xpie1(j,1)=XX1PP(I-1,J)
xpie(j,1)=XXPP(I-1,J)
49 CONTINUE
CALL MUL(AKGD,xpie1,7,7,1,guodu1)
CALL MUL(AKGD,xpie2,7,7,1,guodu2)
do 101 j=1,7
101 FF(J,1)=-AMZL(J,J)*XXG(I)+AMZL(J,J)*(XX1PP(I-1,J)*(6./THITA/DATAT+
* 3.*alfa)+XX2PP(I-1,J)*(3.+alfa*THITA*DATAT/2.))+
* bata*(3.*guodu1(j,1)+THITA*DATAT/2.*guodu2(j,1))
CALL MUL(AKK1,FF,7,7,1,XX)
do 108 j=1,7
XX2PP(I,J)=XX(J)*6./(THITA**3*DATAT**2)-6./(THITA**2*DATAT)
* *XX1PP(I-1,J)+(1.-3./THITA)*XX2PP(I-1,J)
XX1PP(I,J)=XX1PP(I-1,J)+(XX2PP(I,J)+XX2PP(I-1,J))*DATAT/2.
XXPP(I,J)=XXPP(I-1,J)+DATAT*XX1PP(I-1,J)+(2.*XX2PP(I-1,J)+
* XX2PP(I,J))*DATAT**2/6.
108 CONTINUE
200 CONTINUE
OPEN(2,FILE='XX2PP.out')
OPEN(3,FILE='XX1PP.out')
OPEN(4,FILE='XXPP.out')
write(2,111)((XX2PP(I,J),J=1,7),I=1,2720)
write(3,111)((XX1PP(I,J),J=1,7),I=1,2720)
write(4,111)((XXPP(I,J),J=1,7),I=1,2720)
111 format(7e13.5)
close(2)
close(3)
close(4)
END
SUBROUTINE INV(A,N,L,IS,JS)
DIMENSION A(N,N),IS(N),JS(N)
DOUBLE PRECISION A,T,D
L=1
DO 100 K=1,N
D=0.0
DO 10 I=K,N
DO 10 J=K,N
IF (ABS(A(I,J)).GT.D) THEN
D=ABS(A(I,J))
IS(K)=I
JS(K)=J
END IF
10 CONTINUE
IF (D+1.0.EQ.1.0) THEN
L=0
WRITE(*,20)
RETURN
END IF
20 FORMAT(1X,'ERR**NOT INV')
DO 30 J=1,N
T=A(K,J)
A(K,J)=A(IS(K),J)
A(IS(K),J)=T
30 CONTINUE
DO 40 I=1,N
T=A(I,K)
A(I,K)=A(I,JS(K))
A(I,JS(K))=T
40 CONTINUE
A(K,K)=1/A(K,K)
DO 50 J=1,N
IF (J.NE.K) THEN
A(K,J)=A(K,J)*A(K,K)
END IF
50 CONTINUE
DO 70 I=1,N
IF (I.NE.K) THEN
DO 60 J=1,N
IF (J.NE.K) THEN
A(I,J)=A(I,J)-A(I,K)*A(K,J)
END IF
60 CONTINUE
END IF
70 CONTINUE
DO 80 I=1,N
IF (I.NE.K) THEN
A(I,K)=-A(I,K)*A(K,K)
END IF
80 CONTINUE
100 CONTINUE
DO 130 K=N,1,-1
DO 110 J=1,N
T=A(K,J)
A(K,J)=A(JS(K),J)
A(JS(K),J)=T
110 CONTINUE
DO 120 I=1,N
T=A(I,K)
A(I,K)=A(I,IS(K))
A(I,IS(K))=T
120 CONTINUE
130 CONTINUE
RETURN
END
SUBROUTINE MUL(A,B,M,N,K,C)
DIMENSION A(M,N),B(N,K),C(M,K)
DOUBLE PRECISION A,B,C
DO 50 I=1,M
DO 50 J=1,K
C(I,J)=0.0
DO 10 L=1,N
C(I,J)=C(I,J)+A(I,L)*B(L,J)
10 CONTINUE
50 CONTINUE
RETURN
END
[ 本帖最后由 风花雪月 于 2006-10-11 15:16 编辑 ]
页:
[1]