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

function [n,x,w,xh]=cosqb1(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  COSQB1
%***SUBSIDIARY
%***PURPOSE  Compute the unnormalized inverse of COSQF1.
%***LIBRARY   SLATEC (FFTPACK)
%***CATEGORY  J1A3
%***TYPE      SINGLE PRECISION (COSQB1-S)
%***KEYWORDS  FFTPACK, FOURIER TRANSFORM
%***AUTHOR  Swarztrauber, P. N., (NCAR)
%***DESCRIPTION
%
%  subroutine COSQB1 computes the fast Fourier transform of quarter
%  wave data. That is, COSQB1 computes a sequence from its
%  representation in terms of a cosine series with 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  RFFTB
%***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  COSQB1
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  COSQB1
ns2 =fix(fix((n+1)./2));
np2 = fix(n + 2);
for i = 3 : 2: n ;
xim1 = x(i-1) + x(i);
x(i) = x(i) - x(i-1);
x(i-1) = xim1;
end; i = fix(n +1);
x(1) = x(1) + x(1);
modn = fix(rem(n,2));
if( modn==0 )
x(n) = x(n) + x(n);
end;
[n,x,xh]=rfftb(n,x,xh);
for k = 2 : ns2;
kc = fix(np2 - k);
xh(k) = w(k-1).*x(kc) + w(kc-1).*x(k);
xh(kc) = w(k-1).*x(k) - w(kc-1).*x(kc);
end; k = fix(ns2+1);
if( modn==0 )
x(ns2+1) = w(ns2).*(x(ns2+1)+x(ns2+1));
end;
for k = 2 : ns2;
kc = fix(np2 - k);
x(k) = xh(k) + xh(kc);
x(kc) = xh(k) - xh(kc);
end; k = fix(ns2+1);
x(1) = x(1) + x(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 COSQB