马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
- $debug
- IMPLICIT 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 编辑 ] |
发表于 2006-10-10 10:28
提升卡
变色卡
显身卡