| 1 | program test
|
|---|
| 2 | IMPLICIT NONE
|
|---|
| 3 | INCLUDE 'fftw3.f'
|
|---|
| 4 | INTEGER, PARAMETER :: N=200
|
|---|
| 5 | INTEGER(8) :: plan, plan2
|
|---|
| 6 | complex(8) :: in(N), out(N)
|
|---|
| 7 | REAL(8) :: x(N), pi, w(N), v(N), p(0:N/2), f(N), phase, length, dx_coarse, dx_fine, kvec, dx_shift, k_min
|
|---|
| 8 | iNTEGER :: i, mx_coarse, mx_fine
|
|---|
| 9 |
|
|---|
| 10 | length=1d0
|
|---|
| 11 | mx_coarse=N/2
|
|---|
| 12 | mx_fine=N
|
|---|
| 13 | dx_coarse=length/mx_coarse
|
|---|
| 14 | dx_fine=length/mx_fine
|
|---|
| 15 | k_min=2d0*acos(-1d0)/length
|
|---|
| 16 |
|
|---|
| 17 | pi=acos(-1d0)
|
|---|
| 18 | in=0
|
|---|
| 19 | out=0
|
|---|
| 20 | CALL dfftw_plan_dft_1d(plan, N/2, in, out, FFTW_FORWARD, FFTW_ESTIMATE)
|
|---|
| 21 |
|
|---|
| 22 |
|
|---|
| 23 | DO i=1,mx_coarse
|
|---|
| 24 | x(i)=(REAL(i)-.5d0) * dx_coarse
|
|---|
| 25 | w(i)=myfunc(x(i))
|
|---|
| 26 | END DO
|
|---|
| 27 | write(*,*) '# input_coarse'
|
|---|
| 28 | DO i=1,N/2
|
|---|
| 29 | write(*,*) x(i), w(i)
|
|---|
| 30 | END DO
|
|---|
| 31 | in(1:mx_coarse)=w(1:mx_coarse)
|
|---|
| 32 | CALL dfftw_execute(plan)
|
|---|
| 33 | CALL dfftw_destroy_plan(plan)
|
|---|
| 34 |
|
|---|
| 35 | CALL dfftw_plan_dft_1d(plan, N, in, out, FFTW_BACKWARD, FFTW_ESTIMATE)
|
|---|
| 36 |
|
|---|
| 37 | out=out/mx_coarse
|
|---|
| 38 |
|
|---|
| 39 | in=0
|
|---|
| 40 | in(1)=out(1)
|
|---|
| 41 | DO i=2,mx_coarse/2 !N/4
|
|---|
| 42 | !phase = k*dx_shift
|
|---|
| 43 | dx_shift=-dx_fine/2d0
|
|---|
| 44 | kvec=2d0*acos(-1d0)/length*real(i-1)
|
|---|
| 45 | phase=dx_shift*kvec
|
|---|
| 46 | in(i)=out(i)*exp(cmplx(0d0,phase))
|
|---|
| 47 | in(N+2-i)=out(N/2+2-i)*exp(cmplx(0d0,-phase))
|
|---|
| 48 | END DO
|
|---|
| 49 |
|
|---|
| 50 | CALL dfftw_execute(plan)
|
|---|
| 51 | DO i=1,mx_fine
|
|---|
| 52 | x(i)=(REAL(i)-.5d0)*dx_fine
|
|---|
| 53 | w(i)=myfunc(x(i))
|
|---|
| 54 | END DO
|
|---|
| 55 |
|
|---|
| 56 | write(*,*)
|
|---|
| 57 | write(*,*)
|
|---|
| 58 |
|
|---|
| 59 |
|
|---|
| 60 | write(*,*) '# target_fine'
|
|---|
| 61 | DO i=1,N
|
|---|
| 62 | write(*,*) x(i), w(i)
|
|---|
| 63 | END DO
|
|---|
| 64 |
|
|---|
| 65 | write(*,*)
|
|---|
| 66 | write(*,*)
|
|---|
| 67 |
|
|---|
| 68 |
|
|---|
| 69 | write(*,*) '# reconstruct_fine_real'
|
|---|
| 70 | DO i=1,N
|
|---|
| 71 | write(*,*) x(i), real(out(i))
|
|---|
| 72 | END DO
|
|---|
| 73 |
|
|---|
| 74 | write(*,*)
|
|---|
| 75 | write(*,*)
|
|---|
| 76 |
|
|---|
| 77 | write(*,*) '# reconstruct_fine_cmplx'
|
|---|
| 78 | DO i=1,N
|
|---|
| 79 | write(*,*) x(i), aimag(out(i))
|
|---|
| 80 | END DO
|
|---|
| 81 |
|
|---|
| 82 |
|
|---|
| 83 | CONTAINS
|
|---|
| 84 | REAL(8) function myfunc(x0)
|
|---|
| 85 | REAL(8) x0
|
|---|
| 86 | REAL(8), DIMENSION(5) :: wavemodes=(/1d0,2d0,3d0,5d0,9d0/)
|
|---|
| 87 | REAL(8), DIMENSION(5) :: amplitudes=(/5d0,2d0,1d0,2d0,4d0/)
|
|---|
| 88 | INTEGER :: i
|
|---|
| 89 | myfunc=0d0
|
|---|
| 90 | do i=1,5
|
|---|
| 91 | myfunc=myfunc+amplitudes(i)*cos(2d0*acos(-1d0)*wavemodes(i)*x0)
|
|---|
| 92 | end do
|
|---|
| 93 | END function myfunc
|
|---|
| 94 |
|
|---|
| 95 |
|
|---|
| 96 | end program test
|
|---|
