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

function [n,wa,ifac]=rffti1(n,wa,ifac);
persistent arg argh argld fi firstCall gt i ib ido ii ip ipm is j k1 l1 l2 ld nf nfm1 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(ib), ib=0; end;
if isempty(ido), ido=0; end;
if isempty(ii), ii=0; end;
if isempty(ip), ip=0; end;
if isempty(ipm), ipm=0; end;
if isempty(is), is=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(nfm1), nfm1=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(gt), gt=0; end;
if isempty(ntryh), ntryh=zeros(1,4); end;
%***BEGIN PROLOGUE  RFFTI1
%***PURPOSE  Initialize a real and an integer work array for RFFTF1 and
%            RFFTB1.
%***LIBRARY   SLATEC (FFTPACK)
%***CATEGORY  J1A1
%***TYPE      SINGLE PRECISION (RFFTI1-S, CFFTI1-C)
%***KEYWORDS  FFTPACK, FOURIER TRANSFORM
%***AUTHOR  Swarztrauber, P. N., (NCAR)
%***DESCRIPTION
%
%   subroutine RFFTI1 initializes the work arrays WA and IFAC which are
%   used in both RFFTF1 and RFFTB1.  The prime factorization of N and a
%   tabulation of the trigonometric functions are computed and stored in
%   IFAC and WA, respectively.
%
%   Input Argument
%
%   N       the length of the sequence to be transformed.
%
%   Output Arguments
%
%   WA      a real work array which must be dimensioned at least N.
%
%   IFAC    an integer work array which must be dimensioned at least 15.
%
%   The same work arrays can be used for both RFFTF1 and RFFTB1 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 RFFTF1 or RFFTB1.
%
%***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 functions instead of DATA
%               statements.
%   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  RFFTI1
wa_shape=size(wa);wa=reshape(wa,1,[]);
ifac_shape=size(ifac);ifac=reshape(ifac,1,[]);
if firstCall,   ntryh(1) =[4];  end;
if firstCall,  ntryh(2) =[2];  end;
if firstCall,  ntryh(3) =[3];  end;
if firstCall,  ntryh(4)=[5];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  RFFTI1
nl = fix(n);
nf = 0;
j = 0;
gt=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 )
gt=1;
break;
end;
end;
if(gt==1)
break;
end;
end;
ifac(1) = fix(n);
ifac(2) = fix(nf);
tpi = 8..*atan(1.);
argh = tpi./n;
is = 0;
nfm1 = fix(nf - 1);
l1 = 1;
if( nfm1==0 )
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;
return;
end;
for k1 = 1 : nfm1;
ip = fix(ifac(k1+2));
ld = 0;
l2 = fix(l1.*ip);
ido = fix(fix(n./l2));
ipm = fix(ip - 1);
for j = 1 : ipm;
ld = fix(ld + l1);
i = fix(is);
argld = ld.*argh;
fi = 0.;
for ii = 3 : 2: ido ;
i = fix(i + 2);
fi = fi + 1.;
arg = fi.*argld;
wa(i-1) = cos(arg);
wa(i) = sin(arg);
end; ii = fix(ido +1);
is = fix(is + ido);
end; j = fix(ipm+1);
l1 = fix(l2);
end; k1 = fix(nfm1+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 rffti1
%DECK RFFTI