get_fourier_transform Subroutine

pure subroutine get_fourier_transform(kernel, kval, phi, dphi, grad)

   !> Prepared coefficients of the damped pair potential
   type(fourier_kernel), intent(in) :: kernel

   !> 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

   !> Whether the derivative is required
   logical, intent(in) :: grad

   integer :: ipole
   real(wp) :: arg, damp, acc, dacc, krad

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

   if (kval <= 0.0_wp) then
      phi = kernel%zero
      return
   end if

   krad = kval * kernel%radius
   acc = 0.0_wp
   dacc = 0.0_wp
   do ipole = 1, size(kernel%angle)
      arg = kernel%phase(ipole) + krad*kernel%cosine(ipole)
      damp = exp(-krad*kernel%sine(ipole))
      if (ipole <= kernel%npair) damp = 2.0_wp * damp
      acc = acc + sin(arg) * damp
      if (grad) dacc = dacc + cos(arg + kernel%angle(ipole)) * damp
   end do

   phi = kernel%prefactor * acc / kval
   if (grad) dphi = (kernel%prefactor * kernel%radius * dacc - phi) / kval

end subroutine get_fourier_transform