- 积分
- 7263
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-12-20
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
本帖最后由 nickbsb 于 2012-4-19 12:54 编辑
program main
!character*10 a,sa1
common /t/X(200)
!write(*,10)
!10 !format(2x,'n=3630,ih=1210')
integer::n=3630,ih=1210
!read(*,*)n,ih
!write(*,*) 'input data file name:'
!read(*,'(a)') a
!write(*,*) 'input output file name:'
!read(*,'(a)') sa1
open (1,file='e:\aam\zonalaam.grd',form='binary')
read (1) (x(i),i=1,n)
close(1)
open (6,file='e:\aamps\SA1.grd',form='binary')
call tsa(n,ih)
close(6)
!stop
end
SUBROUTINE TSA(IDAY,IH)
C******************************************************************************
C TSA (Time Serial Analysis) subroutine do Spectrum Analysis *
C *
C Variable: X(IDAY) Time Serial Data *
C L(18),AKI(18,4) coefficient *
C******************************************************************************
COMMON /T/X(200)
DIMENSION L(18),AKI(18,4),B(4),
& SP(1000),SQ(1000),R1(1000),R(1000)
PARAMETER (PI=3.141592653589793)
DATA L/ 2, 3, 4, 5, 6, 8, 10, 15, 20,
& 30, 40, 50, 60, 80, 100,200,400,1000/
DATA AKI
& /.010, .038, .074, .111, .145, .206, .256, .349,
& .413, .498, .554, .594, .625, .669, .701, .782,
& .843, .899, .052, .117, .178, .229, .272, .342,
& .394, .484, .543, .616, .663, .695, .720, .775,
& .779, .841, .887, .928, 3.000, 2.605, 2.372, 2.214,
& 2.099, 1.938, 1.831, 1.666, 1.570, 1.459, 1.394, 1.350,
& 1.318, 1.274, 1.243, 1.170, 1.119, 1.075, 4.605, 3.782,
& 3.319, 3.017, 2.802, 2.511, 2.321, 2.038, 1.878, 1.696,
& 1.592, 1.523, 1.473, 1.404, 1.358, 1.247, 1.172, 1.107/
N1=IDAY
N=N1-IH+1
M=INT(N/3-1)
C=0.
DO 10 I=1, IH
10 C=C+X(I)
DO 20 I=1, N1-IH
D=C-X(I)+X(I+IH)
X(I)=C/IH
20 C=D
X(N)=C/IH
C=X(1)
D=X(1)
AMX=0.
DO 30 I=1, N
30 AMX=AMX+X(I)
AMX=AMX/N
CGM=0.
DO 40 I=1, N
40 CGM=CGM+(X(I)-AMX)*(X(I)-AMX)
CGM=SQRT(CGM/N)
DO 50 I=1, N
50 X(I)=(X(I)-AMX)/CGM
DO 70 K=0, M
C=0.
DO 60 J=1, N-K
60 C=C+X(J)*X(J+K)
C=C/FLOAT(N-K)
R1(K)=C
70 R(K)=C*(1.0-1.0*K/M)
WRITE (6,'(/24H Coefficient correlation)')
WRITE (6,'(10F7.3)') (R1(K),K=0,M) !相关系数
DO 90 K=0, M
C=COS(PI*K/M)
U1=0.
U2=0.
DO 80 I=M-1, 1, -1
U0=R(I)+2*C*U1-U2
U2=U1
80 U1=U0
D=U1*C-U2
X(M+K+2)=(R(0)+2*D+R(M)*(-1)**K)/M
90 SP(K)=X(M+K+2)
D=((N-1)*R(1)+1)/SQRT(FLOAT(N-2))
C=0.
DO 100 K=0, M
100 C=C+X(M+K+2)
C=C/(M+1)
WRITE (6,'(/15H White spectrum,/11X,F12.6)') C
U0=(1-R(1)*R(1))*C
U1=1+R(1)*R(1)
U2=2*R(1)
DO 110 K=0, M
110 X(2*M+K+3)=U0/(U1-U2*COS(PI*K/M))
WRITE (6,'(/15H Red spectrum)')
WRITE (6,'(10E7.2)') (X(2*M+K+3),K=0,M) !红谱
K=(2*N-M/2)/M
IF (K.LE.1000) THEN
I=18
DO 150 J=1, 4
C1=FLOAT(L(I))/L(I-1)
C2=FLOAT(K)/L(I-1)
150 B(J)=ALOG(C2)/ALOG(C1)*(AKI(I,J)-AKI(I-1,J))+AKI(I-1,J)
END IF
WRITE (6,'(/13H White - 0.95,/16X,F9.6/)') C*B(3)
DO 170 K=2*M+3, 3*M+3
170 SQ(K-2*M-3)=X(K)*B(3)
WRITE (6,'(A6,9X,A4,12X,A4,11X,A5,8X,A6,9X,A4,12X,A4,11X,A5)')
& 'Period','Cal*','0.95','Error','Period','Cal*','0.95','Error'
DO 190 I=1,(M-1)/2
R(I)=SP(I)-SQ(I)
R1(I)=2.*M/I
R(I+(M+1)/2)=SP(I+(M+1)/2)-SQ(I+(M+1)/2)
R1(I+(M+1)/2)=2.*M/(I+(M+1)/2)
190 WRITE (6,'(F6.2,3E16.6,2X,1H|,2X,F6.2,3E16.6)')
& R1(I), SP(I),SQ(I),R(I),
& R1(I+(M+1)/2), SP(I+(M+1)/2),SQ(I+(M+1)/2),R(I+(M+1)/2)
IF (M/2.EQ.(M+1)/2) THEN
R(M/2)=SP(M/2)-SQ(M/2)
R1(M/2)=2.*M/(M/2)
R(M)=SP(M)-SQ(M)
R1(M)=2.*M/M
WRITE (6,'(F6.3,3E16.6,2X,1H|,2X,F6.3,3E16.6)')
& R1(M/2), SP(M/2),SQ(M/2),R(M/2),R1(M), SP(M),SQ(M),R(M)
ELSE
R((M+1)/2)=SP((M+1)/2)-SQ((M+1)/2)
R1((M+1)/2)=2.*(M+1)/2/((M+1)/2)
WRITE (6,'(F6.3,3E16.6)')
& R1((M+1)/2), SP((M+1)/2),SQ((M+1)/2),R((M+1)/2)
END IF
WRITE (6,'(/17H Period Analysis)') !周期分析
DO 200 I=1,M-1
DO 200 J=I+1,M
IF (R(I).LT.R(J)) THEN
C=R(I)
R(I)=R(J)
R(J)=C
C=R1(I)
R1(I)=R1(J)
R1(J)=C
END IF
200 CONTINUE
DO 210 I=1,M
IF (R(I).LT.0.) RETURN
210 WRITE (6,'(16X,F8.3,E20.6)') R1(I),R(I)
RETURN
END
结果是这样
program exception array bounds exceeds
这个程序我改的对吗
|
-
-
30-60.grd
14.18 KB, 下载次数: 2, 下载积分: 金钱 -5
输入的资料文件
|
窥视卡
雷达卡
发表于 2012-4-19 12:43:53
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
发表于 2012-4-19 12:56:37
楼主