subroutine panel_potential_field(geometry, charge, target, side, potential, field)
type(panel_geometry_type), intent(in) :: geometry
real(dp), intent(in) :: charge, target(3)
integer(i32), intent(in) :: side
real(dp), intent(out) :: potential, field(3)
real(dp) :: projection(3), height, potential_integral, field_integral(3), omega
real(dp) :: endpoint_a(3), endpoint_b(3), radius_a, radius_b, line_integral, distance, ratio
real(dp) :: relative_vertex(3, 3), numerator, denominator, scale, plane_tolerance
integer :: edge, next_edge
logical :: on_surface
if (side < panel_side_normal_minus .or. side > panel_side_normal_plus) then
error stop 'invalid panel evaluation side.'
end if
height = dot_product(target - geometry%vertex(:, 1), geometry%normal)
projection = target - height*geometry%normal
scale = maxval(geometry%edge_length)
plane_tolerance = 64.0_dp*epsilon(1.0_dp)*scale
on_surface = abs(height) <= plane_tolerance
if (on_surface) then
call panel_on_surface_integrals(geometry, projection, potential_integral, field_integral)
potential = k_coulomb*charge/geometry%area*potential_integral
field = k_coulomb*charge/geometry%area*field_integral
if (point_in_triangle(geometry, projection)) then
field = field + real(side, dp)*charge/(2.0_dp*geometry%area*eps0)*geometry%normal
end if
return
end if
omega = 0.0_dp
do edge = 1, 3
relative_vertex(:, edge) = geometry%vertex(:, edge) - target
end do
numerator = dot_product(relative_vertex(:, 1), &
cross_product(relative_vertex(:, 2), relative_vertex(:, 3)))
denominator = product_norms(relative_vertex) + &
dot_product(relative_vertex(:, 1), relative_vertex(:, 2))*norm2(relative_vertex(:, 3)) + &
dot_product(relative_vertex(:, 2), relative_vertex(:, 3))*norm2(relative_vertex(:, 1)) + &
dot_product(relative_vertex(:, 3), relative_vertex(:, 1))*norm2(relative_vertex(:, 2))
omega = -2.0_dp*atan2(numerator, denominator)
potential_integral = -height*omega
field_integral = geometry%normal*omega
do edge = 1, 3
next_edge = merge(edge + 1, 1, edge < 3)
endpoint_a = geometry%vertex(:, edge)
endpoint_b = geometry%vertex(:, next_edge)
radius_a = norm2(target - endpoint_a)
radius_b = norm2(target - endpoint_b)
ratio = geometry%edge_length(edge)/(radius_a + radius_b)
if (ratio >= 1.0_dp) error stop 'panel target lies on an edge or vertex.'
line_integral = log((1.0_dp + ratio)/(1.0_dp - ratio))
distance = dot_product(endpoint_a - projection, geometry%edge_outward(:, edge))
potential_integral = potential_integral + distance*line_integral
field_integral = field_integral + geometry%edge_outward(:, edge)*line_integral
end do
potential = k_coulomb*charge/geometry%area*potential_integral
field = k_coulomb*charge/geometry%area*field_integral
end subroutine panel_potential_field