[ido,ip,l1,idl1,cc,c1,c2,ch,ch2,wa]=radfg(ido,ip,l1,idl1,cc,c1,c2,ch,ch2,wa);

function [ido,ip,l1,idl1,cc,c1,c2,ch,ch2,wa]=radfg(ido,ip,l1,idl1,cc,c1,c2,ch,ch2,wa);
persistent ai1 ai2 ar1 ar1h ar2 ar2h arg dc2 dcp ds2 dsp i ic idij idp2 ik ipp2 ipph is j j2 jc k l lc nbd tpi ; 

if isempty(ai1), ai1=0; end;
if isempty(ai2), ai2=0; end;
if isempty(ar1), ar1=0; end;
if isempty(ar1h), ar1h=0; end;
if isempty(ar2), ar2=0; end;
if isempty(ar2h), ar2h=0; end;
if isempty(arg), arg=0; end;
if isempty(dc2), dc2=0; end;
if isempty(dcp), dcp=0; end;
if isempty(ds2), ds2=0; end;
if isempty(dsp), dsp=0; end;
if isempty(tpi), tpi=0; end;
if isempty(i), i=0; end;
if isempty(ic), ic=0; end;
if isempty(idij), idij=0; end;
if isempty(idp2), idp2=0; end;
if isempty(ik), ik=0; end;
if isempty(ipp2), ipp2=0; end;
if isempty(ipph), ipph=0; end;
if isempty(is), is=0; end;
if isempty(j), j=0; end;
if isempty(j2), j2=0; end;
if isempty(jc), jc=0; end;
if isempty(k), k=0; end;
if isempty(l), l=0; end;
if isempty(lc), lc=0; end;
if isempty(nbd), nbd=0; end;
%***BEGIN PROLOGUE  RADFG
%***SUBSIDIARY
%***PURPOSE  Calculate the fast Fourier transform of subvectors of
%            arbitrary length.
%***LIBRARY   SLATEC (FFTPACK)
%***TYPE      SINGLE PRECISION (RADFG-S)
%***AUTHOR  Swarztrauber, P. N., (NCAR)
%***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)
%   890831  Modified array declarations.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900402  Added TYPE section.  (WRB)
%***end PROLOGUE  RADFG
ch_shape=size(ch);ch=reshape([ch(:).',zeros(1,ceil(numel(ch)./prod([ido,l1])).*prod([ido,l1])-numel(ch))],ido,l1,[]);
cc_shape=size(cc);cc=reshape([cc(:).',zeros(1,ceil(numel(cc)./prod([ido,ip])).*prod([ido,ip])-numel(cc))],ido,ip,[]);
c1_shape=size(c1);c1=reshape([c1(:).',zeros(1,ceil(numel(c1)./prod([ido,l1])).*prod([ido,l1])-numel(c1))],ido,l1,[]);
c2_shape=size(c2);c2=reshape([c2(:).',zeros(1,ceil(numel(c2)./prod([idl1])).*prod([idl1])-numel(c2))],idl1,[]);
ch2_shape=size(ch2);ch2=reshape([ch2(:).',zeros(1,ceil(numel(ch2)./prod([idl1])).*prod([idl1])-numel(ch2))],idl1,[]);
wa_shape=size(wa);wa=reshape(wa,1,[]);
%***FIRST EXECUTABLE STATEMENT  RADFG
tpi = 8..*atan(1.);
arg = tpi./ip;
dcp = cos(arg);
dsp = sin(arg);
ipph =fix(fix((ip+1)./2));
ipp2 = fix(ip + 2);
idp2 = fix(ido + 2);
nbd =fix(fix((ido-1)./2));
if( ido==1 )
for ik = 1 : idl1;
c2(ik,1) = ch2(ik,1);
end; ik = fix(idl1+1);
else;
for ik = 1 : idl1;
ch2(ik,1) = c2(ik,1);
end; ik = fix(idl1+1);
for j = 2 : ip;
for k = 1 : l1;
ch(1,k,j) = c1(1,k,j);
end; k = fix(l1+1);
end; j = fix(ip+1);
if( nbd>l1 )
is = fix(-ido);
for j = 2 : ip;
is = fix(is + ido);
for k = 1 : l1;
idij = fix(is);
%DIR$ IVDEP
for i = 3 : 2: ido ;
idij = fix(idij + 2);
ch(i-1,k,j) = wa(idij-1).*c1(i-1,k,j) + wa(idij).*c1(i,k,j);
ch(i,k,j) = wa(idij-1).*c1(i,k,j) - wa(idij).*c1(i-1,k,j);
end; i = fix(ido +1);
end; k = fix(l1+1);
end; j = fix(ip+1);
else;
is = fix(-ido);
for j = 2 : ip;
is = fix(is + ido);
idij = fix(is);
for i = 3 : 2: ido ;
idij = fix(idij + 2);
for k = 1 : l1;
ch(i-1,k,j) = wa(idij-1).*c1(i-1,k,j) + wa(idij).*c1(i,k,j);
ch(i,k,j) = wa(idij-1).*c1(i,k,j) - wa(idij).*c1(i-1,k,j);
end; k = fix(l1+1);
end; i = fix(ido +1);
end; j = fix(ip+1);
end;
if( nbd<l1 )
for j = 2 : ipph;
jc = fix(ipp2 - j);
for i = 3 : 2: ido ;
for k = 1 : l1;
c1(i-1,k,j) = ch(i-1,k,j) + ch(i-1,k,jc);
c1(i-1,k,jc) = ch(i,k,j) - ch(i,k,jc);
c1(i,k,j) = ch(i,k,j) + ch(i,k,jc);
c1(i,k,jc) = ch(i-1,k,jc) - ch(i-1,k,j);
end; k = fix(l1+1);
end; i = fix(ido +1);
end; j = fix(ipph+1);
else;
for j = 2 : ipph;
jc = fix(ipp2 - j);
for k = 1 : l1;
%DIR$ IVDEP
for i = 3 : 2: ido ;
c1(i-1,k,j) = ch(i-1,k,j) + ch(i-1,k,jc);
c1(i-1,k,jc) = ch(i,k,j) - ch(i,k,jc);
c1(i,k,j) = ch(i,k,j) + ch(i,k,jc);
c1(i,k,jc) = ch(i-1,k,jc) - ch(i-1,k,j);
end; i = fix(ido +1);
end; k = fix(l1+1);
end; j = fix(ipph+1);
end;
end;
for j = 2 : ipph;
jc = fix(ipp2 - j);
for k = 1 : l1;
c1(1,k,j) = ch(1,k,j) + ch(1,k,jc);
c1(1,k,jc) = ch(1,k,jc) - ch(1,k,j);
end; k = fix(l1+1);
end; j = fix(ipph+1);
%
ar1 = 1.;
ai1 = 0.;
for l = 2 : ipph;
lc = fix(ipp2 - l);
ar1h = dcp.*ar1 - dsp.*ai1;
ai1 = dcp.*ai1 + dsp.*ar1;
ar1 = ar1h;
for ik = 1 : idl1;
ch2(ik,l) = c2(ik,1) + ar1.*c2(ik,2);
ch2(ik,lc) = ai1.*c2(ik,ip);
end; ik = fix(idl1+1);
dc2 = ar1;
ds2 = ai1;
ar2 = ar1;
ai2 = ai1;
for j = 3 : ipph;
jc = fix(ipp2 - j);
ar2h = dc2.*ar2 - ds2.*ai2;
ai2 = dc2.*ai2 + ds2.*ar2;
ar2 = ar2h;
for ik = 1 : idl1;
ch2(ik,l) = ch2(ik,l) + ar2.*c2(ik,j);
ch2(ik,lc) = ch2(ik,lc) + ai2.*c2(ik,jc);
end; ik = fix(idl1+1);
end; j = fix(ipph+1);
end; l = fix(ipph+1);
for j = 2 : ipph;
for ik = 1 : idl1;
ch2(ik,1) = ch2(ik,1) + c2(ik,j);
end; ik = fix(idl1+1);
end; j = fix(ipph+1);
%
if( ido<l1 )
for i = 1 : ido;
for k = 1 : l1;
cc(i,1,k) = ch(i,k,1);
end; k = fix(l1+1);
end; i = fix(ido+1);
else;
for k = 1 : l1;
for i = 1 : ido;
cc(i,1,k) = ch(i,k,1);
end; i = fix(ido+1);
end; k = fix(l1+1);
end;
for j = 2 : ipph;
jc = fix(ipp2 - j);
j2 = fix(j + j);
for k = 1 : l1;
cc(ido,j2-2,k) = ch(1,k,j);
cc(1,j2-1,k) = ch(1,k,jc);
end; k = fix(l1+1);
end; j = fix(ipph+1);
if( ido==1 )
ch_shape=zeros(ch_shape);ch_shape(:)=ch(1:numel(ch_shape));ch=ch_shape;
cc_shape=zeros(cc_shape);cc_shape(:)=cc(1:numel(cc_shape));cc=cc_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ch2_shape=zeros(ch2_shape);ch2_shape(:)=ch2(1:numel(ch2_shape));ch2=ch2_shape;
wa_shape=zeros(wa_shape);wa_shape(:)=wa(1:numel(wa_shape));wa=wa_shape;
return;
end;
if( nbd<l1 )
for j = 2 : ipph;
jc = fix(ipp2 - j);
j2 = fix(j + j);
for i = 3 : 2: ido ;
ic = fix(idp2 - i);
for k = 1 : l1;
cc(i-1,j2-1,k) = ch(i-1,k,j) + ch(i-1,k,jc);
cc(ic-1,j2-2,k) = ch(i-1,k,j) - ch(i-1,k,jc);
cc(i,j2-1,k) = ch(i,k,j) + ch(i,k,jc);
cc(ic,j2-2,k) = ch(i,k,jc) - ch(i,k,j);
end; k = fix(l1+1);
end; i = fix(ido +1);
end; j = fix(ipph+1);
else;
for j = 2 : ipph;
jc = fix(ipp2 - j);
j2 = fix(j + j);
for k = 1 : l1;
%DIR$ IVDEP
for i = 3 : 2: ido ;
ic = fix(idp2 - i);
cc(i-1,j2-1,k) = ch(i-1,k,j) + ch(i-1,k,jc);
cc(ic-1,j2-2,k) = ch(i-1,k,j) - ch(i-1,k,jc);
cc(i,j2-1,k) = ch(i,k,j) + ch(i,k,jc);
cc(ic,j2-2,k) = ch(i,k,jc) - ch(i,k,j);
end; i = fix(ido +1);
end; k = fix(l1+1);
end; j = fix(ipph+1);
ch_shape=zeros(ch_shape);ch_shape(:)=ch(1:numel(ch_shape));ch=ch_shape;
cc_shape=zeros(cc_shape);cc_shape(:)=cc(1:numel(cc_shape));cc=cc_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ch2_shape=zeros(ch2_shape);ch2_shape(:)=ch2(1:numel(ch2_shape));ch2=ch2_shape;
wa_shape=zeros(wa_shape);wa_shape(:)=wa(1:numel(wa_shape));wa=wa_shape;
return;
end;
ch_shape=zeros(ch_shape);ch_shape(:)=ch(1:numel(ch_shape));ch=ch_shape;
cc_shape=zeros(cc_shape);cc_shape(:)=cc(1:numel(cc_shape));cc=cc_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ch2_shape=zeros(ch2_shape);ch2_shape(:)=ch2(1:numel(ch2_shape));ch2=ch2_shape;
wa_shape=zeros(wa_shape);wa_shape(:)=wa(1:numel(wa_shape));wa=wa_shape;
end
%DECK RAND