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(i1)


45  phase=dx_shift*kvec


46  in(i)=out(i)*exp(cmplx(0d0,phase))


47  in(N+2i)=out(N/2+2i)*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

