-
Notifications
You must be signed in to change notification settings - Fork 0
/
sinft.f90
49 lines (40 loc) · 894 Bytes
/
sinft.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
module m_sinft
contains
! SINFT
! taken from numerical recipes (see book for further information)
! subroutines called:
! NONE
SUBROUTINE SINFT(Y,N)
use m_realft
common /vocal/ ivocal
REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
real*8 Y(N)
real*8 Y1, Y2, SUM
THETA=3.14159265358979D0/DBLE(N)
WR=1.0D0
WI=0.0D0
WPR=-2.0D0*DSIN(0.5D0*THETA)**2
WPI=DSIN(THETA)
Y(1)=0.0
M=N/2
DO 11 J=1,M
WTEMP=WR
WR=WR*WPR-WI*WPI+WR
WI=WI*WPR+WTEMP*WPI+WI
Y1=WI*(Y(J+1)+Y(N-J+1))
Y2=0.5*(Y(J+1)-Y(N-J+1))
Y(J+1)=Y1+Y2
Y(N-J+1)=Y1-Y2
11 CONTINUE
CALL REALFT(Y,M,+1)
SUM=0.0
Y(1)=0.5*Y(1)
Y(2)=0.0
DO 12 J=1,N-1,2
SUM=SUM+Y(J)
Y(J)=Y(J+1)
Y(J+1)=SUM
12 CONTINUE
RETURN
END subroutine
end module m_sinft