| 1 | ! Transformee de Fourier complexe-complexe 1d |
|---|
| 2 | |
|---|
| 3 | program tjmccfft |
|---|
| 4 | |
|---|
| 5 | implicit none |
|---|
| 6 | |
|---|
| 7 | integer, parameter :: n = 36 |
|---|
| 8 | complex, dimension(0:n-1) :: x, xx |
|---|
| 9 | real, dimension(0:n-1,2) :: rx |
|---|
| 10 | equivalence( x, rx ) |
|---|
| 11 | |
|---|
| 12 | ! Pour stocker les cosinus et les sinus |
|---|
| 13 | ! En fait, pour les routines jm, 100+2*n suffisent |
|---|
| 14 | integer, parameter :: ntable = 100+8*n |
|---|
| 15 | real, dimension(0:ntable-1) :: table |
|---|
| 16 | |
|---|
| 17 | ! En fait, pour les routines jm, 4*n suffisent |
|---|
| 18 | integer, parameter :: nwork = 8*n |
|---|
| 19 | real, dimension(nwork) :: work |
|---|
| 20 | |
|---|
| 21 | integer :: isign |
|---|
| 22 | real :: scale |
|---|
| 23 | integer :: isys |
|---|
| 24 | integer :: i, j |
|---|
| 25 | real :: twopi |
|---|
| 26 | complex :: s |
|---|
| 27 | |
|---|
| 28 | ! On prepare le tableau d'entree |
|---|
| 29 | call random_number( rx ) |
|---|
| 30 | xx = x |
|---|
| 31 | |
|---|
| 32 | scale = 1. |
|---|
| 33 | isys = 0 |
|---|
| 34 | |
|---|
| 35 | isign = 0 |
|---|
| 36 | call ccfft(isign,n,scale,x,x,table,work,isys) |
|---|
| 37 | isign = 1 |
|---|
| 38 | print *,'jmccfft ',n,0,isign,scale |
|---|
| 39 | call ccfft(isign,n,scale,x,x,table,work,isys) |
|---|
| 40 | |
|---|
| 41 | ! On imprime le tableau de sortie |
|---|
| 42 | open(10,file='temp1',status='unknown',form='formatted') |
|---|
| 43 | write(10,'(2e25.12)') x |
|---|
| 44 | |
|---|
| 45 | ! Ce qu'il faut trouver |
|---|
| 46 | open(11,file='temp2',status='unknown',form='formatted') |
|---|
| 47 | twopi = 2 * acos(real(-1)) |
|---|
| 48 | x = xx |
|---|
| 49 | do i = 0,n-1 |
|---|
| 50 | s = 0 |
|---|
| 51 | do j = 0,n-1 |
|---|
| 52 | s = s + cmplx(cos(twopi*i*j/real(n)),isign*sin(twopi*i*j/real(n)))*x(j) |
|---|
| 53 | end do |
|---|
| 54 | write(11,'(2e25.12)') s*scale |
|---|
| 55 | end do |
|---|
| 56 | |
|---|
| 57 | end program tjmccfft |
|---|
