1      SUBROUTINE STCST (A, TCS, MODE, M, ND, MS, S)
2c>> 1997-03-31 STCST Krogh  Increased KEDIM, more sine table checks.
3c>> 1996-01-23 STCST Krogh  Changes to simplify conversion to C.
4C>> 1994-11-11 STCST Krogh  Declared all vars.
5c>> 1994-10-20 STCST Krogh  Changes to use M77CON
6c>> 1989-06-16 STCST FTK Fix error message on MODE, and TCS.
7c>> 1989-06-05 WVS Change length of MODE and TCS from (ND) to (*)
8c>> 1989-05-08 FTK & CLL
9c>> 1989-04-21 FTK & CLL
10c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
11c ALL RIGHTS RESERVED.
12c Based on Government Sponsored Research NAS7-03001.
13c
14c     This subroutine computes trigonometirc (sine-cosine), sine, or
15c     cosine transforms of real data in up to 6 dimensions using the
16c     Cooley-Tukey fast Fourier transform.
17c
18c     Variables in the calling sequence have the following types
19      REAL             A(*), S(*)
20      INTEGER   ND, M(ND), MS, KEDIM
21      CHARACTER*(*) TCS, MODE
22c
23c     Programmed by Fred F. Krogh at the Jet Propulsion Laboratory,
24c     Pasadena, Calif.   August 1, 1969.
25c     Revised for portability by Krogh -- January 29, 1988
26c
27c     Values for A, TCS, MODE, M, ND, and MS must be specified before
28c     calling the subroutine.
29c
30c     In describing the usage the following notation is used
31c     N(K) = 2 ** M(K)
32c     MA = M(1) + M(2) + ... + M(ND)
33c     NA = 2 ** MA
34c
35c     MTCS(K) = M(K)     TCS(K:K) = 'T'
36c             = M(K)+1   otherwise
37c
38c     MX = MAX(MTCS(1), MTCS(2), ..., MTCS(ND))
39c     NX = 2 ** MX
40c
41c     T(L,j,k) is defined differently for different values of TCS(L)
42c
43c       if TCS(L:L) = 'T' and MODE(L:L) = 'S', T(L,j,k)
44c         =1/2                     if k = 0
45c         =(1/2)*(-1)**j           if k = 1
46c         =COS(j*k*PI/N(L))        if k IS EVEN  (k = 2, 4, ..., N(L)-2)
47c         =SIN(j*(k-1)*PI/N(L))    if k IS ODD   (k = 3, 5, ..., N(L)-1)
48c                    and if MODE(L:L) = 'A', T(L,j,k)
49c         = (4/N) * (value of T(L,k,j) defined above)   If j<2
50c         = (2/N) * (value of T(L,k,j) defined above)   Otherwise
51c
52c       if TCS(L:L) = 'C' and MODE(L:L) = 'S', T(L,j,k)
53c         =1/2                     if k = 0
54c         =COS(j*k*PI/N(L))           k = 1, 2, ..., N(L)-1
55c         =(1/2)*(-1)**j           if k = N(L)
56c                    and if MODE(L:L) = 'A', T(L,j,k)
57c         = (2/N) * (value of T(L,j,k) defined above)
58c
59c       if TCS(L:L) = 'S' and MODE(L:L) = 'S', T(L,j,k)
60c         =SIN(j*k*PI/N(L))           k = 0, 1, ..., N(L)-1
61c                    and if MODE(L:L) = 'A', T(L,j,k)
62c         = (2/N) * (value of T(L,j,k) defined above)
63c
64c     D(L) = N(L)     if TCS(L:L) .ne. 'C'
65c          = N(L)+1   if TCS(L:L) = 'C'
66c
67c     The usage is as follows
68c
69c A() on input is an array of function values if one is doing Fourier
70c   analysis, and is an array of Fourier coefficients if one is doing
71c   Fourier synthesis.  On output, these are reversed.  In either case
72c   A() is a real array with dimension A(D(1), D(2), ..., D(ND)).
73c
74c TCS  is a character variable of length ND.  The k-th character must be
75c   'T' or 't' to select the general Trigonometric transform, or
76c   'C' or 'c' to select the Cosine transform, or
77c   'S' or 's' to select the Sine transform.
78c     See the description of T(L,j,k) and M above.
79c
80c MODE  A character variable of length ND.  The k-th character must be
81c   'A' or 'a' to select Analysis in the k-th dimension, or
82c   'S' or 's' to select Synthesis in the k-th dimension.
83c   One may be doing analysis, MODE(k:k) = 'A', with respect to one
84c   dimension and synthesis, MODE(k:k) = 'S', with respect to
85c   another.  A(j1+1, j2+1, ..., jND+1) is replaced by the sum over
86c   0 .le. k1 .le. D(1)-1, 0 .le. k2 .le. D(2)-1, ..., 0 .le. kND .le.
87c   D(ND)-1, of A(k1+1, k2+1, ..., kND+1) * T(1, k1, j1) * T(2, k2, j2)
88c   ... * T(ND, kND, jND), 0 .le. j1 .le. D(1)-1, ..., 0 .le. jND .le.
89c   D(ND)-1.
90
91c M() is a vector used to indicate N(k) = 2**M(k)).  The number of
92c   points in the k-th variable is then given by D(k) (see above).  M
93c   must be specified in such a way that 0 < M(k) < 22
94c   for k = 1, 2, ..., ND.
95c
96c ND is the dimension of (i.e. the number of subscripts in) the
97c    array A.  ND must satisfy 1 .le. ND .le. 6.
98c
99c MS gives the state of the sine table.  If MS > 0, there are NT =
100c    2 ** (MS-2) good entries in the sine table.  On the initial call,
101c    MS must be set to 0, and when the return is made, it will be set
102c    to MX, which is the value of MS required for computing a
103c    transform of size N.  If MS = -1, the sine table will be computed
104c    as for the case MS = 0, and then a return to the user will be made
105c    with MS set as before, but no transform will be computed.  This
106c    option is useful if the user would like access to the sine table
107c    before computing the FFT.
108c    On detected errors the error message subrs are called and
109c    execution stops.  If the user overrides the stop to cause
110c    continuation, then this subr will return with MS = -2.
111c
112c S() is a vector, S(j) = sin(pi*j/2*NT)), j = 1, 2, ..., NT-1, where
113c  NT is defined in the description of MS above.  S is computed by the
114c  subroutine if MX .gt. MS.  (If S is altered, set MS=0 so that S
115c  is recomputed.)
116c     -----------------------------------------------------------------
117c                Notes on COMMON, PARAMETER's, and local variables
118c
119c     NDMAX = the maximum value for ND, and MAXMX = the maximum
120c     permitted for MTCS(1), ..., MTCS(ND)
121c
122c     NF(1) = 1, NF(K+1) = NF(K) * D(K) K = 1, 2, ..., ND
123c
124c     MU is used in the process of eliminating transforms with respect
125c     to the first subscript of transforms with TCS(:) = 'S'.
126c     (This is only necessary if ND.GT.1.)
127c
128c     The dimension of KE must be at least as large as MAXMX-1.
129c     The named common CSFFTC is used for communication between this
130c     subroutine and the subroutine SFFT which computes a one
131c     dimensional complex Fourier transform and computes the sine table.
132c     The use of the variables in CSFFTC is contained in the listing
133c     of SFFT.
134c
135c     The input character variable TCS is mapped to the internal
136c     integer array ITCS() by mapping 'T' to 1, 'C' to 2, 'S' to 3.
137c
138c     -----------------------------------------------------------------
139c--S replaces "?": ?TCST, ?FFT, C?FFTC
140c     Both versions use IERM1
141c     and need ERFIN, IERV1
142c     -----------------------------------------------------------------
143      INTEGER MAXMX, NDMAX
144      PARAMETER (NDMAX = 6)
145
146      INTEGER I, I1, II, IR, ITCS(NDMAX), ITCSK
147      INTEGER J, JDIF, JJ, JK
148      INTEGER K, KDR
149      INTEGER KII, KIN, KK, KKI, KKL, KKN
150      INTEGER L
151      INTEGER MA, MI, MMAX, MSI, MU(NDMAX)
152      INTEGER N, NDD, NDIV, NF(NDMAX+1)
153      INTEGER NI, NI1, NI2, NI2I, NTOT2
154
155      CHARACTER  MSG1*13, MSG2*12
156
157      REAL             FN, FOUR, ONE
158      REAL             SPI4, SUM, T, T1, TP, TPI, TS, TS1, TWO
159      REAL             WI, WR, ZERO
160
161      PARAMETER (ZERO = 0.0E0)
162      PARAMETER (ONE = 1.0E0)
163      PARAMETER (TWO = 2.0E0)
164      PARAMETER (FOUR = 4.0E0)
165      PARAMETER (SPI4 = .70710 67811 86547 52440 08443 62104 8490E0)
166c Common variables
167      PARAMETER (KEDIM=30)
168      LOGICAL NEEDST
169      INTEGER MT, NT, MM, KS, ILAST, KE(KEDIM), KEE(KEDIM+1)
170c Note that KEE(1) is equivalent to ILAST.
171      EQUIVALENCE (KE(1), KEE(2))
172      COMMON /CSFFTC/ NEEDST, MT, NT, MM, KS, ILAST, KE
173      SAVE /CSFFTC/
174      PARAMETER (MAXMX = KEDIM+1)
175      DATA MSG1 / 'MODE(K:K) = ?' /
176      DATA MSG2 / 'TCS(K:K) = ?' /
177c     -----------------------------------------------------------------
178c
179      NDD = ND
180      IF ((NDD .LE. 0) .OR. (NDD .GT. NDMAX)) THEN
181C                               FATAL ERROR, DEFAULT IS TO STOP IN IERM1
182         CALL IERM1 ('STCST', 1, 2, 'BAD ND', 'ND', ND, '.')
183         MS = -2
184         RETURN
185      END IF
186      MA = 0
187      MMAX = 0
188      NDIV = 1
189c Every element in the array A is divided by NDIV before computing
190c the transform.  The value computed for NDIV depends on whether
191c one is doing analysis or synthesis and on the type of
192c transform being computed.
193      DO 10 K = 1, NDD
194         MM = M(K)
195         IF ((MM .LT. 0) .OR. (MM .GT. MAXMX)) GO TO 200
196         MA = MA + MM
197         N = 2 ** MM
198         IF (MODE(K:K) .eq. 'A' .or. MODE(K:K) .eq. 'a') THEN
199            NDIV = NDIV * N
200         ELSE IF (MODE(K:K) .eq. 'S' .or. MODE(K:K) .eq. 's') THEN
201            NDIV = 2 * NDIV
202         else
203            MSG1(13:13) = MODE(K:K)
204            CALL IERM1 ('STCST',2,2, MSG1, 'for K =',  K, '.')
205            MS = -2
206            return
207         END IF
208         ITCSK = (index('TtCcSs', TCS(K:K)) + 1) / 2
209         ITCS(K) = ITCSK
210         if(ITCSK .eq. 0) then
211            MSG2(12:12) = TCS(K:K)
212            CALL IERM1 ('STCST', 3, 2, MSG2, 'for K =', K, '.')
213            return
214         endif
215         NF(1) = 1
216         IF (ITCSK .GE. 2) THEN
217            IF (ITCSK .EQ. 2) N = N + 1
218            NDIV = 2 * NDIV
219            MM = MM + 1
220         END IF
221         NF(K+1) = NF(K) * N
222         IF (MM .GT. MMAX) THEN
223            MMAX = MM
224         END IF
225   10 CONTINUE
226c
227      MSI = MS
228      NEEDST = MMAX .GT. MSI
229
230      if (.NOT. NEEDST) then
231c  Check internal parameters to catch certain user errors.
232         if (MT .lt. KEDIM) then
233            if (MMAX .le. MT + 2) then
234c              Skip sine table computation if all appears O.K.
235               if (MT .le. 0) go to 15
236               if (abs(S(NT/2) - SPI4) .le. 1.E-7) go to 15
237            end if
238         end if
239         NEEDST = .true.
240         call ERMSG('STCST', 3, 1,
241     1      'Invalid sine table (re)computed', '.')
242      end if
243      MS = MMAX
244      MT = MMAX - 2
245      CALL SFFT (A, A, S)
246      IF (MSI .EQ. -1) RETURN
247c                   All setup for computation now
248   15 NTOT2 = NF(NDD+1)
249c
250      FN = ONE / REAL(NDIV)
251c     Divide every element of A by NDIV
252      DO 20 I = 1, NTOT2
253         A(I) = A(I) * FN
254   20 CONTINUE
255c
256c     Beginning of loop for computing multiple sum
257      DO 170 K = 1, NDD
258         ITCSK = ITCS(K)
259         MI = M(K)
260         MM = MI - 1
261         if(MODE(K:K) .eq. 'A' .or. MODE(K:K) .eq. 'a') MI = -MI
262         KDR = NF(K)
263         KS = KDR + KDR
264         ILAST = NF(K+1)
265         IF (ITCSK .EQ. 2) ILAST = ILAST - KDR
266         DO 30 L = 1, MM
267            KEE(L+1) = KEE(L) / 2
268   30    CONTINUE
269c
270         I = 1
271         J = NDD
272   40    DO 50 L = 1, J
273            MU(L) = 0
274            IF ((L .NE. K) .AND. (ITCS(L) .GT. 2)) THEN
275c           Skip the part of the array left empty by the sine transform
276               MU(L) = NF(L)
277               I = I + NF(L)
278            END IF
279   50    CONTINUE
280c
281c        Compute indices associated with the current value of I (and K)
282   60    I1 = I + KDR
283         NI1 = I + NF(K+1)
284         IF (ITCSK .EQ. 2) NI1 = NI1 - KDR
285         NI = NI1 - KDR
286         NI2 = (NI1 + I) / 2
287         NI2I = NI2 + KDR
288         IF (ITCSK .NE. 1) THEN
289c                Doing a cosine or a sine transform -- set MI = 0 and do
290c                calculations not required for sine-cosine transforms
291            MI = 0
292            J = NI
293            SUM = A(I1)
294            T = A(J)
295   70       JK = J - KS
296               IF (JK .GE.I1) THEN
297                  SUM = SUM + A(J)
298                  A(J) = A(JK) - A(J)
299                  J = JK
300                  GO TO 70
301               END IF
302            IF (ITCSK .NE. 2) THEN
303c                                Calculations for the sine transform
304               A(I) = TWO * A(I1)
305               A(I1) = -TWO * T
306               IF (MM.EQ.0) GO TO 90
307               T = TWO * A(NI2)
308               A(NI2) = -TWO * A(NI2I)
309               A(NI2I) = T
310               GO TO 90
311            END IF
312c                               Set for cosine transform
313            A(I1) = A(NI1)
314         END IF
315         IF (MM.EQ.0) GO TO 90
316         if (MI .lt. 0) go to 160
317c        Begin calculations for the sine-cosine transform
318   80    A(NI2) = TWO * A(NI2)
319         A(NI2I) = TWO * A(NI2I)
320   90    T = A(I) + A(I1)
321         A(I1) = A(I) - A(I1)
322         IF (MI .LT. 0) THEN
323            IF (ITCSK .EQ. 1) THEN
324               A(I1) = TWO * A(I1)
325               T = TWO * T
326            END IF
327         END IF
328         A(I) = T
329         J = 0
330         JDIF = 2 ** (MT - MM + 1)
331         KKL = KE(1) - KDR
332         IF (MM .GT. 1) THEN
333            DO 100 KK = KS, KKL, KS
334               KKI = I + KK
335               KII = KKI + KDR
336               KKN = NI1 - KK
337               KIN = KKN + KDR
338               J = J + JDIF
339               WI = S(J)
340               JJ = NT - J
341               WR = S(JJ)
342               T = A(KKI) + A(KKN)
343               TS = A(KKN) - A(KKI)
344               T1 = A(KIN) - A(KII)
345               TS1 = A(KIN) + A(KII)
346               IF (ITCSK .GT. 2) THEN
347c                         The sine-cosine transform must be computed
348c                         differently in the case of the sine transform
349c                         because the input data is stored differently.
350                  TP = WR * T - WI * T1
351                  TPI = WI * T + WR * T1
352                  A(KKI) = TP - TS1
353                  A(KKN) = -TP - TS1
354                  A(KII) = TPI + TS
355                  A(KIN) = TPI - TS
356               ELSE
357                  TP = WR * TS1 + WI * TS
358                  TPI = WI * TS1 - WR * TS
359                  A(KKI) = T + TP
360                  A(KKN) = T - TP
361                  A(KII) = T1 + TPI
362                  A(KIN) = TPI - T1
363               END IF
364  100       CONTINUE
365         ELSE IF (MM .EQ. 0) THEN
366            GO TO 120
367         END IF
368c        End of computing sine-cosine transform
369c
370         IF (MI .LT. 0) GO TO 140
371         IR = I
372         II = IR + KDR
373c
374c        Compute a one-dimensional complex Fourier transform
375  110    CALL SFFT (A(IR), A(II), S)
376         IF (MI .LT. 0) GO TO 80
377         IF (MI .NE. 0) GO TO 140
378  120    IF (ITCSK .EQ. 1) GO TO 140
379         IF (ITCSK .EQ. 3) THEN
380            A(I) = ZERO
381         ELSE
382c                Compute first and last elements of the cosine transform
383            SUM = FOUR * SUM
384            T = TWO * A(I)
385            A(NI1) = T - SUM
386            A(I) = T + SUM
387            IF (MM.GE.0) A(NI2) = TWO * A(NI2)
388         END IF
389         IF (MM .GT. 0) THEN
390c               Extra calculations required by sine and cosine transform
391            J = 0
392            JDIF = JDIF / 2
393            DO 130 KK = KDR, KKL, KDR
394               KKI = I + KK
395               KKN = NI1 - KK
396               J = J + JDIF
397               WI = TWO * S(J)
398               T = A(KKI) + A(KKN)
399               TS = A(KKI) - A(KKN)
400               IF (ITCSK .NE. 2) THEN
401                  T = T / WI
402               ELSE
403                  TS = TS / WI
404               END IF
405               A(KKI) = T + TS
406               A(KKN) = T - TS
407  130       CONTINUE
408         END IF
409c
410c        Logic for deciding which one-dimensional transform to do next
411  140    J = 0
412  150    J = J + 1
413            IF (J .EQ. K) THEN
414               J = J + 1
415               I = I + NF(J) - NF(J-1)
416            END IF
417            IF (J.GT.NDD) GO TO 170
418            MU(J) = MU(J) + NF(J)
419            IF (MU(J).GE.NF(J+1)) GO TO 150
420         I = I + NF(1)
421         J = J - 1
422         IF (J .NE. 0) GO TO 40
423         GO TO 60
424c
425c        Set for Fourier analysis
426  160    II = I
427         IR = II + KDR
428         GO TO 110
429c
430  170 CONTINUE
431c     End of K loop
432      RETURN
433c
434c                               Fatal error, default is to stop in IERM1
435  200 CALL IERM1 ('STCST', 4, 2, 'Require 0 .le. max(M(K)) .le. 31',
436     1  'M', MM, '.')
437      MS = -2
438      RETURN
439c
440      END
441