[n,x,w,xh]=cosqf1(n,x,w,xh);

function [n,x,w,xh]=cosqf1(n,x,w,xh);
persistent i k kc modn np2 ns2 xim1 ; 

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(np2), np2=0; end;
if isempty(ns2), ns2=0; end;
if isempty(xim1), xim1=0; end;
%***BEGIN PROLOGUE  COSQF1
%***SUBSIDIARY
%***PURPOSE  Compute the forward cosine transform with odd wave numbers.
%***LIBRARY   SLATEC (FFTPACK)
%***CATEGORY  J1A3
%***TYPE      SINGLE PRECISION (COSQF1-S)
%***KEYWORDS  FFTPACK, FOURIER TRANSFORM
%***AUTHOR  Swarztrauber, P. N., (NCAR)
%***DESCRIPTION
%
%  subroutine COSQF1 computes the fast Fourier transform of quarter
%  wave data. That is, COSQF1 computes the coefficients in a cosine
%  series representation with only odd wave numbers.  The transform
%  is defined below at Output Parameter X
%
%***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 (*).
%   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  COSQF1
x_shape=size(x);x=reshape(x,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
xh_shape=size(xh);xh=reshape(xh,1,[]);
%***FIRST EXECUTABLE STATEMENT  COSQF1
ns2 =fix(fix((n+1)./2));
np2 = fix(n + 2);
for k = 2 : ns2;
kc = fix(np2 - k);
xh(k) = x(k) + x(kc);
xh(kc) = x(k) - x(kc);
end; k = fix(ns2+1);
modn = fix(rem(n,2));
if( modn==0 )
xh(ns2+1) = x(ns2+1) + x(ns2+1);
end;
for k = 2 : ns2;
kc = fix(np2 - k);
x(k) = w(k-1).*xh(kc) + w(kc-1).*xh(k);
x(kc) = w(k-1).*xh(k) - w(kc-1).*xh(kc);
end; k = fix(ns2+1);
if( modn==0 )
x(ns2+1) = w(ns2).*xh(ns2+1);
end;
[n,x,xh]=rfftf(n,x,xh);
for i = 3 : 2: n ;
xim1 = x(i-1) - x(i);
x(i) = x(i-1) + x(i);
x(i-1) = xim1;
end; i = fix(n +1);
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
xh_shape=zeros(xh_shape);xh_shape(:)=xh(1:numel(xh_shape));xh=xh_shape;
end
%DECK COSQF