compute_laplace_derivatives Subroutine

  subroutine compute_laplace_derivatives(plan, r, deriv)
    type(fmm_plan_type), intent(in) :: plan
    real(dp), intent(in) :: r(3)
    real(dp), intent(out) :: deriv(:)
    real(dp) :: a, coeff, soft2
    real(dp) :: q(plan%nderiv), power(plan%nderiv), accum(plan%nderiv), tmp(plan%nderiv)
    integer(i32) :: idx0, idx, n

    if (size(deriv) /= plan%nderiv) error stop 'Derivative buffer size mismatch.'
    soft2 = plan%options%softening*plan%options%softening
    a = dot_product(r, r) + soft2
    if (a <= tiny(1.0d0)) error stop 'Cannot expand Laplace kernel at r=0.'

    q = 0.0d0
    power = 0.0d0
    accum = 0.0d0
    idx0 = plan%deriv_map(0_i32, 0_i32, 0_i32)
    power(idx0) = 1.0d0
    accum(idx0) = 1.0d0
    if (plan%options%order >= 1_i32) then
      q(plan%deriv_map(1_i32, 0_i32, 0_i32)) = 2.0d0*r(1)/a
      q(plan%deriv_map(0_i32, 1_i32, 0_i32)) = 2.0d0*r(2)/a
      q(plan%deriv_map(0_i32, 0_i32, 1_i32)) = 2.0d0*r(3)/a
    end if
    if (2_i32 <= 2_i32*plan%options%order) then
      q(plan%deriv_map(2_i32, 0_i32, 0_i32)) = q(plan%deriv_map(2_i32, 0_i32, 0_i32)) + 1.0d0/a
      q(plan%deriv_map(0_i32, 2_i32, 0_i32)) = q(plan%deriv_map(0_i32, 2_i32, 0_i32)) + 1.0d0/a
      q(plan%deriv_map(0_i32, 0_i32, 2_i32)) = q(plan%deriv_map(0_i32, 0_i32, 2_i32)) + 1.0d0/a
    end if

    coeff = 1.0d0
    do n = 1_i32, 2_i32*plan%options%order
      call multiply_polynomials(plan, power, q, tmp)
      power = tmp
      coeff = coeff*(-(2.0d0*real(n, dp) - 1.0d0)/(2.0d0*real(n, dp)))
      accum = accum + coeff*power
    end do
    accum = accum/sqrt(a)
    do idx = 1_i32, plan%nderiv
      deriv(idx) = accum(idx)*plan%deriv_factorial(idx)
    end do
  end subroutine compute_laplace_derivatives