Program poisson2d
implicit none
integer ( kind = 4 ), parameter :: Nx = 128
integer ( kind = 4 ), parameter :: Ny = 128
integer ( kind = 4 ), parameter :: Nh = ( Nx / 2 ) + 1
real ( kind = 8 ), parameter :: pi=3.14159265358979323846d0
include "fftw3.f"
integer ( kind = 4 ) i,j
real ( kind = 8 ) Lx,Ly,dx,dy,kx,ky
real ( kind = 8 ) x(Nx),y(Ny),phi(Nx,Ny),rho(Nx,Ny),rho_dum(Nx,Ny)
complex ( kind = 8 ) rhok(Nh,Ny),rhok_dum(Nh,Ny),phik(Nh,Ny),phik_dum(Nh,Ny)
integer ( kind = 8 ) plan_forward,plan_backward
open(unit=10,file='poisson.dat',status='unknown')
Lx = 2.0*pi
Ly = 2.0*pi
dx = Lx/float(Nx)
dy = Ly/float(Ny)
do i = 1, Nx
x(i)=0.0d0+real(i-1)*dx
do j = 1, Ny
y(j)=0.0+real(j-1)*dy
rho(i,j) = 2.0d0*sin(x(i))*cos(y(j))
rho_dum(i,j) = rho(i,j)
end do
end do
call dfftw_plan_dft_r2c_2d_ (plan_forward, Nx, Ny, rho_dum, rhok, FFTW_ESTIMATE)
call dfftw_execute_ (plan_forward)
do i = 1,Nx/2+1
do j = 1,Ny/2
kx = 2.0d0*pi*float(i-1)/Lx
ky = 2.0d0*pi*float(j-1)/Ly
if (i == 1 .and. j == 1) then
phik(i,j) = (0.0d0,0.0d0)
else
phik(i,j) = rhok(i,j)/( kx*kx + ky*ky )
endif
phik_dum(i,j) = phik(i,j)
enddo
do j = Ny/2+1,Ny
kx = 2.0d0*pi*float(i-1)/Lx
ky = 2.0d0*pi*float((j-1)-Ny)/Ly
phik(i,j) = rhok(i,j)/( kx*kx + ky*ky )
phik_dum(i,j) = phik(i,j)
enddo
enddo
call dfftw_plan_dft_c2r_2d_ (plan_backward, Nx, Ny, phik_dum, phi, FFTW_ESTIMATE)
call dfftw_execute_ (plan_backward)
call dfftw_destroy_plan_ (plan_forward)
call dfftw_destroy_plan_ (plan_backward)
do i = 1, Nx
do j = 1, Ny
phi(i,j) = phi(i,j)/(float(Nx)*float(Ny))
write(10,*) x(i),y(j),rho(i,j),phi(i,j)
end do
end do
end program poisson2d