1C
2      SUBROUTINE FFTRC  (A,N,X,IWK,WK)
3C                                  SPECIFICATIONS FOR ARGUMENTS
4      INTEGER            N,IWK(1)
5      REAL*8             A(N),WK(1)
6      COMPLEX*16         X(1)
7C                                  SPECIFICATIONS FOR LOCAL VARIABLES
8      INTEGER            ND2P1,ND2,I,MTWO,M,IMAX,ND4,NP2,K,NMK,J
9      REAL*8             RPI,ZERO,ONE,HALF,THETA,TP,G(2),B(2),Z(2),AI,
10     1                   AR
11      COMPLEX*16         XIMAG,ALPH,BETA,GAM,S1,ZD
12      EQUIVALENCE        (GAM,G(1)),(ALPH,B(1)),(Z(1),AR),(Z(2),AI),
13     1                   (ZD,Z(1))
14      DATA               ZERO/0.0D0/,HALF/0.5D0/,ONE/1.0D0/,IMAX/24/
15      DATA               RPI/3.141592653589793D0/
16C                                  FIRST EXECUTABLE STATEMENT
17      IF (N .NE. 2) GO TO 5
18C                                  N EQUAL TO 2
19      ZD = DCMPLX(A(1),A(2))
20      THETA = AR
21      TP = AI
22      X(2) = DCMPLX(THETA-TP,ZERO)
23      X(1) = DCMPLX(THETA+TP,ZERO)
24      GO TO 9005
25    5 CONTINUE
26C                                  N GREATER THAN 2
27      ND2 = N/2
28      ND2P1 = ND2+1
29C                                  MOVE A TO X
30      J = 1
31      DO 6 I=1,ND2
32         X(I) = DCMPLX(A(J),A(J+1))
33         J = J+2
34    6 CONTINUE
35C                                  COMPUTE THE CENTER COEFFICIENT
36      GAM = DCMPLX(ZERO,ZERO)
37      DO 10 I=1,ND2
38         GAM = GAM + X(I)
39   10 CONTINUE
40      TP = G(1)-G(2)
41      GAM = DCMPLX(TP,ZERO)
42C                                  DETERMINE THE SMALLEST M SUCH THAT
43C                                  N IS LESS THAN OR EQUAL TO 2**M
44      MTWO = 2
45      M = 1
46      DO 15 I=1,IMAX
47         IF (ND2 .LE. MTWO) GO TO 20
48         MTWO = MTWO+MTWO
49         M = M+1
50   15 CONTINUE
51   20 IF (ND2 .EQ. MTWO) GO TO 25
52C                                  N IS NOT A POWER OF TWO, CALL FFTCC
53      CALL FFTCC (X,ND2,IWK,WK)
54      GO TO 30
55C                                  N IS A POWER OF TWO, CALL FFT2C
56   25 CALL FFT2C (X,M,IWK)
57   30 ALPH = X(1)
58      X(1) = B(1) + B(2)
59      ND4 = (ND2+1)/2
60      IF (ND4 .LT. 2) GO TO 40
61      NP2 = ND2 + 2
62      THETA = RPI/ND2
63      TP = THETA
64      XIMAG = DCMPLX(ZERO,ONE)
65C                                  DECOMPOSE THE COMPLEX VECTOR X
66C                                  INTO THE COMPONENTS OF THE TRANSFORM
67C                                  OF THE INPUT DATA.
68      DO 35 K = 2,ND4
69         NMK = NP2 - K
70         S1 = DCONJG(X(NMK))
71         ALPH = X(K) + S1
72         BETA = XIMAG*(S1-X(K))
73         S1 = DCMPLX(DCOS(THETA),DSIN(THETA))
74         X(K) = (ALPH+BETA*S1)*HALF
75         X(NMK) = DCONJG(ALPH-BETA*S1)*HALF
76         THETA = THETA + TP
77   35 CONTINUE
78   40 CONTINUE
79      X(ND2P1) = GAM
80 9005 RETURN
81      END
82
83