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