[n,wa,ifac]=cffti1(n,wa,ifac);

function [n,wa,ifac]=cffti1(n,wa,ifac);
persistent arg argh argld fi firstCall i i1 ib ido idot igo ii ip ipm j k1 l1 l2 ld nf nl nq nr ntry ntryh tpi ; if isempty(firstCall),firstCall=1;end; 

if isempty(arg), arg=0; end;
if isempty(argh), argh=0; end;
if isempty(argld), argld=0; end;
if isempty(fi), fi=0; end;
if isempty(tpi), tpi=0; end;
if isempty(i), i=0; end;
if isempty(i1), i1=0; end;
if isempty(ib), ib=0; end;
if isempty(ido), ido=0; end;
if isempty(idot), idot=0; end;
if isempty(ii), ii=0; end;
if isempty(ip), ip=0; end;
if isempty(ipm), ipm=0; end;
if isempty(j), j=0; end;
if isempty(k1), k1=0; end;
if isempty(l1), l1=0; end;
if isempty(l2), l2=0; end;
if isempty(ld), ld=0; end;
if isempty(nf), nf=0; end;
if isempty(nl), nl=0; end;
if isempty(nq), nq=0; end;
if isempty(nr), nr=0; end;
if isempty(ntry), ntry=0; end;
if isempty(igo), igo=0; end;
if isempty(ntryh), ntryh=zeros(1,4); end;
%***BEGIN PROLOGUE  CFFTI1
%***PURPOSE  Initialize a real and an integer work array for CFFTF1 and
%            CFFTB1.
%***LIBRARY   SLATEC (FFTPACK)
%***CATEGORY  J1A2
%***TYPE      COMPLEX (RFFTI1-S, CFFTI1-C)
%***KEYWORDS  FFTPACK, FOURIER TRANSFORM
%***AUTHOR  Swarztrauber, P. N., (NCAR)
%***DESCRIPTION
%
%  subroutine CFFTI1 initializes the work arrays WA and IFAC which are
%  used in both CFFTF1 and CFFTB1.  The prime factorization of N and a
%  tabulation of the trigonometric functions are computed and stored in
%  IFAC and WA, respectively.
%
%  Input Parameter
%
%  N       the length of the sequence to be transformed
%
%  Output Parameters
%
%  WA      a real work array which must be dimensioned at least 2*N.
%
%  IFAC    an integer work array which must be dimensioned at least 15.
%
%          The same work arrays can be used for both CFFTF1 and CFFTB1
%          as long as N remains unchanged.  Different WA and IFAC arrays
%          are required for different values of N.  The contents of
%          WA and IFAC must not be changed between calls of CFFTF1 or
%          CFFTB1.
%
%***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
%                 Computations (G. Rodrigue, ed.), Academic Press,
%                 1982, pp. 51-83.
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   790601  DATE WRITTEN
%   830401  Modified to use SLATEC library source file format.
%   860115  Modified by Ron Boisvert to adhere to Fortran 77 by
%           (a) changing dummy array size declarations (1) to (*),
%           (b) changing references to intrinsic function FLOAT
%               to REAL, and
%           (c) changing definition of variable TPI by using
%               FORTRAN intrinsic function ATAN instead of a DATA
%               statement.
%   881128  Modified by Dick Valent to meet prologue standards.
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900131  Routine changed from subsidiary to user-callable.  (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  CFFTI1
wa_shape=size(wa);wa=reshape(wa,1,[]);
ifac_shape=size(ifac);ifac=reshape(ifac,1,[]);
if firstCall,   ntryh(1) =[3];  end;
if firstCall,  ntryh(2) =[4];  end;
if firstCall,  ntryh(3) =[2];  end;
if firstCall,  ntryh(4)=[5];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  CFFTI1
nl = fix(n);
nf = 0;
j = 0;
igo=0;
while( true );
j = fix(j + 1);
if( j<=4 )
ntry = fix(ntryh(j));
else;
ntry = fix(ntry + 2);
end;
while( true );
nq = fix(fix(nl./ntry));
nr = fix(nl - ntry.*nq);
if( nr~=0 )
break;
end;
nf = fix(nf + 1);
ifac(nf+2) = fix(ntry);
nl = fix(nq);
if( ntry==2 )
if( nf~=1 )
for i = 2 : nf;
ib = fix(nf - i + 2);
ifac(ib+2) = fix(ifac(ib+1));
end; i = fix(nf+1);
ifac(3) = 2;
end;
end;
if( nl==1 )
igo=1;
break;
end;
end;
if(igo==1)
break;
end;
end;
ifac(1) = fix(n);
ifac(2) = fix(nf);
tpi = 8..*atan(1.);
argh = tpi./n;
i = 2;
l1 = 1;
for k1 = 1 : nf;
ip = fix(ifac(k1+2));
ld = 0;
l2 = fix(l1.*ip);
ido = fix(fix(n./l2));
idot = fix(ido + ido + 2);
ipm = fix(ip - 1);
for j = 1 : ipm;
i1 = fix(i);
wa(i-1) = 1.;
wa(i) = 0.;
ld = fix(ld + l1);
fi = 0.;
argld = ld.*argh;
for ii = 4 : 2: idot ;
i = fix(i + 2);
fi = fi + 1.;
arg = fi.*argld;
wa(i-1) = cos(arg);
wa(i) = sin(arg);
end; ii = fix(idot +1);
if( ip>5 )
wa(i1-1) = wa(i-1);
wa(i1) = wa(i);
end;
end; j = fix(ipm+1);
l1 = fix(l2);
end; k1 = fix(nf+1);
wa_shape=zeros(wa_shape);wa_shape(:)=wa(1:numel(wa_shape));wa=wa_shape;
ifac_shape=zeros(ifac_shape);ifac_shape(:)=ifac(1:numel(ifac_shape));ifac=ifac_shape;
end %subroutine cffti1
%DECK CFFTI