! 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