[n,c,ch,wa,ifac]=cfftb1(n,c,ch,wa,ifac);

function [n,c,ch,wa,ifac]=cfftb1(n,c,ch,wa,ifac);
persistent i idl1 ido idot ip iw ix2 ix3 ix4 k1 l1 l2 n2 na nac nf ; 

if isempty(i), i=0; end;
if isempty(idl1), idl1=0; end;
if isempty(ido), ido=0; end;
if isempty(idot), idot=0; end;
if isempty(ip), ip=0; end;
if isempty(iw), iw=0; end;
if isempty(ix2), ix2=0; end;
if isempty(ix3), ix3=0; end;
if isempty(ix4), ix4=0; end;
if isempty(k1), k1=0; end;
if isempty(l1), l1=0; end;
if isempty(l2), l2=0; end;
if isempty(n2), n2=0; end;
if isempty(na), na=0; end;
if isempty(nac), nac=0; end;
if isempty(nf), nf=0; end;
%***BEGIN PROLOGUE  CFFTB1
%***PURPOSE  Compute the unnormalized inverse of CFFTF1.
%***LIBRARY   SLATEC (FFTPACK)
%***CATEGORY  J1A2
%***TYPE      COMPLEX (RFFTB1-S, CFFTB1-C)
%***KEYWORDS  FFTPACK, FOURIER TRANSFORM
%***AUTHOR  Swarztrauber, P. N., (NCAR)
%***DESCRIPTION
%
%  subroutine CFFTB1 computes the backward complex discrete Fourier
%  transform (the Fourier synthesis).  Equivalently, CFFTB1 computes
%  a complex periodic sequence from its Fourier coefficients.
%  The transform is defined below at output parameter C.
%
%  A call of CFFTF1 followed by a call of CFFTB1 will multiply the
%  sequence by N.
%
%  The arrays WA and IFAC which are used by subroutine CFFTB1 must be
%  initialized by calling subroutine CFFTI1 (N, WA, IFAC).
%
%  Input Parameters
%
%  N       the length of the complex sequence C.  The method is
%          more efficient when N is the product of small primes.
%
%  C       a complex array of length N which contains the sequence
%
%  CH      a real work array of length at least 2*N
%
%  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 WA and IFAC arrays must be initialized by calling
%          subroutine CFFTI1 (N, WA, IFAC), and different WA and IFAC
%          arrays must be used for each different value of N.  This
%          initialization does not have to be repeated so long as N
%          remains unchanged.  Thus subsequent transforms can be
%          obtained faster than the first.  The same WA and IFAC arrays
%          can be used by CFFTF1 and CFFTB1.
%
%  Output Parameters
%
%  C       For J=1,...,N
%
%              C(J)=the sum from K=1,...,N of
%
%                 C(K)*EXP(I*(J-1)*(K-1)*2*PI/N)
%
%                         where I=SQRT(-1)
%
%  NOTE:   WA and IFAC contain initialization calculations which must
%          not be destroyed between calls of subroutine CFFTF1 or CFFTB1
%
%***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
%                 Computations (G. Rodrigue, ed.), Academic Press,
%                 1982, pp. 51-83.
%***ROUTINES CALLED  PASSB, PASSB2, PASSB3, PASSB4, PASSB5
%***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
%           changing dummy array size declarations (1) to (*).
%   881128  Modified by Dick Valent to meet prologue standards.
%   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  CFFTB1
ch_shape=size(ch);ch=reshape(ch,1,[]);
c_shape=size(c);c=reshape(c,1,[]);
wa_shape=size(wa);wa=reshape(wa,1,[]);
ifac_shape=size(ifac);ifac=reshape(ifac,1,[]);
%***FIRST EXECUTABLE STATEMENT  CFFTB1
nf = fix(ifac(2));
na = 0;
l1 = 1;
iw = 1;
for k1 = 1 : nf;
ip = fix(ifac(k1+2));
l2 = fix(ip.*l1);
ido = fix(fix(n./l2));
idot = fix(ido + ido);
idl1 = fix(idot.*l1);
if( ip==4 )
ix2 = fix(iw + idot);
ix3 = fix(ix2 + idot);
if( na~=0 )
[idot,l1,ch,c,dumvar5,dumvar6,dumvar7]=passb4(idot,l1,ch,c,wa(sub2ind(size(wa),max(iw,1)):end),wa(sub2ind(size(wa),max(ix2,1)):end),wa(sub2ind(size(wa),max(ix3,1)):end));   dumvar5i=find((wa(sub2ind(size(wa),max(iw,1)):end))~=(dumvar5));dumvar6i=find((wa(sub2ind(size(wa),max(ix2,1)):end))~=(dumvar6));dumvar7i=find((wa(sub2ind(size(wa),max(ix3,1)):end))~=(dumvar7));   wa(iw-1+dumvar5i)=dumvar5(dumvar5i); wa(ix2-1+dumvar6i)=dumvar6(dumvar6i); wa(ix3-1+dumvar7i)=dumvar7(dumvar7i); 
else;
[idot,l1,c,ch,dumvar5,dumvar6,dumvar7]=passb4(idot,l1,c,ch,wa(sub2ind(size(wa),max(iw,1)):end),wa(sub2ind(size(wa),max(ix2,1)):end),wa(sub2ind(size(wa),max(ix3,1)):end));   dumvar5i=find((wa(sub2ind(size(wa),max(iw,1)):end))~=(dumvar5));dumvar6i=find((wa(sub2ind(size(wa),max(ix2,1)):end))~=(dumvar6));dumvar7i=find((wa(sub2ind(size(wa),max(ix3,1)):end))~=(dumvar7));   wa(iw-1+dumvar5i)=dumvar5(dumvar5i); wa(ix2-1+dumvar6i)=dumvar6(dumvar6i); wa(ix3-1+dumvar7i)=dumvar7(dumvar7i); 
end;
na = fix(1 - na);
elseif( ip==2 ) ;
if( na~=0 )
[idot,l1,ch,c,wa(sub2ind(size(wa),max(iw,1)):end)]=passb2(idot,l1,ch,c,wa(sub2ind(size(wa),max(iw,1)):end));
else;
[idot,l1,c,ch,wa(sub2ind(size(wa),max(iw,1)):end)]=passb2(idot,l1,c,ch,wa(sub2ind(size(wa),max(iw,1)):end));
end;
na = fix(1 - na);
elseif( ip==3 ) ;
ix2 = fix(iw + idot);
if( na~=0 )
[idot,l1,ch,c,dumvar5,dumvar6]=passb3(idot,l1,ch,c,wa(sub2ind(size(wa),max(iw,1)):end),wa(sub2ind(size(wa),max(ix2,1)):end));   dumvar5i=find((wa(sub2ind(size(wa),max(iw,1)):end))~=(dumvar5));dumvar6i=find((wa(sub2ind(size(wa),max(ix2,1)):end))~=(dumvar6));   wa(iw-1+dumvar5i)=dumvar5(dumvar5i); wa(ix2-1+dumvar6i)=dumvar6(dumvar6i); 
else;
[idot,l1,c,ch,dumvar5,dumvar6]=passb3(idot,l1,c,ch,wa(sub2ind(size(wa),max(iw,1)):end),wa(sub2ind(size(wa),max(ix2,1)):end));   dumvar5i=find((wa(sub2ind(size(wa),max(iw,1)):end))~=(dumvar5));dumvar6i=find((wa(sub2ind(size(wa),max(ix2,1)):end))~=(dumvar6));   wa(iw-1+dumvar5i)=dumvar5(dumvar5i); wa(ix2-1+dumvar6i)=dumvar6(dumvar6i); 
end;
na = fix(1 - na);
elseif( ip~=5 ) ;
if( na~=0 )
ch_orig=ch; ch_orig=ch; c_orig=c;    [nac,idot,ip,l1,idl1,ch,dumvar7,dumvar8,c,dumvar10,wa(sub2ind(size(wa),max(iw,1)):end)]=passb(nac,idot,ip,l1,idl1,ch,ch,ch,c,c,wa(sub2ind(size(wa),max(iw,1)):end));    c(dumvar10~=c_orig)=dumvar10(dumvar10~=c_orig); ch(dumvar8~=ch_orig)=dumvar8(dumvar8~=ch_orig); ch(dumvar7~=ch_orig)=dumvar7(dumvar7~=ch_orig);
else;
c_orig=c; c_orig=c; ch_orig=ch;    [nac,idot,ip,l1,idl1,c,dumvar7,dumvar8,ch,dumvar10,wa(sub2ind(size(wa),max(iw,1)):end)]=passb(nac,idot,ip,l1,idl1,c,c,c,ch,ch,wa(sub2ind(size(wa),max(iw,1)):end));    ch(dumvar10~=ch_orig)=dumvar10(dumvar10~=ch_orig); c(dumvar8~=c_orig)=dumvar8(dumvar8~=c_orig); c(dumvar7~=c_orig)=dumvar7(dumvar7~=c_orig);
end;
if( nac~=0 )
na = fix(1 - na);
end;
else;
ix2 = fix(iw + idot);
ix3 = fix(ix2 + idot);
ix4 = fix(ix3 + idot);
if( na~=0 )
[idot,l1,ch,c,dumvar5,dumvar6,dumvar7,dumvar8]=passb5(idot,l1,ch,c,wa(sub2ind(size(wa),max(iw,1)):end),wa(sub2ind(size(wa),max(ix2,1)):end),wa(sub2ind(size(wa),max(ix3,1)):end),wa(sub2ind(size(wa),max(ix4,1)):end));   dumvar5i=find((wa(sub2ind(size(wa),max(iw,1)):end))~=(dumvar5));dumvar6i=find((wa(sub2ind(size(wa),max(ix2,1)):end))~=(dumvar6));dumvar7i=find((wa(sub2ind(size(wa),max(ix3,1)):end))~=(dumvar7));dumvar8i=find((wa(sub2ind(size(wa),max(ix4,1)):end))~=(dumvar8));   wa(iw-1+dumvar5i)=dumvar5(dumvar5i); wa(ix2-1+dumvar6i)=dumvar6(dumvar6i); wa(ix3-1+dumvar7i)=dumvar7(dumvar7i); wa(ix4-1+dumvar8i)=dumvar8(dumvar8i); 
else;
[idot,l1,c,ch,dumvar5,dumvar6,dumvar7,dumvar8]=passb5(idot,l1,c,ch,wa(sub2ind(size(wa),max(iw,1)):end),wa(sub2ind(size(wa),max(ix2,1)):end),wa(sub2ind(size(wa),max(ix3,1)):end),wa(sub2ind(size(wa),max(ix4,1)):end));   dumvar5i=find((wa(sub2ind(size(wa),max(iw,1)):end))~=(dumvar5));dumvar6i=find((wa(sub2ind(size(wa),max(ix2,1)):end))~=(dumvar6));dumvar7i=find((wa(sub2ind(size(wa),max(ix3,1)):end))~=(dumvar7));dumvar8i=find((wa(sub2ind(size(wa),max(ix4,1)):end))~=(dumvar8));   wa(iw-1+dumvar5i)=dumvar5(dumvar5i); wa(ix2-1+dumvar6i)=dumvar6(dumvar6i); wa(ix3-1+dumvar7i)=dumvar7(dumvar7i); wa(ix4-1+dumvar8i)=dumvar8(dumvar8i); 
end;
na = fix(1 - na);
end;
l1 = fix(l2);
iw = fix(iw +(ip-1).*idot);
end; k1 = fix(nf+1);
if( na==0 )
ch_shape=zeros(ch_shape);ch_shape(:)=ch(1:numel(ch_shape));ch=ch_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
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;
n2 = fix(n + n);
for i = 1 : n2;
c(i) = ch(i);
end; i = fix(n2+1);
ch_shape=zeros(ch_shape);ch_shape(:)=ch(1:numel(ch_shape));ch=ch_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
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
%DECK CFFTB