get_fourier_transform Subroutine

pure subroutine get_fourier_transform(term, kval, phi, dphi)

   !> Term of the damped pair potential
   type(fourier_term), intent(in) :: term

   !> Wave number
   real(wp), intent(in) :: kval

   !> Fourier transform of the term
   real(wp), intent(out) :: phi

   !> Derivative of the Fourier transform w.r.t. the wave number
   real(wp), intent(out) :: dphi

   integer :: ipole, npole, nexp
   real(wp) :: pre, theta, ct, st, arg, damp, acc, dacc, krad

   phi = 0.0_wp
   dphi = 0.0_wp
   if (abs(term%prefactor) <= 0.0_wp) return

   nexp = term%m - term%alpha + 2
   pre = 4.0_wp*pi*pi / term%alpha * term%prefactor * term%radius**nexp

   if (kval <= 0.0_wp) then
      phi = pre * term%radius / sin(pi*(term%m + 3)/real(term%alpha, wp))
      return
   end if

   npole = term%alpha / 2
   krad = kval * term%radius
   acc = 0.0_wp
   dacc = 0.0_wp
   do ipole = 0, npole - 1
      theta = pi * (2*ipole + 1) / real(term%alpha, wp)
      ct = cos(theta)
      st = sin(theta)
      arg = 0.5_wp*pi + nexp*theta + krad*ct
      damp = exp(-krad*st)
      acc = acc + sin(arg) * damp
      dacc = dacc + cos(arg + theta) * damp
   end do

   phi = pre * acc / kval
   dphi = (pre * term%radius * dacc - phi) / kval

end subroutine get_fourier_transform