panel_singular_potential_oracle Subroutine

  subroutine panel_singular_potential_oracle(geometry, charge, target, order, potential)
    type(panel_geometry_type), intent(in) :: geometry
    real(dp), intent(in) :: charge, target(3)
    integer(i32), intent(in) :: order
    real(dp), intent(out) :: potential
    real(dp), allocatable :: node(:), weight(:)
    real(dp) :: edge1(3), edge2(3), direction(3), source(3), radius, jacobian, integral
    real(dp) :: u, v, sub_area2
    integer :: edge, next_edge, iu, iv

    call gauss_legendre_unit(order, node, weight)
    integral = 0.0_dp
    do edge = 1, 3
      next_edge = merge(edge + 1, 1, edge < 3)
      edge1 = geometry%vertex(:, edge) - target
      edge2 = geometry%vertex(:, next_edge) - target
      sub_area2 = sqrt(sum(cross_product(edge1, edge2)**2))
      do iu = 1, order
        u = node(iu)
        do iv = 1, order
          v = node(iv)
          direction = (1.0_dp - v)*edge1 + v*edge2
          source = target + u*direction
          radius = sqrt(sum((target - source)**2))
          jacobian = sub_area2*u*weight(iu)*weight(iv)
          integral = integral + jacobian/radius
        end do
      end do
    end do
    potential = k_coulomb*charge/geometry%area*integral
  end subroutine panel_singular_potential_oracle