-
Notifications
You must be signed in to change notification settings - Fork 2
/
four1.f
59 lines (55 loc) · 1.33 KB
/
four1.f
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
50
51
52
53
54
55
56
57
58
c FOUR1
c this routine is taken directly out of numerical recipes
c see the book for further information
c subroutines called:
c NONE
SUBROUTINE FOUR1(DATA,NN,ISIGN)
common /vocal/ ivocal
REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
DIMENSION DATA(*)
N=2*NN
J=1
DO 11 I=1,N,2
IF(J.GT.I)THEN
TEMPR=DATA(J)
TEMPI=DATA(J+1)
DATA(J)=DATA(I)
DATA(J+1)=DATA(I+1)
DATA(I)=TEMPR
DATA(I+1)=TEMPI
ENDIF
M=N/2
1 IF ((M.GE.2).AND.(J.GT.M)) THEN
J=J-M
M=M/2
GO TO 1
ENDIF
J=J+M
11 CONTINUE
MMAX=2
2 IF (N.GT.MMAX) THEN
ISTEP=2*MMAX
THETA=6.28318530717959D0/(ISIGN*MMAX)
WPR=-2.D0*DSIN(0.5D0*THETA)**2
WPI=DSIN(THETA)
WR=1.D0
WI=0.D0
DO 13 M=1,MMAX,2
DO 12 I=M,N,ISTEP
J=I+MMAX
TEMPR=SNGL(WR)*DATA(J)-SNGL(WI)*DATA(J+1)
TEMPI=SNGL(WR)*DATA(J+1)+SNGL(WI)*DATA(J)
DATA(J)=DATA(I)-TEMPR
DATA(J+1)=DATA(I+1)-TEMPI
DATA(I)=DATA(I)+TEMPR
DATA(I+1)=DATA(I+1)+TEMPI
12 CONTINUE
WTEMP=WR
WR=WR*WPR-WI*WPI+WR
WI=WI*WPR+WTEMP*WPI+WI
13 CONTINUE
MMAX=ISTEP
GO TO 2
ENDIF
RETURN
END