[x,fval,cheb12,cheb24]=dqcheb(x,fval,cheb12,cheb24);

function [x,fval,cheb12,cheb24]=dqcheb(x,fval,cheb12,cheb24);
%***BEGIN PROLOGUE  DQCHEB
%***SUBSIDIARY
%***PURPOSE  This routine computes the CHEBYSHEV series expansion
%            of degrees 12 and 24 of a function using A
%            FAST FOURIER TRANSFORM METHOD
%            F(X) = SUM(K=1,..,13) (CHEB12(K)*T(K-1,X)),
%            F(X) = SUM(K=1,..,25) (CHEB24(K)*T(K-1,X)),
%            Where T(K,X) is the CHEBYSHEV POLYNOMIAL OF DEGREE K.
%***LIBRARY   SLATEC
%***TYPE      doubleprecision (QCHEB-S, DQCHEB-D)
%***KEYWORDS  CHEBYSHEV SERIES EXPANSION, FAST FOURIER TRANSFORM
%***AUTHOR  Piessens, Robert
%             Applied Mathematics and Programming Division
%             K. U. Leuven
%           de Doncker, Elise
%             Applied Mathematics and Programming Division
%             K. U. Leuven
%***DESCRIPTION
%
%        Chebyshev Series Expansion
%        Standard Fortran Subroutine
%        doubleprecision version
%
%        PARAMETERS
%          ON ENTRY
%           X      - doubleprecision
%                    Vector of dimension 11 containing the
%                    Values COS(K*PI/24), K = 1, ..., 11
%
%           FVAL   - doubleprecision
%                    Vector of dimension 25 containing the
%                    function values at the points
%                    (B+A+(B-A)*COS(K*PI/24))/2, K = 0, ...,24,
%                    where (A,B) is the approximation interval.
%                    FVAL(1) and FVAL(25) are divided by two
%                    (these values are destroyed at output).
%
%          ON RETURN
%           CHEB12 - doubleprecision
%                    Vector of dimension 13 containing the
%                    CHEBYSHEV coefficients for degree 12
%
%           CHEB24 - doubleprecision
%                    Vector of dimension 25 containing the
%                    CHEBYSHEV Coefficients for degree 24
%
%***SEE ALSO  DQC25C, DQC25F, DQC25S
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   810101  DATE WRITTEN
%   830518  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%***end PROLOGUE  DQCHEB
%
persistent alam alam1 alam2 i j part1 part2 part3 v ; 

if isempty(alam), alam=0; end;
if isempty(alam1), alam1=0; end;
if isempty(alam2), alam2=0; end;
if isempty(part1), part1=0; end;
if isempty(part2), part2=0; end;
if isempty(part3), part3=0; end;
if isempty(v), v=zeros(1,12); end;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
%
%
%***FIRST EXECUTABLE STATEMENT  DQCHEB
for i = 1 : 12;
j = fix(26 - i);
v(i) = fval(i) - fval(j);
fval(i) = fval(i) + fval(j);
end; i = fix(12+1);
alam1 = v(1) - v(9);
alam2 = x(6).*(v(3)-v(7)-v(11));
cheb12(4) = alam1 + alam2;
cheb12(10) = alam1 - alam2;
alam1 = v(2) - v(8) - v(10);
alam2 = v(4) - v(6) - v(12);
alam = x(3).*alam1 + x(9).*alam2;
cheb24(4) = cheb12(4) + alam;
cheb24(22) = cheb12(4) - alam;
alam = x(9).*alam1 - x(3).*alam2;
cheb24(10) = cheb12(10) + alam;
cheb24(16) = cheb12(10) - alam;
part1 = x(4).*v(5);
part2 = x(8).*v(9);
part3 = x(6).*v(7);
alam1 = v(1) + part1 + part2;
alam2 = x(2).*v(3) + part3 + x(10).*v(11);
cheb12(2) = alam1 + alam2;
cheb12(12) = alam1 - alam2;
alam = x(1).*v(2) + x(3).*v(4) + x(5).*v(6) + x(7).*v(8) + x(9).*v(10)+ x(11).*v(12);
cheb24(2) = cheb12(2) + alam;
cheb24(24) = cheb12(2) - alam;
alam = x(11).*v(2) - x(9).*v(4) + x(7).*v(6) - x(5).*v(8) + x(3).*v(10)- x(1).*v(12);
cheb24(12) = cheb12(12) + alam;
cheb24(14) = cheb12(12) - alam;
alam1 = v(1) - part1 + part2;
alam2 = x(10).*v(3) - part3 + x(2).*v(11);
cheb12(6) = alam1 + alam2;
cheb12(8) = alam1 - alam2;
alam = x(5).*v(2) - x(9).*v(4) - x(1).*v(6) - x(11).*v(8) + x(3).*v(10)+ x(7).*v(12);
cheb24(6) = cheb12(6) + alam;
cheb24(20) = cheb12(6) - alam;
alam = x(7).*v(2) - x(3).*v(4) - x(11).*v(6) + x(1).*v(8) - x(9).*v(10)- x(5).*v(12);
cheb24(8) = cheb12(8) + alam;
cheb24(18) = cheb12(8) - alam;
for i = 1 : 6;
j = fix(14 - i);
v(i) = fval(i) - fval(j);
fval(i) = fval(i) + fval(j);
end; i = fix(6+1);
alam1 = v(1) + x(8).*v(5);
alam2 = x(4).*v(3);
cheb12(3) = alam1 + alam2;
cheb12(11) = alam1 - alam2;
cheb12(7) = v(1) - v(5);
alam = x(2).*v(2) + x(6).*v(4) + x(10).*v(6);
cheb24(3) = cheb12(3) + alam;
cheb24(23) = cheb12(3) - alam;
alam = x(6).*(v(2)-v(4)-v(6));
cheb24(7) = cheb12(7) + alam;
cheb24(19) = cheb12(7) - alam;
alam = x(10).*v(2) - x(6).*v(4) + x(2).*v(6);
cheb24(11) = cheb12(11) + alam;
cheb24(15) = cheb12(11) - alam;
for i = 1 : 3;
j = fix(8 - i);
v(i) = fval(i) - fval(j);
fval(i) = fval(i) + fval(j);
end; i = fix(3+1);
cheb12(5) = v(1) + x(8).*v(3);
cheb12(9) = fval(1) - x(8).*fval(3);
alam = x(4).*v(2);
cheb24(5) = cheb12(5) + alam;
cheb24(21) = cheb12(5) - alam;
alam = x(8).*fval(2) - fval(4);
cheb24(9) = cheb12(9) + alam;
cheb24(17) = cheb12(9) - alam;
cheb12(1) = fval(1) + fval(3);
alam = fval(2) + fval(4);
cheb24(1) = cheb12(1) + alam;
cheb24(25) = cheb12(1) - alam;
cheb12(13) = v(1) - v(3);
cheb24(13) = cheb12(13);
alam = 0.1d+01./0.6d+01;
for i = 2 : 12;
cheb12(i) = cheb12(i).*alam;
end; i = fix(12+1);
alam = 0.5d+00.*alam;
cheb12(1) = cheb12(1).*alam;
cheb12(13) = cheb12(13).*alam;
for i = 2 : 24;
cheb24(i) = cheb24(i).*alam;
end; i = fix(24+1);
cheb24(1) = 0.5d+00.*alam.*cheb24(1);
cheb24(25) = 0.5d+00.*alam.*cheb24(25);
end
%DECK DQDOTA