[n,x,wsave]=cost(n,x,wsave);

function [n,x,wsave]=cost(n,x,wsave);
persistent c1 i k kc modn nm1 np1 ns2 t1 t2 tx2 x1h x1p3 xi xim2 ; 

if isempty(c1), c1=0; end;
if isempty(t1), t1=0; end;
if isempty(t2), t2=0; end;
if isempty(tx2), tx2=0; end;
if isempty(x1h), x1h=0; end;
if isempty(x1p3), x1p3=0; end;
if isempty(xi), xi=0; end;
if isempty(xim2), xim2=0; end;
if isempty(i), i=0; end;
if isempty(k), k=0; end;
if isempty(kc), kc=0; end;
if isempty(modn), modn=0; end;
if isempty(nm1), nm1=0; end;
if isempty(np1), np1=0; end;
if isempty(ns2), ns2=0; end;
%***BEGIN PROLOGUE  COST
%***PURPOSE  Compute the cosine transform of a real, even sequence.
%***LIBRARY   SLATEC (FFTPACK)
%***CATEGORY  J1A3
%***TYPE      SINGLE PRECISION (COST-S)
%***KEYWORDS  COSINE FOURIER TRANSFORM, FFTPACK
%***AUTHOR  Swarztrauber, P. N., (NCAR)
%***DESCRIPTION
%
%  subroutine COST computes the discrete Fourier cosine transform
%  of an even sequence X(I).  The transform is defined below at output
%  parameter X.
%
%  COST is the unnormalized inverse of itself since a call of COST
%  followed by another call of COST will multiply the input sequence
%  X by 2*(N-1).  The transform is defined below at output parameter X.
%
%  The array WSAVE which is used by subroutine COST must be
%  initialized by calling subroutine COSTI(N,WSAVE).
%
%  Input Parameters
%
%  N       the length of the sequence X.  N must be greater than 1.
%          The method is most efficient when N-1 is a product of
%          small primes.
%
%  X       an array which contains the sequence to be transformed
%
%  WSAVE   a work array which must be dimensioned at least 3*N+15
%          in the program that calls COST.  The WSAVE array must be
%          initialized by calling subroutine COSTI(N,WSAVE), and a
%          different WSAVE array 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.
%
%  Output Parameters
%
%  X       For I=1,...,N
%
%             X(I) = X(1)+(-1)**(I-1)*X(N)
%
%               + the sum from K=2 to K=N-1
%
%                 2*X(K)*COS((K-1)*(I-1)*PI/(N-1))
%
%               A call of COST followed by another call of
%               COST will multiply the sequence X by 2*(N-1).
%               Hence COST is the unnormalized inverse
%               of itself.
%
%  WSAVE   contains initialization calculations which must not be
%          destroyed between calls of COST.
%
%***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
%                 Computations (G. Rodrigue, ed.), Academic Press,
%                 1982, pp. 51-83.
%***ROUTINES CALLED  RFFTF
%***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 (*)
%   861211  REVISION DATE from Version 3.2
%   881128  Modified by Dick Valent to meet prologue standards.
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  COST
x_shape=size(x);x=reshape(x,1,[]);
wsave_shape=size(wsave);wsave=reshape(wsave,1,[]);
%***FIRST EXECUTABLE STATEMENT  COST
nm1 = fix(n - 1);
np1 = fix(n + 1);
ns2 = fix(fix(n./2));
if( n<2 )
elseif( n==2 ) ;
x1h = x(1) + x(2);
x(2) = x(1) - x(2);
x(1) = x1h;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
wsave_shape=zeros(wsave_shape);wsave_shape(:)=wsave(1:numel(wsave_shape));wsave=wsave_shape;
return;
elseif( n>3 ) ;
c1 = x(1) - x(n);
x(1) = x(1) + x(n);
for k = 2 : ns2;
kc = fix(np1 - k);
t1 = x(k) + x(kc);
t2 = x(k) - x(kc);
c1 = c1 + wsave(kc).*t2;
t2 = wsave(k).*t2;
x(k) = t1 - t2;
x(kc) = t1 + t2;
end; k = fix(ns2+1);
modn = fix(rem(n,2));
if( modn~=0 )
x(ns2+1) = x(ns2+1) + x(ns2+1);
end;
[nm1,x,wsave(sub2ind(size(wsave),max(n+1,1)):end)]=rfftf(nm1,x,wsave(sub2ind(size(wsave),max(n+1,1)):end));
xim2 = x(2);
x(2) = c1;
for i = 4 : 2: n ;
xi = x(i);
x(i) = x(i-2) - x(i-1);
x(i-1) = xim2;
xim2 = xi;
end; i = fix(n +1);
if( modn~=0 )
x(n) = xim2;
end;
else;
x1p3 = x(1) + x(3);
tx2 = x(2) + x(2);
x(2) = x(1) - x(3);
x(1) = x1p3 + tx2;
x(3) = x1p3 - tx2;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
wsave_shape=zeros(wsave_shape);wsave_shape(:)=wsave(1:numel(wsave_shape));wsave=wsave_shape;
return;
end;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
wsave_shape=zeros(wsave_shape);wsave_shape(:)=wsave(1:numel(wsave_shape));wsave=wsave_shape;
end
%DECK COSTI