- CFFTMI multiple complex initialization
- CFFTMB multiple complex backward
- CFFTMF multiple complex forward

- COSTMI multiple real cosine initialization
- COSTMB multiple real cosine backward
- COSTMF multiple real cosine forward

- SINTMI multiple real sine initialization
- SINTMB multiple real sine backward
- SINTMF multiple real sine forward

- COSQ1I 1D real quarter-cosine initialization
- COSQ1B 1D real quarter-cosine backward
- COSQ1F 1D real quarter-cosine forward

- COSQMI multiple real quarter-cosine initialization
- COSQMB multiple real quarter-cosine backward
- COSQMF multiple real quarter-cosine forward

C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SUBROUTINE CFFT1I (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine CFFT1I initializes array WSAVE for use in its companion routines CFFT1B and CFFT1F. Routine CFFT1I must be called before the first call to CFFT1B or CFFT1F, and after whenever the value of integer N changes. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines CFFT1B or CFFT1F. IER = 0 successful exit = 2 input parameter LENSAV not big enough

C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SUBROUTINE CFFT1B (N, INC, C, LENC, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENC, LENSAV, LENWRK, IER COMPLEX C(LENC) REAL WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine CFFT1B computes the one-dimensional Fourier transform of a single periodic sequence within a complex array. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequence from spectral to physical space. This transform is normalized since a call to CFFT1B followed by a call to CFFT1F (or vice-versa) reproduces the original array subject to algorithm constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array C, of two consecutive elements within the sequence to be transformed. C Complex array of length LENC containing the sequence to be transformed. LENC Integer dimension of C array. LENC must be at least INC*(N-1) + 1. WSAVE Real work array with dimension LENSAV. WSAVE's contents must be initialized with a call to subroutine CFFT1I before the first call to routine CFFT1F or CFFT1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to CFFT1F and CFFT1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least 2*N. Output Arguments C For index J*INC+1 where J=0,...,N-1, C(J*INC+1) = N-1 SUM C(K*INC+1)*EXP(I*J*K*2*PI/N) K=0 where I=SQRT(-1). At other indices, the output value of C does not differ from input. IER = 0 successful exit = 1 input parameter LENC not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SUBROUTINE CFFT1F (N, INC, C, LENC, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENC, LENSAV, LENWRK, IER COMPLEX C(LENC) REAL WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine CFFT1F computes the one-dimensional Fourier transform of a single periodic sequence within a complex array. This transform is referred to as the forward transform or Fourier analysis, transforming the sequence from physical to spectral space. This transform is normalized since a call to CFFT1F followed by a call to CFFT1B (or vice-versa) reproduces the original array subject to algorithm constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array C, of two consecutive elements within the sequence to be transformed. C Complex array of length LENC containing the sequence to be transformed. LENC Integer dimension of C array. LENC must be at least INC*(N-1) + 1. WSAVE Real work array with dimension LENSAV. WSAVE's contents must be initialized with a call to subroutine CFFT1I before the first call to routine CFFT1F or CFFT1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to CFFT1F and CFFT1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least 2*N. Output Arguments C For index J*INC+1 where J=0,...,N-1 (that is, for the Jth element of the sequence), C(J*INC+1) = N-1 SUM C(K*INC+1)*EXP(-I*J*K*2*PI/N) K=0 where I=SQRT(-1). At other indices, the output value of C does not differ from input. IER = 0 successful exit = 1 input parameter LENC not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

SUBROUTINE CFFT2I (L, M, WSAVE, LENSAV, IER) INTEGER L, M, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 routine CFFT2I initializes real array WSAVE for use in its companion routines CFFT2F and CFFT2B for computing two- dimensional fast Fourier transforms of complex data. Prime factorizations of L and M, together with tabulations of the trigonometric functions, are computed and stored in array WSAVE. Input Arguments L Integer number of elements to be transformed in the first dimension. The transform is most efficient when L is a product of small primes. M Integer number of elements to be transformed in the second dimension. The transform is most efficient when M is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*(L+M) + INT(LOG(REAL(L))/LOG(2.)) + INT(LOG(REAL(M))/LOG(2.)) + 8. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of L and M, and also containing certain trigonometric values which will be used in routines CFFT2B or CFFT2F. WSAVE Real work array with dimension LENSAV. The WSAVE array must be initialized with a call to subroutine CFFT2I before the first call to CFFT2B or CFFT2F, and thereafter whenever the values of L, M or the contents of array WSAVE change. Using different WSAVE arrays for different transform lengths or types in the same program may reduce computation costs because the array contents can be re-used. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

SUBROUTINE CFFT2B (LDIM, L, M, C, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER L, M, LDIM, LENSAV, LENWRK, IER COMPLEX C(LDIM,M) REAL WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine CFFT2B computes the two-dimensional discrete Fourier transform of a complex periodic array. This transform is known as the backward transform or Fourier synthesis, transforming from spectral to physical space. Routine CFFT2B is normalized, in that a call to CFFT2B followed by a call to CFFT2F (or vice-versa) reproduces the original array subject to algorithm constraints, roundoff error, etc. Input Arguments LDIM Integer first dimension of two-dimensional complex array C. L Integer number of elements to be transformed in the first dimension of the two-dimensional complex array C. The value of L must be less than or equal to that of LDIM. The transform is most efficient when L is a product of small primes. M Integer number of elements to be transformed in the second dimension of the two-dimensional complex array C. The transform is most efficient when M is a product of small primes. C Complex array of two dimensions containing the (L,M) subarray to be transformed. C's first dimension is LDIM, its second dimension must be at least M. WSAVE Real work array with dimension LENSAV. WSAVE's contents must be initialized with a call to subroutine CFFT2I before the first call to routine CFFT2F or CFFT2B with transform lengths L and M. WSAVE's contents may be re-used for subsequent calls to CFFT2F and CFFT2B with the same transform lengths L and M. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*(L+M) + INT(LOG(REAL(L))/LOG(2.)) + INT(LOG(REAL(M))/LOG(2.)) + 8. WORK Real work array. LENWRK Integer dimension of WORK array. LENWRK must be at least 2*L*M. Output Arguments C Complex output array. For purposes of exposition, assume the index ranges of array C are defined by C(0:L-1,0:M-1). For I=0,...,L-1 and J=0,...,M-1, the C(I,J)'s are given in the traditional aliased form by L-1 M-1 C(I,J) = SUM SUM C(L1,M1)* L1=0 M1=0 EXP(SQRT(-1)*2*PI*(I*L1/L + J*M1/M)) And in unaliased form, the C(I,J)'s are given by LF MF C(I,J) = SUM SUM C(L1,M1,K1)* L1=LS M1=MS EXP(SQRT(-1)*2*PI*(I*L1/L + J*M1/M)) where LS= -L/2 and LF=L/2-1 if L is even; LS=-(L-1)/2 and LF=(L-1)/2 if L is odd; MS= -M/2 and MF=M/2-1 if M is even; MS=-(M-1)/2 and MF=(M-1)/2 if M is odd; and C(L1,M1) = C(L1+L,M1) if L1 is zero or negative; C(L1,M1) = C(L1,M1+M) if M1 is zero or negative; The two forms give different results when used to interpolate between elements of the sequence. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 5 input parameter L > LDIM = 20 input error returned by lower level routine

SUBROUTINE CFFT2F (LDIM, L, M, C, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER L, M, LDIM, LENSAV, LENWRK, IER COMPLEX C(LDIM,M) REAL WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine CFFT2F computes the two-dimensional discrete Fourier transform of a complex periodic array. This transform is known as the forward transform or Fourier analysis, transforming from physical to spectral space. Routine CFFT2F is normalized, in that a call to CFFT2F followed by a call to CFFT2B (or vice-versa) reproduces the original array subject to algorithm constraints, roundoff error, etc. Input Arguments LDIM Integer first dimension of two-dimensional complex array C. L Integer number of elements to be transformed in the first dimension of the two-dimensional complex array C. The value of L must be less than or equal to that of LDIM. The transform is most efficient when L is a product of small primes. M Integer number of elements to be transformed in the second dimension of the two-dimensional complex array C. The transform is most efficient when M is a product of small primes. C Complex array of two dimensions containing the (L,M) subarray to be transformed. C's first dimension is LDIM, its second dimension must be at least M. WSAVE Real work array with dimension LENSAV. WSAVE's contents must be initialized with a call to subroutine CFFT2I before the first call to routine CFFT2F or CFFT2B with transform lengths L and M. WSAVE's contents may be re-used for subsequent calls to CFFT2F and CFFT2B having those same transform lengths. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*(L+M) + INT(LOG(REAL(L))/LOG(2.)) + INT(LOG(REAL(M))/LOG(2.)) + 8. WORK Real work array. LENWRK Integer dimension of WORK array. LENWRK must be at least 2*L*M. Output Arguments C Complex output array. For purposes of exposition, assume the index ranges of array C are defined by C(0:L-1,0:M-1). For I=0,...,L-1 and J=0,...,M-1, the C(I,J)'s are given in the traditional aliased form by L-1 M-1 C(I,J) = 1/(L*M)*SUM SUM C(L1,M1)* L1=0 M1=0 EXP(-SQRT(-1)*2*PI*(I*L1/L + J*M1/M)) And in unaliased form, the C(I,J)'s are given by LF MF C(I,J) = 1/(L*M)*SUM SUM C(L1,M1)* L1=LS M1=MS EXP(-SQRT(-1)*2*PI*(I*L1/L + J*M1/M)) where LS= -L/2 and LF=L/2-1 if L is even; LS=-(L-1)/2 and LF=(L-1)/2 if L is odd; MS= -M/2 and MF=M/2-1 if M is even; MS=-(M-1)/2 and MF=(M-1)/2 if M is odd; and C(L1,M1) = C(L1+L,M1) if L1 is zero or negative; C(L1,M1) = C(L1,M1+M) if M1 is zero or negative; The two forms give different results when used to interpolate between elements of the sequence. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 5 input parameter L > LDIM = 20 input error returned by lower level routine

SUBROUTINE CFFTMI (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine CFFTMI initializes array WSAVE for use in its companion routines CFFTMB and CFFTMF. Routine CFFTMI must be called before the first call to CFFTMB or CFFTMF, and after whenever the value of integer N changes. Input Arguments N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines CFFTMB or CFFTMF. IER = 0 successful exit = 2 input parameter LENSAV not big enough

SUBROUTINE CFFTMB (LOT, JUMP, N, INC, C, LENC, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENC, LENSAV, LENWRK, IER COMPLEX C(LENC) REAL WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine CFFTMB computes the one-dimensional Fourier transform of multiple periodic sequences within a complex array. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to CFFTMF followed by a call to CFFTMB (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array C. JUMP Integer increment between the locations, in array C, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array C, of two consecutive elements within the same sequence to be transformed. C Complex array containing LOT sequences, each having length N, to be transformed. C can have any number of dimensions, but the total number of locations must be at least LENC. LENC Integer dimension of C array. LENC must be at least (LOT-1)*JUMP + INC*(N-1) + 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine CFFTMI before the first call to routine CFFTMF or CFFTMB for a given transform length N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least 2*LOT*N. Output Arguments C For index L*JUMP+J*INC+1 where J=0,...,N-1 and L=0,...,LOT-1, (that is, for the Jth element of the Lth sequence), C(L*JUMP+J*INC+1) = N-1 SUM C(L*JUMP+K*INC+1)*EXP(I*J*K*2*PI/N) K=0 where I=SQRT(-1). At other indices, the output value of C does not differ from input. IER = 0 successful exit = 1 input parameter LENC not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once.

SUBROUTINE CFFTMF (LOT, JUMP, N, INC, C, LENC, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENC, LENSAV, LENWRK, IER COMPLEX C(LENC) REAL WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine CFFTMF computes the one-dimensional Fourier transform of multiple periodic sequences within a complex array. This transform is referred to as the forward transform or Fourier analysis, transforming the sequences from physical to spectral space. This transform is normalized since a call to CFFTMF followed by a call to CFFTMB (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array C. JUMP Integer increment between the locations, in array C, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array C, of two consecutive elements within the same sequence to be transformed. C Complex array containing LOT sequences, each having length N, to be transformed. C can have any number of dimensions, but the total number of locations must be at least LENC. LENC Integer dimension of C array. LENC must be at least (LOT-1)*JUMP + INC*(N-1) + 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine CFFTMI before the first call to routine CFFTMF or CFFTMB for a given transform length N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least 2*LOT*N. Output Arguments C For index L*JUMP + J*INC +1 where J=0,...,N-1 and L=0,...,LOT-1, (that is, for the Jth element of the Lth sequence), C(L*JUMP+J*INC+1) = N-1 SUM C(L*JUMP+K*INC+1)*EXP(-I*J*K*2*PI/N) K=0 where I=SQRT(-1). At other indices, the output value of C does not differ from input. IER = 0 successful exit = 1 input parameter LENC not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once.

SUBROUTINE RFFT1I (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine RFFT1I initializes array WSAVE for use in its companion routines RFFT1B and RFFT1F. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines RFFT1B or RFFT1F. IER = 0 successful exit = 2 input parameter LENSAV not big enough

SUBROUTINE RFFT1B (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine RFFT1B computes the one-dimensional Fourier transform of a periodic sequence within a real array. This is referred to as the backward transform or Fourier synthesis, transforming the sequence from spectral to physical space. This transform is normalized since a call to RFFT1B followed by a call to RFFT1F (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1) + 1. WSAVE Real work array o length LENSAV. WSAVE's contents must be initialized with a call to subroutine RFFT1I before the first call to routine RFFT1F or RFFT1B for a given transform length N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(N-1)*INC). The output values of R are written over the input values. If N is even, set NH=N/2-1; then for J=0,...,N-1 R(J*INC) = R(0) + [(-1)**J*R((N-1)*INC)] NH + SUM R((2*N1-1)*INC)*COS(J*N1*2*PI/N) N1=1 NH + SUM R(2*N1*INC)*SIN(J*N1*2*PI/N) N1=1 If N is odd, set NH=(N-1)/2 and define R as above, except remove the expression in square brackets []. IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough

SUBROUTINE RFFT1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine RFFT1F computes the one-dimensional Fourier transform of a periodic sequence within a real array. This is referred to as the forward transform or Fourier analysis, transforming the sequence from physical to spectral space. This transform is normalized since a call to RFFT1F followed by a call to RFFT1B (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1) + 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine RFFT1I before the first call to routine RFFT1F or RFFT1B for a given transform length N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(N-1)*INC). Then N-1 R(0) = SUM R(N1*INC)/N N1=0 If N is even, set NH=N/2-1; if N is odd set NH=(N-1)/2; then for J=1,...,NH R((2*J-1)*INC) = N-1 2.*SUM (R(N1*INC)*COS(J*N1*2*PI/N)/N N1=0 and R(2*J*INC) = N-1 2.*SUM (R(N1*INC)*SIN(J*N1*2*PI/N)/N N1=0 Also if N is even then R((N-1)*INC) = N-1 SUM (-1)**N1*R(N1*INC)/N N1=0 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough

SUBROUTINE RFFT2I (L, M, WSAVE, LENSAV, IER) INTEGER L, M, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 routine RFFT2I initializes real array WSAVE for use in its companion routines RFFT2F and RFFT2B for computing the two- dimensional fast Fourier transform of real data. Prime factorizations of L and M, together with tabulations of the trigonometric functions, are computed and stored in array WSAVE. RFFT2I must be called prior to the first call to RFFT2F or RFFT2B. Separate WSAVE arrays are required for different values of L or M. Input Arguments L Integer number of elements to be transformed in the first dimension. The transform is most efficient when L is a product of small primes. M Integer number of elements to be transformed in the second dimension. The transform is most efficient when M is a product of small primes. LENSAV Integer number of elements in the WSAVE array. LENSAV must be at least L + 3*M + INT(LOG(REAL(L))/LOG(2.)) + 2*INT(LOG(REAL(M))/LOG(2.)) +12. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of L and M, and also containing certain trigonometric values which will be used in routines RFFT2B or RFFT2F. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

backward fast Fourier transform

SUBROUTINE RFFT2B (LDIM, L, M, R, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LDIM, L, M, LENSAV, LENWRK, IER REAL R(LDIM,M), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine RFFT2B computes the two-dimensional discrete Fourier transform of the complex Fourier coefficients a real periodic array. This transform is known as the backward transform or Fourier synthesis, transforming from spectral to physical space. Routine RFFT2B is normalized: a call to RFFT2B followed by a call to RFFT2F (or vice-versa) reproduces the original array subject to algorithmic sonstraints, roundoff error, etc. Input Arguments LDIM Integer first dimension of the two-dimensional real array R, which must be at least L. L Integer number of elements to be transformed in the first dimension of the two-dimensional real array R. The value of L must be less than or equal to that of LDIM. The transform is most efficient when L is a product of small primes. M Integer number of elements to be transformed in the second dimension of the two-dimensional real array R. The transform is most efficient and accurate when M is a product of small primes. R A real L by M array containing the spectral coefficients of a real L by M array that are stored as described in the documentation of subroutine rfft2f. THe first dimension is LDIM which must be at least L. The second dimension must be at least M. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine RFFT2I before the first call to routine RFFT2F or RFFT2B. WSAVE's contents may be re-used for subsequent calls to RFFT2F and RFFT2B with the same L and M. LENSAV Integer number of elements in the WSAVE array. LENSAV must be at least L + 3*M + INT(LOG(REAL(L))/LOG(2.)) + 2*INT(LOG(REAL(M))/LOG(2.)) +12. WORK Real array of dimension LENWRK, where LENWRK is defined below. WORK provides workspace, and its contents need not be saved between calls to routines RFFT2B and RFFT2F. LENWRK Integer number of elements in the WORK array. LENWRK must be at least (L+1)*M. Output Arguments R A real L by M array. If the full transform c is reconstructed using subroutine r2c(ldim,lcdim,l,m,r,c) then for i=0,...,l-1 and j=0,...,m-1 L-1 M-1 R(I,J) = SUM SUM C(L1,M1) L1=0 M1=0 *EXP(SQRT(-1)*2*PI*(I*L1/L+J*M1/M)) or using the conjugate symmetry c(i,j)=c(l-1,m-j) this can be written in terms of c(i,j), i=0,...,(L+1)/2 and j=0,...,m as: M-1 R(I,J) = REAL[ SUM C(0,M1)*EXP(SQRT(-1)*2*PI*J*M1/M ] M1=0 (L+1)/2-1 M-1 + 2*REAL[ SUM SUM C(L1,M1)* L1=1 M1=0 *EXP(SQRT(-1)*2*PI*(I*L1/L+J*M1/M)) ] If L is even then add M-1 + REAL[ SUM (-1)**I*C(L/2,M1)*EXP(SQRT(-1)*2*PI*J*M1/M) ] M1=0 c(i,j) = a(i,j)+i*b(I,J) are contained in the real output array r(i,j) except for c(0,m-j) and c(l,m-j) which can be obtained from c(0,m-j) = c(0,j) and c(l,j) = c(l.m-j) = c(l,j) IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 6 input parameter LDIM is less than 2*INT((L+1)/2) = 20 input error returned by lower level routine

forward fast Fourier transform

SUBROUTINE RFFT2F (LDIM, L, M, R, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LDIM, L, M, LENSAV, LENWRK, IER REAL R(LDIM,M), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine RFFT2F computes the two-dimensional discrete Fourier transform of a real periodic array. This transform is known as the forward transform or Fourier analysis, transforming from physical to spectral space. Routine RFFT2F is normalized: a call to RFFT2F followed by a call to RFFT2B (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LDIM Integer first dimension of the two-dimensional real array R, which must be at least L. L Integer number of elements to be transformed in the first dimension of the two-dimensional real array R. The value of L must be less than or equal to that of LDIM. The transform is most efficient when L is a product of small primes. M Integer number of elements to be transformed in the second dimension of the two-dimensional real array R. The transform is most efficient when M is a product of small primes. R Real array of two dimensions containing the L-by-M subarray to be transformed. R's first dimension is LDIM and its second dimension must be at least as large as M. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine RFFT2I before the first call to routine RFFT2F or RFFT2B. WSAVE's contents may be re-used for subsequent calls to RFFT2F and RFFT2B as long as L and M remain unchanged. LENSAV Integer number of elements in the WSAVE array. LENSAV must be at least L + 3*M + INT(LOG(REAL(L))/LOG(2.)) + 2*INT(LOG(REAL(M))/LOG(2.)) +12. WORK Real array of dimension LENWRK which is defined below. WORK provides workspace, and its contents need not be saved between calls to routines RFFT2F and RFFT2B. LENWRK Integer number of elements in the WORK array. LENWRK must be at least (L+1)*M. Output Arguments R Real output array of two dimensions. The full complex transform of r(i,j) is given by: L-1 M-1 C(I,J) = 1/(L*M)*SUM SUM R(L1,M1)* L1=0 M1=0 EXP(-SQRT(-1)*2*PI*(I*L1/L + J*M1/M)) The complex transform of a real array has conjugate symmetry. That is, c(i,j) = congugate c(l-i,m-j) so only half the transform is computed and packed back into the original array R. Examples: Let a(i,j) = re[c(i,j)] and b(i,j) = Im[c(i,j)] then following the forward transform For l=m=4 a(0,0) a(0,1) b(0,1) a(0,2) r(i,j) = a(1,0) a(1,1) a(1,2) a(1,3) b(1,0) b(1,1) b(1,2) b(1,3) a(2,0) a(2,1) b(2,1) a(2,2) For l=m=5 a(0,0) a(0,1) b(0,1) a(0,2),b(0,2) a(1,0) a(1,1) a(1,2) a(1,3),a(1,4) r(i,j) = b(1,0) b(1,1) b(1,2) b(1,3),b(1,4) a(2,0) a(2,1) a(2,2) a(2,3),a(2,4) b(2,0) b(2,1) b(2,2) b(2,3),b(2,4) The remaining c(i,j) for i=int((L+1)/2)+1,..,L and m=0,...m-1 can be obtained via the conjugate symmetry, which also implies that c(0,j) = conjugate c(0,m-j) and for even l, c(l/2,0) = conjugate c(l/2,m-j). The full complex transform c(i,j), i=1,...,L and j=1,...,M can also be extracted using subroutine r2c(ldim,lcdim,l,m,r,c) where lcdim is the first dimension of the complex array c, which must be greater than or equal to l. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 6 input parameter LDIM is less than 2*INT((L+1)/2) = 20 input error returned by lower level routine

SUBROUTINE RFFTMI (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine RFFTMI initializes array WSAVE for use in its companion routines RFFTMB and RFFTMF. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines RFFTMB or RFFTMF. IER = 0 successful exit = 2 input parameter LENSAV not big enough

SUBROUTINE RFFTMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine RFFTMB computes the one-dimensional Fourier transform of multiple periodic sequences within a real array. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to RFFTMB followed by a call to RFFTMF (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1) + 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine RFFTMI before the first call to routine RFFTMF or RFFTMB for a given transform length N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(LOT-1)*JUMP+(N-1)*INC). The output values of R are written over the input values. If N is even, set NH=N/2-1; then for I=0,...,LOT-1 and J=0,...,N-1 R(I*JUMP+J*INC) = R(I*JUMP) + [(-1)**J*R(I*JUMP+(N-1)*INC)] NH + SUM R(I*JUMP+(2*N1-1)*INC)*COS(J*N1*2*PI/N) N1=1 NH + SUM R(I*JUMP+2*N1*INC)*SIN(J*N1*2*PI/N) N1=1 If N is odd, set NH=(N-1)/2 and define R as above, except remove the expression in square brackets []. IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once.

SUBROUTINE RFFTMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine RFFTMF computes the one-dimensional Fourier transform of multiple periodic sequences within a real array. This transform is referred to as the forward transform or Fourier analysis, transforming the sequences from physical to spectral space. This transform is normalized since a call to RFFTMF followed by a call to RFFTMB (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1) + 1. WSAVE Real work array o length LENSAV. WSAVE's contents must be initialized with a call to subroutine RFFTMI before the first call to routine RFFTMF or RFFTMB for a given transform length N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(LOT-1)*JUMP+(N-1)*INC). Then for I=0,...,LOT-1 N-1 R(I*JUMP) = SUM R(I*JUMP+N1*INC)/N N1=0 If N is even, set NH=N/2-1; if N is odd set NH=(N-1)/2; then for J=1,...,NH R(I*JUMP+(2*J-1)*INC) = N-1 2.*SUM (R(I*JUMP+N1*INC)*COS(J*N1*2*PI/N)/N N1=0 and R(I*JUMP+2*J*INC) = N-1 2.*SUM (R(I*JUMP+N1*INC)*SIN(J*N1*2*PI/N)/N N1=0 Also if N is even then R(I*JUMP+(N-1)*INC) = N-1 SUM (-1)**N1*R(I*JUMP+N1*INC)/N N1=0 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once.

SUBROUTINE COST1I (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine COST1I initializes array WSAVE for use in its companion routines COST1F and COST1B. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N-1 is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines COST1B or COST1F. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

SUBROUTINE COST1B (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COST1B computes the one-dimensional Fourier transform of an even sequence within a real array. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequence from spectral to physical space. This transform is normalized since a call to COST1B followed by a call to COST1F (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N-1 is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine COST1I before the first call to routine COST1F or COST1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COST1F and COST1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least N-1. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(N-1)*INC). The output values of R are written over the input values. For J=0,...,N-1 R(J*INC) = N-1 SUM R(N1*INC)*COS(J*N1*PI/(N-1)) N1=0 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

SUBROUTINE COST1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COST1F computes the one-dimensional Fourier transform of an even sequence within a real array. This transform is referred to as the forward transform or Fourier analysis, transforming the sequence from physical to spectral space. This transform is normalized since a call to COST1F followed by a call to COST1B (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N-1 is a product of small primes. INC Integer increment between the locations, in array R of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine COST1I before the first call to routine COST1F or COST1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COST1F and COST1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least N-1. Output Arguments R Real output array R. For purposes of this exposition, assume R's range of indices is given by R(0:(N-1)*INC) and that the input values of R are denoted by X. The output values of R are written over the input values X as follows: CASE N=1 -------- R(0) = X(0) CASE N>1 -------- For J=0 R(0) = 0.5*X(0)/(N-1) N-2 + SUM R(N1*INC)/(N-1) N1=1 + 0.5*X((N-1)*INC)/(N-1) For J=1,...,N-2 R(J*INC) = R(0)/(N-1) N-2 + SUM 2.0*(X(N1*INC)*COS(J*N1*PI/(N-1)))/(N-1) N1=1 + ((-1)**J)*X((N-1)*INC)/(N-1) R((N-1)*INC) = 0.5*X(0)/(N-1) N-2 + SUM R(N1*INC)*((-1)**N1)/(N-1) N1=1 + 0.5*((-1)**(N-1))*X((N-1)*INC)/(N-1) IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

SUBROUTINE COSTMI (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine COSTMI initializes array WSAVE for use in its companion routines COSTMF and COSTMB. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4 Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines COSTMB or COSTMF. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

SUBROUTINE COSTMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COSTMB computes the one-dimensional Fourier transform of multiple even sequences within a real array. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to COSTMB followed by a call to COSTMF (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N-1 is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine COSTMI before the first call to routine COSTMF or COSTMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COSTMF and COSTMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*(N+1). Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(LOT-1)*JUMP+(N-1)*INC). The output values of R are written over the input values. For I=0,...,LOT-1 and J=0,...,N-1 R(I*JUMP+J*INC) = N-1 SUM R(I*JUMP+N1*INC)*COS(J*N1*PI/(N-1)) N1=0 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. = 20 input error returned by lower level routine The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent, otherwise at least one array element mistakenly is transformed more than once.

SUBROUTINE COSTMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COSTMF computes the one-dimensional Fourier transform of multiple even sequences within a real array. This transform is referred to as the forward transform or Fourier analysis, transforming the sequences from physical to spectral space. This transform is normalized since a call to COSTMF followed by a call to COSTMB (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N-1 is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine COSTMI before the first call to routine COSTMF or COSTMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COSTMF and COSTMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*(N+1). Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(LOT-1)*JUMP+(N-1)*INC). The output values of R are written over the input values. For I=0,...,LOT-1 R(I*JUMP) = 0.5*X(I*JUMP)/(N-1) N-2 + SUM R(I*JUMP+*N1*INC)/(N-1) N1=1 + 0.5*X(I*JUMP+(N-1)*INC)/(N-1) For I=0,...,LOT-1 and J=1,...,N-2 R(I*JUMP+J*INC) = R(I*JUMP)/(N-1) N-2 + SUM 2.0*(X(I*JUMP+*N1*INC)*COS(J*N1*PI/(N-1)))/(N-1) N1=1 + ((-1)**J)*X(I*JUMP+(N-1)*INC)/(N-1) For I=0,...,LOT-1 R(I*JUMP+(N-1)*INC) = 0.5*X(I*JUMP)/(N-1) N-2 + SUM R(I*JUMP+*N1*INC)*((-1)**N1)/(N-1) N1=1 + 0.5*((-1)**(N-1))*X(I*JUMP+(N-1)*INC)/(N-1) IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. = 20 input error returned by lower level routine The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent, otherwise at least one array element mistakenly is transformed more than once.

SUBROUTINE SINT1I (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine SINT1I initializes array WSAVE for use in its companion routines SINT1F and SINT1B. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N+1 is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines SINT1B or SINT1F. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

SUBROUTINE SINT1B (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINT1B computes the one-dimensional Fourier transform of an odd sequence within a real array. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequence from spectral to physical space. This transform is normalized since a call to SINT1B followed by a call to SINT1F (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N+1 is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINT1I before the first call to routine SINT1F or SINT1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINT1F and SINT1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. Must be at least 2*N+2. Output Arguments R Real output array. For purposes of exposition, assume R's range of indices is given by R(INC:N*INC). The output values of R are written over the input values. For J=1,...,N R(J*INC) = N SUM R(N1*INC)*SIN(J*N1*PI/(N+1)) N1=1 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

SUBROUTINE SINT1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINT1F computes the one-dimensional Fourier transform of an odd sequence within a real array. This transform is referred to as the forward transform or Fourier analysis, transforming the sequence from physical to spectral space. This transform is normalized since a call to SINT1F followed by a call to SINT1B (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N+1 is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINT1I before the first call to routine SINT1F or SINT1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINT1F and SINT1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. Must be at least 2*N+2. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(INC:(N-1)*INC). The output values of R are written over the input values. For J=1,...,N R(J*INC) = N SUM 2.*R(N1*INC)*SIN(J*N1*PI/(N+1))/(N+1) N1=1 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

SUBROUTINE SINTMI (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine SINTMI initializes array WSAVE for use in its companion routines SINTMF and SINTMB. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines SINTMB or SINTMF. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

SUBROUTINE SINTMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINTMB computes the one-dimensional Fourier transform of multiple odd sequences within a real array. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to SINTMB followed by a call to SINTMF (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences. N Integer length of each sequence to be transformed. The transform is most efficient when N+1 is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINTMI before the first call to routine SINTMF or SINTMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINTMF and SINTMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*(2*N+4). Output Arguments R Real output array. For purposes of exposition, assume R's range of indices is given by R(INC:(LOT-1)*JUMP+N*INC). The output values of R are written over the input values. For I=0,...,LOT-1 and J=1,...,N R(I*JUMP+J*INC) = N SUM R(I*JUMP+*N1*INC)*SIN(J*N1*PI/(N+1)) N1=1 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. = 20 input error returned by lower level routine The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once.

SUBROUTINE SINTMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINTMF computes the one-dimensional Fourier transform of multiple odd sequences within a real array. This transform is referred to as the forward transform or Fourier analysis, transforming the sequences from physical to spectral space. This transform is normalized since a call to SINTMF followed by a call to SINTMB (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N+1 is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINTMI before the first call to routine SINTMF or SINTMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINTMF and SINTMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*(2*N+4). Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(LOT-1)*JUMP+(N-1)*INC). The output values of R are written over the input values. For I=0,...,LOT-1 and J=1,...,N R(I*JUMP+J*INC) = N SUM 2.*R(I*JUMP+*N1*INC)*SIN(J*N1*PI/(N+1))/(N+1) N1=1 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once.

SUBROUTINE COSQ1I (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine COSQ1I initializes array WSAVE for use in its companion routines COSQ1F and COSQ1B. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines COSQ1B or COSQ1F. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

SUBROUTINE COSQ1B (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COSQ1B computes the one-dimensional Fourier transform of a sequence which is a cosine series with odd wave numbers. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequence from spectral to physical space. This transform is normalized since a call to COSQ1B followed by a call to COSQ1F (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer number of elements to be transformed in the sequence. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine COSQ1I before the first call to routine COSQ1F or COSQ1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COSQ1F and COSQ1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least N. Output Arguments R Real output array. For purposes of exposition, assume R's range of indices is given by R(0:(N-1)*INC). The output values of R are written over the input values. For J=0,...,N-1 R(J*INC) = N-1 SUM R(N1*INC)*COS(J*(2*N1+1)*PI/(2*N)) N1=0 WSAVE Contains values initialized by subroutine COSQ1I that must not be destroyed between calls to routine COSQ1F or COSQ1B. IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

SUBROUTINE COSQ1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COSQ1F computes the one-dimensional Fourier transform of a sequence which is a cosine series with odd wave numbers. This transform is referred to as the forward transform or Fourier analysis, transforming the sequence from physical to spectral space. This transform is normalized since a call to COSQ1F followed by a call to COSQ1B (or vice-versa) reproduces the original array subject to algorithmic constrain, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array with dimension LENSAV. WSAVE's contents must be initialized with a call to subroutine COSQ1I before the first call to routine COSQ1F or COSQ1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COSQ1F and COSQ1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(N-1)*INC). The output values of R are written over the input values. For J=0,...,N-1 R(J*INC) = R(0)/N N-1 + SUM 2.*R(N1*INC)*COS((2*J+1)*N1*PI/(2*N))/N N1=1 WSAVE Contains values initialized by subroutine COSQ1I that must not be destroyed between calls to routine COSQ1F or COSQ1B. IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

SUBROUTINE COSQMI (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine COSQMI initializes array WSAVE for use in its companion routines COSQMF and COSQMB. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines COSQMB or COSQMF. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

SUBROUTINE COSQMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COSQMB computes the one-dimensional Fourier transform of multiple sequences, each of which is a cosine series with odd wave numbers. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to COSQMB followed by a call to COSQMF (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array with dimension LENSAV. WSAVE's contents must be initialized with a call to subroutine COSQMI before the first call to routine COSQMF or COSQMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COSQMF and COSQMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*N. Output Arguments R Real output array. For purposes of exposition, assume R's range of indices is given by R(0:(LOT-1)*JUMP+(N-1)*INC). The output values of R are written over the input values. For I=0,...,LOT-1 and J=0,...,N-1 R(I*JUMP+J*INC) = N-1 SUM R(I*JUMP+N1*INC)*COS(J*(2*N1+1)*PI/(2*N)) N1=0 WSAVE Contains values initialized by subroutine COSQMI that must not be destroyed between calls to routine COSQMF or COSQMB. IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. = 20 input error returned by lower level routine The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent, otherwise at least one array element mistakenly is transformed more than once.

SUBROUTINE COSQMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COSQMF computes the one-dimensional Fourier transform of multiple sequences within a real array, where each of the sequences is a cosine series with odd wave numbers. This transform is referred to as the forward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to COSQMF followed by a call to COSQMB (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array o length LENSAV. WSAVE's contents must be initialized with a call to subroutine COSQMI before the first call to routine COSQMF or COSQMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COSQMF and COSQMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(LOT-1)*JUMP+(N-1)*INC). The output values of R are written over the input values. For I=0,...,LOT-1 and J=0,...,N-1 R(I*JUMP+J*INC) = R(I*JUMP)/N N-1 + SUM 2.*R(I*JUMP+*N1*INC)*COS((2*J+1)*N1*PI/(2*N))/N N1=1 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. = 20 input error returned by lower level routine The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent, otherwise at least one array element mistakenly is transformed more than once.

SUBROUTINE SINQ1I (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine SINQ1I initializes array WSAVE for use in its companion routines SINQ1F and SINQ1B. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines SINQ1B or SINQ1F. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

SUBROUTINE SINQ1B (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINQ1B computes the one-dimensional Fourier transform of a sequence which is a sine series with odd wave numbers. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequence from spectral to physical space. This transform is normalized since a call to SINQ1B followed by a call to SINQ1F (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINQ1I before the first call to routine SINQ1F or SINQ1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINQ1F and SINQ1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(INC:N*INC). The output values of R are written over the input values. For J=1,...,N R(J*INC) = N SUM R(N1*INC)*SIN(J*(2*N1-1)*PI/(2*N)) N1=1 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

SUBROUTINE SINQ1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINQ1F computes the one-dimensional Fourier transform of a sequence which is a sine series of odd wave numbers. This transform is referred to as the forward transform or Fourier analysis, transforming the sequence from physical to spectral space. This transform is normalized since a call to SINQ1F followed by a call to SINQ1B (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINQ1I before the first call to routine SINQ1F or SINQ1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINQ1F and SINQ1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(INC:N*INC). The output values of R are written over the input values. For J=1,...,N R(J*INC) = N-1 + SUM (2.*R(N1*INC)*SIN(((2*J-1)*N1*PI/(2*N)))/N N1=1 + ((-1)**(J+1))*R(N*INC)/N IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

SUBROUTINE SINQMI (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine SINQMI initializes array WSAVE for use in its companion routines SINQMF and SINQMB. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines SINQMB or SINQMF. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

SUBROUTINE SINQMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINQMB computes the one-dimensional Fourier transform of multiple sequences within a real array, where each of the sequences is a sine series with odd wave numbers. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to SINQMB followed by a call to SINQMF (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINQMI before the first call to routine SINQMF or SINQMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINQMF and SINQMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(INC:(LOT-1)*JUMP+N*INC). The output values of R are written over the input values. For I=0,...,LOT-1 and J=1,...,N R(I*JUMP+J*INC) = N SUM R(I*JUMP+N1*INC)*SIN(J*(2*N1-1)*PI/(2*N)) N1=1 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. = 20 input error returned by lower level routine The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once.

SUBROUTINE SINQMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINQMF computes the one-dimensional Fourier transform of multiple sequences within a real array, where each sequence is a sine series with odd wave numbers. This transform is referred to as the forward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to SINQMF followed by a call to SINQMB (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINQMI before the first call to routine SINQMF or SINQMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINQMF and SINQMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(INC:(LOT-1)*JUMP+N*INC). The output values of R are written over the input values. For I=0,...,LOT-1 and J=1,...,N R(I*JUMP+J*INC) = N-1 + SUM (2.*R(I*JUMP+*N1*INC)*SIN(((2*J-1)*N1*PI/(2*N)))/N N1=1 + ((-1)**(J+1))*R(I*JUMP+N*INC)/N IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. = 20 input error returned by lower level routine The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once. FFTPACK5, Copyright (C) 2004-2011, Computational Information Systems Laboratory, University Corporation for Atmospheric Research