该楼层疑似违规已被系统折叠 隐藏此楼查看此楼
PROGRAM MAIN
IMPLICIT NONE
INTEGER I, NUM, NMAX, M, MB, ME, MS, PN
PARAMETER (NMAX=2000)
parameter (PN=1500)
REAL*8 T(0:NMAX-1),X(0:NMAX-1), PERIOD, VV, PB, PS, RP, IP
integer tmax,maxn,tmin,minn
REAL*8 fmax,fmin,tt,df,w(0:PN-1),ss,dd,td(0:NMAX-1)
OPEN (UNIT=2,FILE='DATA.DAT',STATUS='OLD',ACCESS='SEQUENTIAL',
$ FORM='FORMATTED')
NUM=0
20 READ (2,*,END=30) T(NUM), X(NUM)
NUM=NUM+1
GOTO 20
30 CLOSE (2)
****************************normalize********************************
ss=0.
do i=0,num-1
ss=ss+X(i)
enddo
ss=ss/dble(num)
dd=0.
do i=0,num-1
dd=dd+(x(i)-ss)**2
enddo
dd=sqrt(dd/dble(num-1))
do i=0,num-1
x(i)=(x(i)-ss)
enddo
do i=0,num-2
td(i)=t(i+1)-t(i)
enddo
TMAX=MAXN(NUM,T)
TMIN=MINN(NUM,T)
tt=(t(tmax)-t(tmin))
fmax=dble(num-1)/tt*2.
fmin=1./tt/2.
tmin=minn(num-1,td)
fmax=1./td(tmin)/2.
ps =(fmax-fmin)/dble(pn-1)
* pb=1./fmin
* df=(1./t(tmin)-1./t(tmax))/dble(pn-1)
write(*,*) tt,fmax,fmin,ps,pb
do i=0,pn-1
w(i)=(fmin+ps*dble(i))
enddo
OPEN (UNIT=3,FILE='OUT.DAT',STATUS='UNKNOWN',ACCESS='SEQUENTIAL',
$ FORM='FORMATTED')
DO I=0,PN-2
PERIOD=W(I)
* WRITE(*,*) " Frequency:",PERIOD
CALL DCDFT(NUM,T,X,PERIOD,RP,IP)
PERIOD = PERIOD*365.
WRITE(3,*) PERIOD, RP, IP
ENDDO
CLOSE(3)
1001 FORMAT (F16.6,F16.10)
END
SUBROUTINE DCDFT(N, T, F, FREQ, RP, IP)
INTEGER N, I
REAL*8 SUMFLUX, A(0:2), C(0:2)
REAL*8 T(0:N-1), F(0:N-1), H2(0:N-1), H1(0:N-1), H0(0:N-1)
REAL*8 LH0(0:N-1), LH1(0:N-1), LH2(0:N-1)
REAL*8 INNER, FREQ, FPOWER, temp, IP, RP
DO I=0,N-1
H0(I) = 1.
H1(I) = DCOS(2*3.14159265*FREQ*T(I))
H2(I) = DSIN(2*3.14159265*FREQ*T(I))
ENDDO
SUMFLUX = 0.
DO I=0,N-1
SUMFLUX = SUMFLUX+F(I)
ENDDO
DO I=0,N-1
F(I) = F(I)-SUMFLUX/DBLE(N)
ENDDO
A(0) = DSQRT(1./DBLE(N))
A(1) = INNER(N, H1, H1)-(A(0)*INNER(N, H0, H1))**2
A(1) = DSQRT(1./A(1))
A(2) = INNER(N, H2, H2)-(A(0)*INNER(N, H0, H2))**2
$ -(A(1)*INNER(N, H1, H2))**2
A(2) = A(2)-A(1)**2*A(0)**4*(INNER(N, H0, H1))**2
$ *(INNER(N, H0, H2))**2
A(2) = A(2)+2.*(A(0)*A(1))**2*INNER(N, H0, H1)
$ *INNER(N, H0, H2)*INNER(N, H1, H2)
A(2) = DSQRT(1./A(2))
DO I=0,N-1
LH0(I) = A(0)*H0(I)
ENDDO
DO I=0,N-1
LH1(I) = A(1)*H1(I)-A(1)*LH0(I)*INNER(N, LH0, H1)
ENDDO
DO I=0,N-1
LH2(I) = A(2)*H2(I)-A(2)*LH0(I)*INNER(N, LH0, H2)
$ -A(2)*LH1(I)*INNER(N, LH1, H2)
ENDDO
C(0) = INNER(N, F, LH0)
C(1) = INNER(N, F, LH1)
C(2) = INNER(N, F, LH2)
RP=c(1)*(a(0)*sqrt(2.))
IP=c(2)*(a(0)*sqrt(2.))
* write(*,*) "c(0)=", c(0)
END SUBROUTINE
REAL*8 FUNCTION INNER(N, A, B)
INTEGER N
REAL*8 A(0:N-1), B(0:N-1)
INNER = 0.
DO I=0,N-1
INNER = INNER+A(I)*B(I)
ENDDO
END FUNCTION
INTEGER FUNCTION MAXN(N, A)
INTEGER N
REAL*8 A(0:N-1), TEMP
MAXN=0
TEMP=A(0)
DO I=1,N-1
IF(A(I).GT.TEMP) THEN
MAXN=I
TEMP=A(I)
ENDIF
ENDDO
END FUNCTION
INTEGER FUNCTION MINN(N, A)
INTEGER N
REAL*8 A(0:N-1), TEMP
MINN=0
TEMP=A(0)
DO I=1,N-1
IF(A(I).LT.TEMP) THEN
MINN=I
TEMP=A(I)
ENDIF
ENDDO
END FUNCTION
