[n,wsave]=costi(n,wsave);

function [n,wsave]=costi(n,wsave);
persistent dt fk k kc nm1 np1 ns2 pi ; 

if isempty(dt), dt=0; end;
if isempty(fk), fk=0; end;
if isempty(pi), pi=0; end;
if isempty(k), k=0; end;
if isempty(kc), kc=0; end;
if isempty(nm1), nm1=0; end;
if isempty(np1), np1=0; end;
if isempty(ns2), ns2=0; end;
%***BEGIN PROLOGUE  COSTI
%***PURPOSE  Initialize a work array for COST.
%***LIBRARY   SLATEC (FFTPACK)
%***CATEGORY  J1A3
%***TYPE      SINGLE PRECISION (COSTI-S)
%***KEYWORDS  COSINE FOURIER TRANSFORM, FFTPACK
%***AUTHOR  Swarztrauber, P. N., (NCAR)
%***DESCRIPTION
%
%  subroutine COSTI initializes the array WSAVE which is used in
%  subroutine COST.  The prime factorization of N together with
%  a tabulation of the trigonometric functions are computed and
%  stored in WSAVE.
%
%  Input Parameter
%
%  N       the length of the sequence to be transformed.  The method
%          is most efficient when N-1 is a product of small primes.
%
%  Output Parameter
%
%  WSAVE   a work array which must be dimensioned at least 3*N+15.
%          Different WSAVE arrays are required for different values
%          of N.  The contents of WSAVE must not be changed 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  RFFTI
%***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 PI 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)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  COSTI
wsave_shape=size(wsave);wsave=reshape(wsave,1,[]);
%***FIRST EXECUTABLE STATEMENT  COSTI
if( n<=3 )
wsave_shape=zeros(wsave_shape);wsave_shape(:)=wsave(1:numel(wsave_shape));wsave=wsave_shape;
return;
end;
pi = 4..*atan(1.);
nm1 = fix(n - 1);
np1 = fix(n + 1);
ns2 = fix(fix(n./2));
dt = pi./nm1;
fk = 0.;
for k = 2 : ns2;
kc = fix(np1 - k);
fk = fk + 1.;
wsave(k) = 2..*sin(fk.*dt);
wsave(kc) = 2..*cos(fk.*dt);
end; k = fix(ns2+1);
[nm1,wsave(sub2ind(size(wsave),max(n+1,1)):end)]=rffti(nm1,wsave(sub2ind(size(wsave),max(n+1,1)):end));
wsave_shape=zeros(wsave_shape);wsave_shape(:)=wsave(1:numel(wsave_shape));wsave=wsave_shape;
end
%DECK COT