/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:48 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "stcst.h"
#include <stdlib.h>
#include <string.h>
/* PARAMETER translations */
#define FOUR 4.0e0
#define KEDIM 30
#define MAXMX (KEDIM + 1)
#define NDMAX 6
#define ONE 1.0e0
#define SPI4 .7071067811865475244008443621048490e0
#define TWO 2.0e0
#define ZERO 0.0e0
/* end of PARAMETER translations */
/* COMMON translations */
struct t_csfftc {
LOGICAL32 needst;
long int mt, nt, mm, ks, ilast, ke[KEDIM];
} csfftc;
/* end of COMMON translations */
void /*FUNCTION*/ stcst(
float a[],
char *tcs,
char *mode,
long m[],
long nd,
long *ms,
float s[])
{
char _c0[2];
long int _d_l, _d_m, _do0, _do1, _do2, _do3, i, i1, ii, ir, itcs[NDMAX],
itcsk, j, jdif, jj, jk, k, kdr, kii, kin, kk, kki, kkl, kkn,
l, ma, mi, mmax, msi, mu[NDMAX], n, ndd, ndiv, nf[NDMAX + 1],
ni, ni1, ni2, ni2i, ntot2;
float fn, sum, t, t1, tp, tpi, ts, ts1, wi, wr;
static char msg1[14] = "MODE(K:K) = ?";
static char msg2[13] = "TCS(K:K) = ?";
long int *const kee = (long*)((long*)csfftc.ke + -1);
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const A = &a[0] - 1;
long *const Itcs = &itcs[0] - 1;
long *const Ke = &csfftc.ke[0] - 1;
long *const Kee = &kee[0] - 1;
long *const M = &m[0] - 1;
long *const Mu = &mu[0] - 1;
long *const Nf = &nf[0] - 1;
float *const S = &s[0] - 1;
/* end of OFFSET VECTORS */
/*>> 1997-03-31 STCST Krogh Increased KEDIM, more sine table checks.
*>> 1996-01-23 STCST Krogh Changes to simplify conversion to C.
*>> 1994-11-11 STCST Krogh Declared all vars.
*>> 1994-10-20 STCST Krogh Changes to use M77CON
*>> 1989-06-16 STCST FTK Fix error message on MODE, and TCS.
*>> 1989-06-05 WVS Change length of MODE and TCS from (ND) to (*)
*>> 1989-05-08 FTK & CLL
*>> 1989-04-21 FTK & CLL
* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
* ALL RIGHTS RESERVED.
* Based on Government Sponsored Research NAS7-03001.
*
* This subroutine computes trigonometirc (sine-cosine), sine, or
* cosine transforms of real data in up to 6 dimensions using the
* Cooley-Tukey fast Fourier transform.
*
* Variables in the calling sequence have the following types */
/* Programmed by Fred F. Krogh at the Jet Propulsion Laboratory,
* Pasadena, Calif. August 1, 1969.
* Revised for portability by Krogh -- January 29, 1988
*
* Values for A, TCS, MODE, M, ND, and MS must be specified before
* calling the subroutine.
*
* In describing the usage the following notation is used
* N(K) = 2 ** M(K)
* MA = M(1) + M(2) + ... + M(ND)
* NA = 2 ** MA
*
* MTCS(K) = M(K) TCS(K:K) = 'T'
* = M(K)+1 otherwise
*
* MX = MAX(MTCS(1), MTCS(2), ..., MTCS(ND))
* NX = 2 ** MX
*
* T(L,j,k) is defined differently for different values of TCS(L)
*
* if TCS(L:L) = 'T' and MODE(L:L) = 'S', T(L,j,k)
* =1/2 if k = 0
* =(1/2)*(-1)**j if k = 1
* =COS(j*k*PI/N(L)) if k IS EVEN (k = 2, 4, ..., N(L)-2)
* =SIN(j*(k-1)*PI/N(L)) if k IS ODD (k = 3, 5, ..., N(L)-1)
* and if MODE(L:L) = 'A', T(L,j,k)
* = (4/N) * (value of T(L,k,j) defined above) If j<2
* = (2/N) * (value of T(L,k,j) defined above) Otherwise
*
* if TCS(L:L) = 'C' and MODE(L:L) = 'S', T(L,j,k)
* =1/2 if k = 0
* =COS(j*k*PI/N(L)) k = 1, 2, ..., N(L)-1
* =(1/2)*(-1)**j if k = N(L)
* and if MODE(L:L) = 'A', T(L,j,k)
* = (2/N) * (value of T(L,j,k) defined above)
*
* if TCS(L:L) = 'S' and MODE(L:L) = 'S', T(L,j,k)
* =SIN(j*k*PI/N(L)) k = 0, 1, ..., N(L)-1
* and if MODE(L:L) = 'A', T(L,j,k)
* = (2/N) * (value of T(L,j,k) defined above)
*
* D(L) = N(L) if TCS(L:L) .ne. 'C'
* = N(L)+1 if TCS(L:L) = 'C'
*
* The usage is as follows
*
* A() on input is an array of function values if one is doing Fourier
* analysis, and is an array of Fourier coefficients if one is doing
* Fourier synthesis. On output, these are reversed. In either case
* A() is a real array with dimension A(D(1), D(2), ..., D(ND)).
*
* TCS is a character variable of length ND. The k-th character must be
* 'T' or 't' to select the general Trigonometric transform, or
* 'C' or 'c' to select the Cosine transform, or
* 'S' or 's' to select the Sine transform.
* See the description of T(L,j,k) and M above.
*
* MODE A character variable of length ND. The k-th character must be
* 'A' or 'a' to select Analysis in the k-th dimension, or
* 'S' or 's' to select Synthesis in the k-th dimension.
* One may be doing analysis, MODE(k:k) = 'A', with respect to one
* dimension and synthesis, MODE(k:k) = 'S', with respect to
* another. A(j1+1, j2+1, ..., jND+1) is replaced by the sum over
* 0 .le. k1 .le. D(1)-1, 0 .le. k2 .le. D(2)-1, ..., 0 .le. kND .le.
* D(ND)-1, of A(k1+1, k2+1, ..., kND+1) * T(1, k1, j1) * T(2, k2, j2)
* ... * T(ND, kND, jND), 0 .le. j1 .le. D(1)-1, ..., 0 .le. jND .le.
* D(ND)-1. */
/* M() is a vector used to indicate N(k) = 2**M(k)). The number of
* points in the k-th variable is then given by D(k) (see above). M
* must be specified in such a way that 0 < M(k) < 22
* for k = 1, 2, ..., ND.
*
* ND is the dimension of (i.e. the number of subscripts in) the
* array A. ND must satisfy 1 .le. ND .le. 6.
*
* MS gives the state of the sine table. If MS > 0, there are NT =
* 2 ** (MS-2) good entries in the sine table. On the initial call,
* MS must be set to 0, and when the return is made, it will be set
* to MX, which is the value of MS required for computing a
* transform of size N. If MS = -1, the sine table will be computed
* as for the case MS = 0, and then a return to the user will be made
* with MS set as before, but no transform will be computed. This
* option is useful if the user would like access to the sine table
* before computing the FFT.
* On detected errors the error message subrs are called and
* execution stops. If the user overrides the stop to cause
* continuation, then this subr will return with MS = -2.
*
* S() is a vector, S(j) = sin(pi*j/2*NT)), j = 1, 2, ..., NT-1, where
* NT is defined in the description of MS above. S is computed by the
* subroutine if MX .gt. MS. (If S is altered, set MS=0 so that S
* is recomputed.)
* -----------------------------------------------------------------
* Notes on COMMON, PARAMETER's, and local variables
*
* NDMAX = the maximum value for ND, and MAXMX = the maximum
* permitted for MTCS(1), ..., MTCS(ND)
* */
/* NF(1) = 1, NF(K+1) = NF(K) * D(K) K = 1, 2, ..., ND
*
* MU is used in the process of eliminating transforms with respect
* to the first subscript of transforms with TCS(:) = 'S'.
* (This is only necessary if ND.GT.1.)
*
* The dimension of KE must be at least as large as MAXMX-1.
* The named common CSFFTC is used for communication between this
* subroutine and the subroutine SFFT which computes a one
* dimensional complex Fourier transform and computes the sine table.
* The use of the variables in CSFFTC is contained in the listing
* of SFFT.
*
* The input character variable TCS is mapped to the internal
* integer array ITCS() by mapping 'T' to 1, 'C' to 2, 'S' to 3.
*
* -----------------------------------------------------------------
*--S replaces "?": ?TCST, ?FFT, C?FFTC
* Both versions use IERM1
* and need ERFIN, IERV1
* ----------------------------------------------------------------- */
/* Common variables */
/* Note that KEE(1) is equivalent to ILAST. */
/* -----------------------------------------------------------------
* */
ndd = nd;
if ((ndd <= 0) || (ndd > NDMAX))
{
/* FATAL ERROR, DEFAULT IS TO STOP IN IERM1 */
ierm1( "STCST", 1, 2, "BAD ND", "ND", nd, '.' );
*ms = -2;
return;
}
ma = 0;
mmax = 0;
ndiv = 1;
/* Every element in the array A is divided by NDIV before computing
* the transform. The value computed for NDIV depends on whether
* one is doing analysis or synthesis and on the type of
* transform being computed. */
for (k = 1; k <= ndd; k++)
{
csfftc.mm = M[k];
if ((csfftc.mm < 0) || (csfftc.mm > MAXMX))
goto L_200;
ma += csfftc.mm;
n = ipow(2,csfftc.mm);
if (mode[k - 1] == 'A' || mode[k - 1] == 'a')
{
ndiv *= n;
}
else if (mode[k - 1] == 'S' || mode[k - 1] == 's')
{
ndiv *= 2;
}
else
{
msg1[12] = mode[k - 1];
ierm1( "STCST", 2, 2, msg1, "for K =", k, '.' );
*ms = -2;
return;
}
itcsk = (istrstr( "TtCcSs", STR1(_c0,tcs[k - 1]) ) + 1)/2;
Itcs[k] = itcsk;
if (itcsk == 0)
{
msg2[11] = tcs[k - 1];
ierm1( "STCST", 3, 2, msg2, "for K =", k, '.' );
return;
}
Nf[1] = 1;
if (itcsk >= 2)
{
if (itcsk == 2)
n += 1;
ndiv *= 2;
csfftc.mm += 1;
}
Nf[k + 1] = Nf[k]*n;
if (csfftc.mm > mmax)
{
mmax = csfftc.mm;
}
}
msi = *ms;
csfftc.needst = mmax > msi;
if (!csfftc.needst)
{
/* Check internal parameters to catch certain user errors. */
if (csfftc.mt < KEDIM)
{
if (mmax <= csfftc.mt + 2)
{
/* Skip sine table computation if all appears O.K. */
if (csfftc.mt <= 0)
goto L_15;
if (fabsf( S[csfftc.nt/2] - SPI4 ) <= 1.e-7)
goto L_15;
}
}
csfftc.needst = TRUE;
ermsg( "STCST", 3, 1, "Invalid sine table (re)computed", '.' );
}
*ms = mmax;
csfftc.mt = mmax - 2;
sfft( a, a, s );
if (msi == -1)
return;
/* All setup for computation now */
L_15:
ntot2 = Nf[ndd + 1];
fn = ONE/(float)( ndiv );
/* Divide every element of A by NDIV */
for (i = 1; i <= ntot2; i++)
{
A[i] *= fn;
}
/* Beginning of loop for computing multiple sum */
for (k = 1; k <= ndd; k++)
{
itcsk = Itcs[k];
mi = M[k];
csfftc.mm = mi - 1;
if (mode[k - 1] == 'A' || mode[k - 1] == 'a')
mi = -mi;
kdr = Nf[k];
csfftc.ks = kdr + kdr;
csfftc.ilast = Nf[k + 1];
if (itcsk == 2)
csfftc.ilast -= kdr;
for (l = 1; l <= csfftc.mm; l++)
{
Kee[l + 1] = Kee[l]/2;
}
i = 1;
j = ndd;
L_40:
for (l = 1; l <= j; l++)
{
Mu[l] = 0;
if ((l != k) && (Itcs[l] > 2))
{
/* Skip the part of the array left empty by the sine transform */
Mu[l] = Nf[l];
i += Nf[l];
}
}
/* Compute indices associated with the current value of I (and K) */
L_60:
i1 = i + kdr;
ni1 = i + Nf[k + 1];
if (itcsk == 2)
ni1 -= kdr;
ni = ni1 - kdr;
ni2 = (ni1 + i)/2;
ni2i = ni2 + kdr;
if (itcsk != 1)
{
/* Doing a cosine or a sine transform -- set MI = 0 and do
* calculations not required for sine-cosine transforms */
mi = 0;
j = ni;
sum = A[i1];
t = A[j];
L_70:
jk = j - csfftc.ks;
if (jk >= i1)
{
sum += A[j];
A[j] = A[jk] - A[j];
j = jk;
goto L_70;
}
if (itcsk != 2)
{
/* Calculations for the sine transform */
A[i] = TWO*A[i1];
A[i1] = -TWO*t;
if (csfftc.mm == 0)
goto L_90;
t = TWO*A[ni2];
A[ni2] = -TWO*A[ni2i];
A[ni2i] = t;
goto L_90;
}
/* Set for cosine transform */
A[i1] = A[ni1];
}
if (csfftc.mm == 0)
goto L_90;
if (mi < 0)
goto L_160;
/* Begin calculations for the sine-cosine transform */
L_80:
A[ni2] *= TWO;
A[ni2i] *= TWO;
L_90:
t = A[i] + A[i1];
A[i1] = A[i] - A[i1];
if (mi < 0)
{
if (itcsk == 1)
{
A[i1] *= TWO;
t *= TWO;
}
}
A[i] = t;
j = 0;
jdif = ipow(2,csfftc.mt - csfftc.mm + 1);
kkl = Ke[1] - kdr;
if (csfftc.mm > 1)
{
for (kk = csfftc.ks, _do0=DOCNT(kk,kkl,_do1 = csfftc.ks); _do0 > 0; kk += _do1, _do0--)
{
kki = i + kk;
kii = kki + kdr;
kkn = ni1 - kk;
kin = kkn + kdr;
j += jdif;
wi = S[j];
jj = csfftc.nt - j;
wr = S[jj];
t = A[kki] + A[kkn];
ts = A[kkn] - A[kki];
t1 = A[kin] - A[kii];
ts1 = A[kin] + A[kii];
if (itcsk > 2)
{
/* The sine-cosine transform must be computed
* differently in the case of the sine transform
* because the input data is stored differently. */
tp = wr*t - wi*t1;
tpi = wi*t + wr*t1;
A[kki] = tp - ts1;
A[kkn] = -tp - ts1;
A[kii] = tpi + ts;
A[kin] = tpi - ts;
}
else
{
tp = wr*ts1 + wi*ts;
tpi = wi*ts1 - wr*ts;
A[kki] = t + tp;
A[kkn] = t - tp;
A[kii] = t1 + tpi;
A[kin] = tpi - t1;
}
}
}
else if (csfftc.mm == 0)
{
goto L_120;
}
/* End of computing sine-cosine transform
* */
if (mi < 0)
goto L_140;
ir = i;
ii = ir + kdr;
/* Compute a one-dimensional complex Fourier transform */
L_110:
sfft( &A[ir], &A[ii], s );
if (mi < 0)
goto L_80;
if (mi != 0)
goto L_140;
L_120:
if (itcsk == 1)
goto L_140;
if (itcsk == 3)
{
A[i] = ZERO;
}
else
{
/* Compute first and last elements of the cosine transform */
sum *= FOUR;
t = TWO*A[i];
A[ni1] = t - sum;
A[i] = t + sum;
if (csfftc.mm >= 0)
A[ni2] *= TWO;
}
if (csfftc.mm > 0)
{
/* Extra calculations required by sine and cosine transform */
j = 0;
jdif /= 2;
for (kk = kdr, _do2=DOCNT(kk,kkl,_do3 = kdr); _do2 > 0; kk += _do3, _do2--)
{
kki = i + kk;
kkn = ni1 - kk;
j += jdif;
wi = TWO*S[j];
t = A[kki] + A[kkn];
ts = A[kki] - A[kkn];
if (itcsk != 2)
{
t /= wi;
}
else
{
ts /= wi;
}
A[kki] = t + ts;
A[kkn] = t - ts;
}
}
/* Logic for deciding which one-dimensional transform to do next */
L_140:
j = 0;
L_150:
j += 1;
if (j == k)
{
j += 1;
i += Nf[j] - Nf[j - 1];
}
if (j > ndd)
goto L_170;
Mu[j] += Nf[j];
if (Mu[j] >= Nf[j + 1])
goto L_150;
i += Nf[1];
j -= 1;
if (j != 0)
goto L_40;
goto L_60;
/* Set for Fourier analysis */
L_160:
ii = i;
ir = ii + kdr;
goto L_110;
L_170:
;
}
/* End of K loop */
return;
/* Fatal error, default is to stop in IERM1 */
L_200:
ierm1( "STCST", 4, 2, "Require 0 .le. max(M(K)) .le. 31", "M",
csfftc.mm, '.' );
*ms = -2;
return;
} /* end of function */