program test
   IMPLICIT NONE
   INCLUDE 'fftw3.f'
   INTEGER, PARAMETER :: N=200
   INTEGER(8) :: plan, plan2
   complex(8) :: in(N), out(N)
   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
   iNTEGER :: i, mx_coarse, mx_fine

   length=1d0
   mx_coarse=N/2
   mx_fine=N
   dx_coarse=length/mx_coarse
   dx_fine=length/mx_fine
   k_min=2d0*acos(-1d0)/length

   pi=acos(-1d0)
   in=0
   out=0
   CALL dfftw_plan_dft_1d(plan, N/2, in, out, FFTW_FORWARD, FFTW_ESTIMATE)  


   DO i=1,mx_coarse
      x(i)=(REAL(i)-.5d0) * dx_coarse
      w(i)=myfunc(x(i))
   END DO
   write(*,*) '# input_coarse'
   DO i=1,N/2
      write(*,*) x(i), w(i)
   END DO
   in(1:mx_coarse)=w(1:mx_coarse)
   CALL dfftw_execute(plan)
   CALL dfftw_destroy_plan(plan)

   CALL dfftw_plan_dft_1d(plan, N, in, out, FFTW_BACKWARD, FFTW_ESTIMATE)  

   out=out/mx_coarse

   in=0
   in(1)=out(1)
   DO i=2,mx_coarse/2 !N/4
      !phase = k*dx_shift
      dx_shift=-dx_fine/2d0
      kvec=2d0*acos(-1d0)/length*real(i-1)
      phase=dx_shift*kvec
      in(i)=out(i)*exp(cmplx(0d0,phase))
      in(N+2-i)=out(N/2+2-i)*exp(cmplx(0d0,-phase))
   END DO

   CALL dfftw_execute(plan)
   DO i=1,mx_fine
      x(i)=(REAL(i)-.5d0)*dx_fine
      w(i)=myfunc(x(i))
   END DO

   write(*,*) 
   write(*,*) 


   write(*,*) '# target_fine'
   DO i=1,N
      write(*,*) x(i), w(i)
   END DO

   write(*,*) 
   write(*,*) 


   write(*,*) '# reconstruct_fine_real'
   DO i=1,N
      write(*,*) x(i), real(out(i))
   END DO

   write(*,*) 
   write(*,*) 

   write(*,*) '# reconstruct_fine_cmplx'
   DO i=1,N
      write(*,*) x(i), aimag(out(i))
   END DO


CONTAINS
   REAL(8) function myfunc(x0)
      REAL(8) x0
      REAL(8), DIMENSION(5) :: wavemodes=(/1d0,2d0,3d0,5d0,9d0/)
      REAL(8), DIMENSION(5) :: amplitudes=(/5d0,2d0,1d0,2d0,4d0/)
      INTEGER :: i
      myfunc=0d0
      do i=1,5
         myfunc=myfunc+amplitudes(i)*cos(2d0*acos(-1d0)*wavemodes(i)*x0)
      end do
   END function myfunc


end program test