qcheb.f 5.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160
  1. *DECK QCHEB
  2. SUBROUTINE QCHEB (X, FVAL, CHEB12, CHEB24)
  3. C***BEGIN PROLOGUE QCHEB
  4. C***SUBSIDIARY
  5. C***PURPOSE This routine computes the CHEBYSHEV series expansion
  6. C of degrees 12 and 24 of a function using A
  7. C FAST FOURIER TRANSFORM METHOD
  8. C F(X) = SUM(K=1,..,13) (CHEB12(K)*T(K-1,X)),
  9. C F(X) = SUM(K=1,..,25) (CHEB24(K)*T(K-1,X)),
  10. C Where T(K,X) is the CHEBYSHEV POLYNOMIAL OF DEGREE K.
  11. C***LIBRARY SLATEC
  12. C***TYPE SINGLE PRECISION (QCHEB-S, DQCHEB-D)
  13. C***KEYWORDS CHEBYSHEV SERIES EXPANSION, FAST FOURIER TRANSFORM
  14. C***AUTHOR Piessens, Robert
  15. C Applied Mathematics and Programming Division
  16. C K. U. Leuven
  17. C de Doncker, Elise
  18. C Applied Mathematics and Programming Division
  19. C K. U. Leuven
  20. C***DESCRIPTION
  21. C
  22. C Chebyshev Series Expansion
  23. C Standard Fortran Subroutine
  24. C Real version
  25. C
  26. C PARAMETERS
  27. C ON ENTRY
  28. C X - Real
  29. C Vector of dimension 11 containing the
  30. C Values COS(K*PI/24), K = 1, ..., 11
  31. C
  32. C FVAL - Real
  33. C Vector of dimension 25 containing the
  34. C function values at the points
  35. C (B+A+(B-A)*COS(K*PI/24))/2, K = 0, ...,24,
  36. C where (A,B) is the approximation interval.
  37. C FVAL(1) and FVAL(25) are divided by two
  38. C (these values are destroyed at output).
  39. C
  40. C ON RETURN
  41. C CHEB12 - Real
  42. C Vector of dimension 13 containing the
  43. C CHEBYSHEV coefficients for degree 12
  44. C
  45. C CHEB24 - Real
  46. C Vector of dimension 25 containing the
  47. C CHEBYSHEV Coefficients for degree 24
  48. C
  49. C***SEE ALSO QC25C, QC25F, QC25S
  50. C***ROUTINES CALLED (NONE)
  51. C***REVISION HISTORY (YYMMDD)
  52. C 810101 DATE WRITTEN
  53. C 830518 REVISION DATE from Version 3.2
  54. C 891214 Prologue converted to Version 4.0 format. (BAB)
  55. C 900328 Added TYPE section. (WRB)
  56. C***END PROLOGUE QCHEB
  57. C
  58. REAL ALAM,ALAM1,ALAM2,CHEB12,CHEB24,
  59. 1 FVAL,PART1,PART2,PART3,V,X
  60. INTEGER I,J
  61. C
  62. DIMENSION CHEB12(13),CHEB24(25),FVAL(25),V(12),X(11)
  63. C
  64. C***FIRST EXECUTABLE STATEMENT QCHEB
  65. DO 10 I=1,12
  66. J = 26-I
  67. V(I) = FVAL(I)-FVAL(J)
  68. FVAL(I) = FVAL(I)+FVAL(J)
  69. 10 CONTINUE
  70. ALAM1 = V(1)-V(9)
  71. ALAM2 = X(6)*(V(3)-V(7)-V(11))
  72. CHEB12(4) = ALAM1+ALAM2
  73. CHEB12(10) = ALAM1-ALAM2
  74. ALAM1 = V(2)-V(8)-V(10)
  75. ALAM2 = V(4)-V(6)-V(12)
  76. ALAM = X(3)*ALAM1+X(9)*ALAM2
  77. CHEB24(4) = CHEB12(4)+ALAM
  78. CHEB24(22) = CHEB12(4)-ALAM
  79. ALAM = X(9)*ALAM1-X(3)*ALAM2
  80. CHEB24(10) = CHEB12(10)+ALAM
  81. CHEB24(16) = CHEB12(10)-ALAM
  82. PART1 = X(4)*V(5)
  83. PART2 = X(8)*V(9)
  84. PART3 = X(6)*V(7)
  85. ALAM1 = V(1)+PART1+PART2
  86. ALAM2 = X(2)*V(3)+PART3+X(10)*V(11)
  87. CHEB12(2) = ALAM1+ALAM2
  88. CHEB12(12) = ALAM1-ALAM2
  89. ALAM = X(1)*V(2)+X(3)*V(4)+X(5)*V(6)+X(7)*V(8)
  90. 1 +X(9)*V(10)+X(11)*V(12)
  91. CHEB24(2) = CHEB12(2)+ALAM
  92. CHEB24(24) = CHEB12(2)-ALAM
  93. ALAM = X(11)*V(2)-X(9)*V(4)+X(7)*V(6)-X(5)*V(8)
  94. 1 +X(3)*V(10)-X(1)*V(12)
  95. CHEB24(12) = CHEB12(12)+ALAM
  96. CHEB24(14) = CHEB12(12)-ALAM
  97. ALAM1 = V(1)-PART1+PART2
  98. ALAM2 = X(10)*V(3)-PART3+X(2)*V(11)
  99. CHEB12(6) = ALAM1+ALAM2
  100. CHEB12(8) = ALAM1-ALAM2
  101. ALAM = X(5)*V(2)-X(9)*V(4)-X(1)*V(6)
  102. 1 -X(11)*V(8)+X(3)*V(10)+X(7)*V(12)
  103. CHEB24(6) = CHEB12(6)+ALAM
  104. CHEB24(20) = CHEB12(6)-ALAM
  105. ALAM = X(7)*V(2)-X(3)*V(4)-X(11)*V(6)+X(1)*V(8)
  106. 1 -X(9)*V(10)-X(5)*V(12)
  107. CHEB24(8) = CHEB12(8)+ALAM
  108. CHEB24(18) = CHEB12(8)-ALAM
  109. DO 20 I=1,6
  110. J = 14-I
  111. V(I) = FVAL(I)-FVAL(J)
  112. FVAL(I) = FVAL(I)+FVAL(J)
  113. 20 CONTINUE
  114. ALAM1 = V(1)+X(8)*V(5)
  115. ALAM2 = X(4)*V(3)
  116. CHEB12(3) = ALAM1+ALAM2
  117. CHEB12(11) = ALAM1-ALAM2
  118. CHEB12(7) = V(1)-V(5)
  119. ALAM = X(2)*V(2)+X(6)*V(4)+X(10)*V(6)
  120. CHEB24(3) = CHEB12(3)+ALAM
  121. CHEB24(23) = CHEB12(3)-ALAM
  122. ALAM = X(6)*(V(2)-V(4)-V(6))
  123. CHEB24(7) = CHEB12(7)+ALAM
  124. CHEB24(19) = CHEB12(7)-ALAM
  125. ALAM = X(10)*V(2)-X(6)*V(4)+X(2)*V(6)
  126. CHEB24(11) = CHEB12(11)+ALAM
  127. CHEB24(15) = CHEB12(11)-ALAM
  128. DO 30 I=1,3
  129. J = 8-I
  130. V(I) = FVAL(I)-FVAL(J)
  131. FVAL(I) = FVAL(I)+FVAL(J)
  132. 30 CONTINUE
  133. CHEB12(5) = V(1)+X(8)*V(3)
  134. CHEB12(9) = FVAL(1)-X(8)*FVAL(3)
  135. ALAM = X(4)*V(2)
  136. CHEB24(5) = CHEB12(5)+ALAM
  137. CHEB24(21) = CHEB12(5)-ALAM
  138. ALAM = X(8)*FVAL(2)-FVAL(4)
  139. CHEB24(9) = CHEB12(9)+ALAM
  140. CHEB24(17) = CHEB12(9)-ALAM
  141. CHEB12(1) = FVAL(1)+FVAL(3)
  142. ALAM = FVAL(2)+FVAL(4)
  143. CHEB24(1) = CHEB12(1)+ALAM
  144. CHEB24(25) = CHEB12(1)-ALAM
  145. CHEB12(13) = V(1)-V(3)
  146. CHEB24(13) = CHEB12(13)
  147. ALAM = 0.1E+01/0.6E+01
  148. DO 40 I=2,12
  149. CHEB12(I) = CHEB12(I)*ALAM
  150. 40 CONTINUE
  151. ALAM = 0.5E+00*ALAM
  152. CHEB12(1) = CHEB12(1)*ALAM
  153. CHEB12(13) = CHEB12(13)*ALAM
  154. DO 50 I=2,24
  155. CHEB24(I) = CHEB24(I)*ALAM
  156. 50 CONTINUE
  157. CHEB24(1) = 0.5E+00*ALAM*CHEB24(1)
  158. CHEB24(25) = 0.5E+00*ALAM*CHEB24(25)
  159. RETURN
  160. END