Télécharger rfft1b.eso

Retour à la liste

Numérotation des lignes :

rfft1b
  1. C RFFT1B SOURCE BP208322 18/10/08 21:15:17 9952
  2. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3. C
  4. C FFTPACK 5.1
  5. C
  6. C Authors: Paul N. Swarztrauber and Richard A. Valent
  7. C
  8. c FFTPACK 5.1 routine RFFT1F computes the one-dimensional Fourier
  9. c transform of a periodic sequence within a real array. This
  10. c is referred to as the BACKWARD transform or Fourier analysis,
  11. c transforming the sequence from SPECTRAL to PHYSICAL space.
  12. c
  13. c
  14. c This transform is normalized since a call to RFFT1F followed
  15. c by a call to RFFT1B (or vice-versa) reproduces the original
  16. c array subject to algorithmic constraints, roundoff error, etc.
  17. c
  18. c Input Arguments
  19. c
  20. c N Integer length of the sequence to be transformed. The
  21. c transform is most efficient when N is a product of
  22. c small primes.
  23. c
  24. c INC Integer increment between the locations, in array R,
  25. c of two consecutive elements within the sequence.
  26. c
  27. c R Real array of length LENR containing the sequence to be
  28. c transformed.
  29. c
  30. c LENR Integer dimension of R array. LENR must be at least
  31. c INC*(N-1) + 1.
  32. c
  33. c WSAVE Real work array of length LENSAV. WSAVE's contents must
  34. c be initialized with a call to subroutine RFFT1I before the
  35. c first call to routine RFFT1F or RFFT1B for a given transform
  36. c length N.
  37. c
  38. c
  39. c LENSAV Integer dimension of WSAVE array. LENSAV must be at least
  40. c N + INT(LOG (REAL(N))/LOG(2.)) +4.
  41. c
  42. c
  43. c WORK Real work array of dimension LENWRK.
  44. c
  45. c LENWRK Integer dimension of WORK array. LENWRK must be at N.
  46. c
  47. c
  48. c Output Arguments
  49. c
  50. c R Real output array R. For purposes of exposition,
  51. c assume R's range of indices is given by
  52. c R(0:(N-1)*INC).
  53. c
  54. c Then
  55. c
  56. c N-1
  57. c R(0) = SUM R(N1*INC)/N
  58. c N1=0
  59. c
  60. c If N is even, set NH=N/2-1; if N is odd set NH=(N-1)/2;
  61. c then for J=1,...,NH
  62. c
  63. c R((2*J-1)*INC) =
  64. c
  65. c N-1
  66. c 2.*SUM (R(N1*INC)*COS(J*N1*2*PI/N)/N
  67. c N1=0
  68. c
  69. c and
  70. c
  71. c R(2*J*INC) =
  72. c
  73. c N-1
  74. c 2.*SUM (R(N1*INC)*SIN(J*N1*2*PI/N)/N
  75. c N1=0
  76. c
  77. c Also if N is even then
  78. c
  79. c R((N-1)*INC) =
  80. c
  81. c N-1
  82. c SUM (-1)**N1*R(N1*INC)/N
  83. c N1=0
  84. c
  85. c
  86. c IER Integer error return
  87. c = 0 successful exit
  88. c = 1 input parameter LENR not big enough
  89. c = 2 input parameter LENSAV not big enough
  90. c = 3 input parameter LENWRK not big enough
  91. c
  92. c
  93. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  94.  
  95. SUBROUTINE RFFT1B ( N, INC, R, LENR, WSAVE, LENSAV,
  96. 1 WORK, LENWRK, IER)
  97. IMPLICIT INTEGER(I-N)
  98. IMPLICIT REAL*8(A-H,O-Z)
  99. INTEGER N, INC, LENR, LENSAV, LENWRK, IER
  100. REAL*8 R(LENR), WSAVE(LENSAV) ,WORK(LENWRK)
  101. C
  102. IER = 0
  103. C
  104. IF (LENR .LT. INC*(N-1) + 1) THEN
  105. IER = 1
  106. c CALL XERFFT ('RFFT1B ', 6)
  107. call erreur(5)
  108. ELSEIF (LENSAV .LT. N + INT(LOG(REAL(N))/LOG(2.)) +4) THEN
  109. IER = 2
  110. c CALL XERFFT ('RFFT1B ', 8)
  111. call erreur(5)
  112. ELSEIF (LENWRK .LT. N) THEN
  113. IER = 3
  114. c CALL XERFFT ('RFFT1B ', 10)
  115. call erreur(5)
  116. ENDIF
  117. C
  118. IF (N .EQ. 1) RETURN
  119. C
  120. CALL RFFTB1 (N,INC,R,WORK,WSAVE,WSAVE(N+1))
  121. RETURN
  122. END
  123.  
  124.  
  125.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales