integer(c_int) function beach_kernel_force_on_charges( &
handle, ntarget, target_pos_ptr, target_q_ptr, origin_ptr, force_ptr, torque_ptr &
) bind(C, name='beach_kernel_force_on_charges') result(status)
type(c_ptr), value :: handle
integer(c_int), value :: ntarget
type(c_ptr), value :: target_pos_ptr
type(c_ptr), value :: target_q_ptr
type(c_ptr), value :: origin_ptr
type(c_ptr), value :: force_ptr
type(c_ptr), value :: torque_ptr
type(field_kernel_handle), pointer :: kernel
real(c_double), pointer :: target_pos(:, :)
real(c_double), pointer :: target_q(:)
real(c_double), pointer :: origin(:)
real(c_double), pointer :: force(:)
real(c_double), pointer :: torque(:)
real(dp), allocatable :: e(:, :)
real(dp) :: f_i(3), r_rel(3)
integer(i32) :: i
status = require_charged_kernel(handle, kernel)
if (status /= beach_kernel_ok) return
if (ntarget < 0_c_int .or. .not. c_associated(target_pos_ptr) .or. .not. c_associated(target_q_ptr) .or. &
.not. c_associated(origin_ptr) .or. .not. c_associated(force_ptr) .or. .not. c_associated(torque_ptr)) then
status = beach_kernel_invalid_argument
return
end if
call c_f_pointer(force_ptr, force, [3])
call c_f_pointer(torque_ptr, torque, [3])
force = 0.0_c_double
torque = 0.0_c_double
if (ntarget == 0_c_int) return
call c_f_pointer(target_pos_ptr, target_pos, [3, int(ntarget)])
call c_f_pointer(target_q_ptr, target_q, [int(ntarget)])
call c_f_pointer(origin_ptr, origin, [3])
allocate (e(3, int(ntarget)))
call eval_points(kernel%plan, kernel%state, target_pos, e)
e = k_coulomb*e
do i = 1_i32, int(ntarget, i32)
f_i = real(target_q(i), dp)*e(:, i)
r_rel = target_pos(:, i) - real(origin, dp)
force = force + real(f_i, c_double)
torque(1) = torque(1) + real(r_rel(2)*f_i(3) - r_rel(3)*f_i(2), c_double)
torque(2) = torque(2) + real(r_rel(3)*f_i(1) - r_rel(1)*f_i(3), c_double)
torque(3) = torque(3) + real(r_rel(1)*f_i(2) - r_rel(2)*f_i(1), c_double)
end do
deallocate (e)
status = beach_kernel_ok
end function beach_kernel_force_on_charges