FFTW based One Dimensional Serial Poisson Solver

! Rupak, IPR, 02.06.2017

! Discussion.

! This is an One Dimensional FFTW based Serial Poisson Solver.

!

program poisson1d

    integer ( kind = 4 ), parameter :: N = 1024

    include "fftw3.f"

    double precision:: x, rho, phi

    dimension x(N), rho(N), phi(N)

    double precision:: rho_k_mod

    dimension rho_k_mod(N/2+1)

    double complex:: rho_k

    dimension rho_k(N/2+1)

    double precision:: L,dx

    double precision, parameter:: pi=3.14159265358979323846

    integer:: i

    integer*8:: plan

    open(unit=10,file='FFT.dat',status='unknown')

    open(unit=12,file='DATA.dat',status='unknown')   

    L = 2.0d0*pi

    dx = L/dfloat(N)

    do i=1,N

      x(i) = 0.0d0 + dfloat(i-1)*dx

      rho(i)= - 1.0d0 + dsin(x(i))

    enddo

    call dfftw_plan_dft_r2c_1d(plan,N,rho,rho_k,FFTW_ESTIMATE)

    call dfftw_execute_dft_r2c(plan,rho,rho_k)

    call dfftw_destroy_plan(plan)

    do i=1,N/2+1

      rho_k_mod(i) = abs(rho_k(i))/dfloat(N/2+1)

      write(10,*) i-1, rho_k_mod(i)

    enddo

    do i= 2,N/2+1

      k = 2.0d0*pi*dfloat(i-1)/L

      rho_k(i) = - rho_k(i)/(k*k)

      rho_k(i) = rho_k(i)

    enddo

    call dfftw_plan_dft_c2r_1d(plan,N,rho_k,phi,FFTW_ESTIMATE)

    call dfftw_execute_dft_c2r(plan,rho_k,phi)

    call dfftw_destroy_plan(plan)

    do i=1,N

      phi(i) = phi(i)/dfloat(N)

      write(12,*) x(i), rho(i), phi(i)

    enddo

    close(10)

    close(12)

end program poisson1d