panel_oracle_potential_field Subroutine

  subroutine panel_oracle_potential_field(geometry, charge, target, order, potential, field)
    type(panel_geometry_type), intent(in) :: geometry
    real(dp), intent(in) :: charge, target(3)
    integer(i32), intent(in) :: order
    real(dp), intent(out) :: potential, field(3)
    real(dp), allocatable :: node(:), weight(:)
    real(dp) :: source(3), displacement(3), radius2, jacobian, potential_integral, field_integral(3)
    real(dp) :: u, v, edge1(3), edge2(3), direction(3)
    integer :: iu, iv

    call gauss_legendre_unit(order, node, weight)
    edge1 = geometry%vertex(:, 2) - geometry%vertex(:, 1)
    edge2 = geometry%vertex(:, 3) - geometry%vertex(:, 1)
    potential_integral = 0.0_dp
    field_integral = 0.0_dp
    do iu = 1, order
      u = node(iu)
      do iv = 1, order
        v = node(iv)
        direction = (1.0_dp - v)*edge1 + v*edge2
        source = geometry%vertex(:, 1) + u*direction
        displacement = target - source
        radius2 = sum(displacement*displacement)
        jacobian = 2.0_dp*geometry%area*u*weight(iu)*weight(iv)
        potential_integral = potential_integral + jacobian/sqrt(radius2)
        field_integral = field_integral + jacobian*displacement/(radius2*sqrt(radius2))
      end do
    end do
    potential = k_coulomb*charge/geometry%area*potential_integral
    field = k_coulomb*charge/geometry%area*field_integral
  end subroutine panel_oracle_potential_field