[ido,l1,cc,ch,wa1]=radf2(ido,l1,cc,ch,wa1);

function [ido,l1,cc,ch,wa1]=radf2(ido,l1,cc,ch,wa1);
persistent i ic idp2 k ti2 tr2 ; 

if isempty(ti2), ti2=0; end;
if isempty(tr2), tr2=0; end;
if isempty(i), i=0; end;
if isempty(ic), ic=0; end;
if isempty(idp2), idp2=0; end;
if isempty(k), k=0; end;
%***BEGIN PROLOGUE  RADF2
%***SUBSIDIARY
%***PURPOSE  Calculate the fast Fourier transform of subvectors of
%            length two.
%***LIBRARY   SLATEC (FFTPACK)
%***TYPE      SINGLE PRECISION (RADF2-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
%           changing dummy array size declarations (1) to (*).
%   881128  Modified by Dick Valent to meet prologue standards.
%   890831  Modified array declarations.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900402  Added TYPE section.  (WRB)
%***end PROLOGUE  RADF2
ch_shape=size(ch);ch=reshape([ch(:).',zeros(1,ceil(numel(ch)./prod([ido,2])).*prod([ido,2])-numel(ch))],ido,2,[]);
cc_orig=cc;cc_shape=[ido,l1,2];cc=reshape([cc_orig(1:min(prod(cc_shape),numel(cc_orig))),zeros(1,max(0,prod(cc_shape)-numel(cc_orig)))],cc_shape);
wa1_shape=size(wa1);wa1=reshape(wa1,1,[]);
%***FIRST EXECUTABLE STATEMENT  RADF2
for k = 1 : l1;
ch(1,1,k) = cc(1,k,1) + cc(1,k,2);
ch(ido,2,k) = cc(1,k,1) - cc(1,k,2);
end; k = fix(l1+1);
if( ido<2 )
ch_shape=zeros(ch_shape);ch_shape(:)=ch(1:numel(ch_shape));ch=ch_shape;
cc_orig(1:prod(cc_shape))=cc;cc=cc_orig;
wa1_shape=zeros(wa1_shape);wa1_shape(:)=wa1(1:numel(wa1_shape));wa1=wa1_shape;
return;
end;
if( ido~=2 )
idp2 = fix(ido + 2);
if(fix((ido-1)./2)<l1 )
for i = 3 : 2: ido ;
ic = fix(idp2 - i);
%DIR$ IVDEP
for k = 1 : l1;
tr2 = wa1(i-2).*cc(i-1,k,2) + wa1(i-1).*cc(i,k,2);
ti2 = wa1(i-2).*cc(i,k,2) - wa1(i-1).*cc(i-1,k,2);
ch(i,1,k) = cc(i,k,1) + ti2;
ch(ic,2,k) = ti2 - cc(i,k,1);
ch(i-1,1,k) = cc(i-1,k,1) + tr2;
ch(ic-1,2,k) = cc(i-1,k,1) - tr2;
end; k = fix(l1+1);
end; i = fix(ido +1);
else;
for k = 1 : l1;
%DIR$ IVDEP
for i = 3 : 2: ido ;
ic = fix(idp2 - i);
tr2 = wa1(i-2).*cc(i-1,k,2) + wa1(i-1).*cc(i,k,2);
ti2 = wa1(i-2).*cc(i,k,2) - wa1(i-1).*cc(i-1,k,2);
ch(i,1,k) = cc(i,k,1) + ti2;
ch(ic,2,k) = ti2 - cc(i,k,1);
ch(i-1,1,k) = cc(i-1,k,1) + tr2;
ch(ic-1,2,k) = cc(i-1,k,1) - tr2;
end; i = fix(ido +1);
end; k = fix(l1+1);
end;
if( rem(ido,2)==1 )
ch_shape=zeros(ch_shape);ch_shape(:)=ch(1:numel(ch_shape));ch=ch_shape;
cc_orig(1:prod(cc_shape))=cc;cc=cc_orig;
wa1_shape=zeros(wa1_shape);wa1_shape(:)=wa1(1:numel(wa1_shape));wa1=wa1_shape;
return;
end;
end;
for k = 1 : l1;
ch(1,2,k) = -cc(ido,k,2);
ch(ido,1,k) = cc(ido,k,1);
end; k = fix(l1+1);
ch_shape=zeros(ch_shape);ch_shape(:)=ch(1:numel(ch_shape));ch=ch_shape;
cc_orig(1:prod(cc_shape))=cc;cc=cc_orig;
wa1_shape=zeros(wa1_shape);wa1_shape(:)=wa1(1:numel(wa1_shape));wa1=wa1_shape;
end
%DECK RADF3