/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:45 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include
#include "fcrt.h"
#include "scft.h"
#include
#include
/* PARAMETER translations */
#define KEDIM 30
#define MAXM KEDIM
#define NDMAX 6
#define ONE 1.e0
#define SPI4 .7071067811865475244008443621048490e0
/* 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*/ scft(
float a[],
char *mode,
long m[],
long nd,
long *ms,
float s[])
{
long int _d_l, _d_m, _do0, _do1, i, ii, iic, ir, irc, j, jl, k,
ksm, l, ma, mmax, msi, ndd, nf[NDMAX + 1], ntot2;
float fn;
static char cerror[14] = "MODE(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 Kee = &kee[0] - 1;
long *const M = &m[0] - 1;
long *const Nf = &nf[0] - 1;
float *const S = &s[0] - 1;
/* end of OFFSET VECTORS */
/*>> 1997-03-31 SCFT Krogh Increased KEDIM, more sine table checks.
*>> 1994-10-20 SCFT Krogh Change CERROR to simplify conversion to C.
*>> 1994-10-20 SCFT Krogh Changes to use M77CON
*>> 1994-04-19 SCFT CLL Edited to make DP & SP files similar.
*>> 1989-06-16 FTK Fix error message on MODE.
*>> 1989-06-05 WVS Change length of MODE from *(ND) to *(*).
*>> 1989-05-07 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 complex Fourier transforms in up to 6
* dimensions using the Cooley-Tukey fast Fourier transform.
*
* Variables in the calling sequence have the following types. */
/* Programmed by Fred T. Krogh at the Jet Propulsion Laboratory,
* Pasadena, Calif. August 1, 1969.
*
* Values for A, 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)
* WK = EXP(2*PI*I/N(K)), where I = sqrt(-1) and PI = 3.14159...
* MA = M(1) + M(2) + ... + M(ND)
* NA = 2**MA
* MX = MAX(M(1), M(2), ..., M(ND))
* NX = 2**MX
*
* The usage is as follows
*
* A() on input is an array of function values, X, if one is doing
* Fourier analysis,
* and is an array of Fourier coefficients, C, if one is
* doing Fourier synthesis. On output, A contains C if doing analysis,
* and contains X if doing synthesis. In our description here, we
* assume that A has been declared by the user as a complex array with
* dimension A(N(1),N(2),...,N(ND)). (Of course this is not necessary
* as long as one declares enough storage for A, and keeps in mind that
* the array is treated in this subroutine as if it had been specified
* in this way.) It is assumed that the imaginary part of a complex
* number is stored in the cell immediately following its real part.
*
* MODE A character variable of length at least 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.
* The subroutine sets A(j1+1, j2+1, ..., jND+1) to the sum
* over 0 .le. k1 .le. N(1)-1, 0 .le. k2 .le. N(2)-1, ...,
* 0 .le. kND .le. N(ND)-1 of A(k1+1, k2+1, ..., kND+1) * T(1,j1,k1) *
* T(2,j2,k2) * ... * T(ND,jND,kND), 0 .le. j1 .le. N(1)-1, ...,
* 0 .le. jnd .le. N(ND)-1
* with T(L,j,k) =
* WL ** (j*k) if MODE(L:L) = 'S' or 's'
* (1/N(L)) * WL ** (-j*k) if MODE(L:L) = 'A' or 'a'
*
* M() is a vector used to indicate N(k) = 2**M(k), the number of
* complex points in the k-th variable. M must be such that MX < 21.
* If M(L)=0, no change is made with respect to the L-th dimension.
*
* 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 NX. 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.
*
* 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 MAXM = the maximum
* permitted for M(1), ..., M(ND)
* The dimension of KE must be at least as large as MAXM.
* 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.
*
* NF(1) = 2, NF(K+1) = 2**(1+M(1)+...+M(K)), K = 1,..., ND
* ------------------------------------------------------------------
*--S replaces "?": ?CFT, ?FFT, C?FFTC
* Both versions use IERM1
* ------------------------------------------------------------------ */
/* 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( "SCFT", 1, 2, "BAD ND", "ND", nd, '.' );
*ms = -2;
return;
}
Nf[1] = 2;
ma = 0;
mmax = 0;
fn = ONE;
for (k = 1; k <= ndd; k++)
{
csfftc.mm = M[k];
mmax = max( csfftc.mm, mmax );
ma += csfftc.mm;
i = ipow(2,csfftc.mm);
if ((csfftc.mm < 0) || (csfftc.mm > MAXM))
{
/* Fatal error, default is to stop in IERM1 */
ierm1( "SCFT", 2, 2, "Need 0 .le. M(K) .le. 30", "M(K)"
, csfftc.mm, '.' );
*ms = -2;
return;
}
Nf[k + 1] = Nf[k]*i;
if (mode[k - 1] == 'A' || mode[k - 1] == 'a')
{
fn /= i;
}
else if (mode[k - 1] != 'S' && mode[k - 1] != 's')
{
cerror[12] = mode[k - 1];
ierm1( "SCFT", 3, 2, cerror, "for K =", k, '.' );
*ms = -2;
return;
}
}
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_20;
if (fabsf( S[csfftc.nt/2] - SPI4 ) <= 1.e-7)
goto L_20;
}
}
csfftc.needst = TRUE;
ermsg( "DRFT1", 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_20:
ntot2 = Nf[ndd + 1];
if (fn != ONE)
{
for (i = 1; i <= ntot2; i++)
{
A[i] *= fn;
}
}
/* Beginning of loop for computing multiple sum */
for (k = 1; k <= ndd; k++)
{
if (mode[k - 1] == 'A' || mode[k - 1] == 'a')
{
irc = 1;
iic = 0;
}
else
{
irc = 0;
iic = 1;
}
csfftc.mm = M[k];
csfftc.ks = Nf[k];
ksm = csfftc.ks - 1;
csfftc.ilast = Nf[k + 1];
for (l = 1; l <= csfftc.mm; l++)
{
Kee[l + 1] = Kee[l]/2;
}
for (j = 1, _do0=DOCNT(j,ntot2,_do1 = csfftc.ilast); _do0 > 0; j += _do1, _do0--)
{
jl = j + ksm;
for (i = j; i <= jl; i += 2)
{
ir = i + irc;
ii = i + iic;
sfft( &A[ir], &A[ii], s );
}
}
}
return;
/* End of loop for computing multiple sum
* */
} /* end of function */