subroutine dqcheb(x,fval,cheb12,cheb24) c***begin prologue dqcheb c***refer to dqc25c,dqc25f,dqc25s c***routines called (none) c***revision date 830518 (yymmdd) c***keywords chebyshev series expansion, fast fourier transform c***author piessens,robert,appl. math. & progr. div. - k.u.leuven c de doncker,elise,appl. math. & progr. div. - k.u.leuven c***purpose this routine computes the chebyshev series expansion c of degrees 12 and 24 of a function using a c fast fourier transform method c f(x) = sum(k=1,..,13) (cheb12(k)*t(k-1,x)), c f(x) = sum(k=1,..,25) (cheb24(k)*t(k-1,x)), c where t(k,x) is the chebyshev polynomial of degree k. c***description c c chebyshev series expansion c standard fortran subroutine c double precision version c c parameters c on entry c x - double precision c vector of dimension 11 containing the c values cos(k*pi/24), k = 1, ..., 11 c c fval - double precision c vector of dimension 25 containing the c function values at the points c (b+a+(b-a)*cos(k*pi/24))/2, k = 0, ...,24, c where (a,b) is the approximation interval. c fval(1) and fval(25) are divided by two c (these values are destroyed at output). c c on return c cheb12 - double precision c vector of dimension 13 containing the c chebyshev coefficients for degree 12 c c cheb24 - double precision c vector of dimension 25 containing the c chebyshev coefficients for degree 24 c c***end prologue dqcheb c double precision alam,alam1,alam2,cheb12,cheb24,fval,part1,part2, * part3,v,x integer i,j c dimension cheb12(13),cheb24(25),fval(25),v(12),x(11) c c***first executable statement dqcheb do 10 i=1,12 j = 26-i v(i) = fval(i)-fval(j) fval(i) = fval(i)+fval(j) 10 continue 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 do 20 i=1,6 j = 14-i v(i) = fval(i)-fval(j) fval(i) = fval(i)+fval(j) 20 continue 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 do 30 i=1,3 j = 8-i v(i) = fval(i)-fval(j) fval(i) = fval(i)+fval(j) 30 continue 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 do 40 i=2,12 cheb12(i) = cheb12(i)*alam 40 continue alam = 0.5d+00*alam cheb12(1) = cheb12(1)*alam cheb12(13) = cheb12(13)*alam do 50 i=2,24 cheb24(i) = cheb24(i)*alam 50 continue cheb24(1) = 0.5d+00*alam*cheb24(1) cheb24(25) = 0.5d+00*alam*cheb24(25) return end