!> C ABI wrapper for the physical periodic zero mode.
module bem_periodic_zero_mode_c
  use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
  use, intrinsic :: iso_c_binding, only: c_associated, c_double, c_f_pointer, c_int, c_loc, c_ptr
  use bem_kinds, only: dp, i32, i64
  use bem_periodic_zero_mode_eval, only: eval_periodic_zero_mode, zero_mode_trace_minus, &
                                         zero_mode_trace_plus, zero_mode_trace_principal_value
  use bem_periodic_zero_mode_plan, only: build_periodic_zero_mode_height_plan, periodic_zero_mode_ok, &
                                         periodic_zero_mode_plan_type, periodic_zero_mode_state_type, &
                                         refresh_periodic_zero_mode_state
  implicit none
  private

  integer(c_int), parameter, public :: beach_zero_mode_ok = 0_c_int
  integer(c_int), parameter, public :: beach_zero_mode_invalid_handle = 1_c_int
  integer(c_int), parameter, public :: beach_zero_mode_invalid_argument = 2_c_int
  integer(c_int), parameter, public :: beach_zero_mode_not_ready = 3_c_int

  type :: periodic_zero_mode_handle
    type(periodic_zero_mode_plan_type) :: plan
    type(periodic_zero_mode_state_type) :: state
    logical :: built = .false.
    logical :: charged = .false.
  end type periodic_zero_mode_handle

  public :: beach_zero_mode_create
  public :: beach_zero_mode_destroy
  public :: beach_zero_mode_build
  public :: beach_zero_mode_update
  public :: beach_zero_mode_eval

contains

  integer(c_int) function beach_zero_mode_create(handle) &
    bind(C, name='beach_zero_mode_create') result(status)
    type(c_ptr), intent(out) :: handle
    type(periodic_zero_mode_handle), pointer :: zero

    allocate (zero)
    zero%built = .false.
    zero%charged = .false.
    handle = c_loc(zero)
    status = beach_zero_mode_ok
  end function beach_zero_mode_create

  integer(c_int) function beach_zero_mode_destroy(handle) &
    bind(C, name='beach_zero_mode_destroy') result(status)
    type(c_ptr), value :: handle
    type(periodic_zero_mode_handle), pointer :: zero

    status = get_zero_mode(handle, zero)
    if (status /= beach_zero_mode_ok) return

    deallocate (zero)
    status = beach_zero_mode_ok
  end function beach_zero_mode_destroy

  integer(c_int) function beach_zero_mode_build(handle, nsrc, source_heights_ptr, area_xy) &
    bind(C, name='beach_zero_mode_build') result(status)
    type(c_ptr), value :: handle
    integer(c_int), value :: nsrc
    type(c_ptr), value :: source_heights_ptr
    real(c_double), value :: area_xy
    type(periodic_zero_mode_handle), pointer :: zero
    real(c_double), pointer :: source_heights(:, :)
    integer(i32) :: plan_status
    character(len=128) :: message

    status = get_zero_mode(handle, zero)
    if (status /= beach_zero_mode_ok) return
    if (.not. count_is_addressable(nsrc, 3_i64, 2_i64) .or. nsrc == 0_c_int .or. &
        .not. c_associated(source_heights_ptr) .or. .not. ieee_is_finite(area_xy) .or. &
        area_xy <= 0.0_c_double) then
      status = beach_zero_mode_invalid_argument
      return
    end if

    call c_f_pointer(source_heights_ptr, source_heights, [3, int(nsrc)])
    if (any(.not. ieee_is_finite(source_heights))) then
      status = beach_zero_mode_invalid_argument
      return
    end if

    call build_periodic_zero_mode_height_plan( &
      real(source_heights, dp), real(area_xy, dp), zero%plan, plan_status, message &
      )
    if (plan_status /= periodic_zero_mode_ok) then
      zero%built = .false.
      zero%charged = .false.
      status = beach_zero_mode_invalid_argument
      return
    end if

    zero%built = .true.
    zero%charged = .false.
    status = beach_zero_mode_ok
  end function beach_zero_mode_build

  integer(c_int) function beach_zero_mode_update( &
    handle, nsrc, charge_ptr, e_bottom, z_gauge, phi_gauge &
    ) bind(C, name='beach_zero_mode_update') result(status)
    type(c_ptr), value :: handle
    integer(c_int), value :: nsrc
    type(c_ptr), value :: charge_ptr
    real(c_double), value :: e_bottom, z_gauge, phi_gauge
    type(periodic_zero_mode_handle), pointer :: zero
    real(c_double), pointer :: charge(:)

    status = get_zero_mode(handle, zero)
    if (status /= beach_zero_mode_ok) return
    if (.not. zero%built) then
      status = beach_zero_mode_not_ready
      return
    end if
    if (.not. count_is_addressable(nsrc, 1_i64, 0_i64) .or. &
        int(nsrc, i64) /= int(zero%plan%nelem, i64) .or. .not. c_associated(charge_ptr) .or. &
        .not. ieee_is_finite(e_bottom) .or. .not. ieee_is_finite(z_gauge) .or. &
        .not. ieee_is_finite(phi_gauge)) then
      status = beach_zero_mode_invalid_argument
      return
    end if

    call c_f_pointer(charge_ptr, charge, [int(nsrc)])
    if (any(.not. ieee_is_finite(charge))) then
      status = beach_zero_mode_invalid_argument
      return
    end if

    call refresh_periodic_zero_mode_state( &
      zero%plan, real(charge, dp), real(e_bottom, dp), real(z_gauge, dp), real(phi_gauge, dp), zero%state &
      )
    zero%charged = .true.
    status = beach_zero_mode_ok
  end function beach_zero_mode_update

  integer(c_int) function beach_zero_mode_eval(handle, ntarget, z_ptr, trace, phi_ptr, ez_ptr) &
    bind(C, name='beach_zero_mode_eval') result(status)
    type(c_ptr), value :: handle
    integer(c_int), value :: ntarget
    type(c_ptr), value :: z_ptr
    integer(c_int), value :: trace
    type(c_ptr), value :: phi_ptr, ez_ptr
    type(periodic_zero_mode_handle), pointer :: zero
    real(c_double), pointer :: z(:), phi(:), ez(:)
    real(dp) :: phi_value, field_value
    integer(i32) :: i

    status = require_charged_zero_mode(handle, zero)
    if (status /= beach_zero_mode_ok) return
    if (.not. count_is_addressable(ntarget, 3_i64, 0_i64) .or. &
        .not. c_associated(z_ptr) .or. .not. c_associated(phi_ptr) .or. .not. c_associated(ez_ptr) .or. &
        .not. trace_is_valid(trace)) then
      status = beach_zero_mode_invalid_argument
      return
    end if
    if (ntarget == 0_c_int) return

    call c_f_pointer(z_ptr, z, [int(ntarget)])
    call c_f_pointer(phi_ptr, phi, [int(ntarget)])
    call c_f_pointer(ez_ptr, ez, [int(ntarget)])
    if (any(.not. ieee_is_finite(z))) then
      status = beach_zero_mode_invalid_argument
      return
    end if

    do i = 1_i32, int(ntarget, i32)
      call eval_periodic_zero_mode( &
        zero%plan, zero%state, real(z(i), dp), int(trace, i32), phi_value, field_value &
        )
      phi(i) = real(phi_value, c_double)
      ez(i) = real(field_value, c_double)
    end do
    status = beach_zero_mode_ok
  end function beach_zero_mode_eval

  integer(c_int) function get_zero_mode(handle, zero) result(status)
    type(c_ptr), intent(in) :: handle
    type(periodic_zero_mode_handle), pointer, intent(out) :: zero

    if (.not. c_associated(handle)) then
      nullify (zero)
      status = beach_zero_mode_invalid_handle
      return
    end if

    call c_f_pointer(handle, zero)
    status = beach_zero_mode_ok
  end function get_zero_mode

  integer(c_int) function require_charged_zero_mode(handle, zero) result(status)
    type(c_ptr), intent(in) :: handle
    type(periodic_zero_mode_handle), pointer, intent(out) :: zero

    status = get_zero_mode(handle, zero)
    if (status /= beach_zero_mode_ok) return
    if (.not. zero%built .or. .not. zero%charged) status = beach_zero_mode_not_ready
  end function require_charged_zero_mode

  pure logical function count_is_addressable(count, values_per_item, shape_margin) result(addressable)
    integer(c_int), intent(in) :: count
    integer(i64), intent(in) :: values_per_item, shape_margin
    integer(i64) :: shape_limit

    shape_limit = min(int(huge(0), i64), int(huge(0_i32), i64))
    addressable = .false.
    if (count < 0_c_int .or. values_per_item <= 0_i64) return
    if (shape_margin < 0_i64 .or. shape_margin > shape_limit) return
    addressable = int(count, i64) <= (shape_limit - shape_margin)/values_per_item
  end function count_is_addressable

  pure logical function trace_is_valid(trace) result(valid)
    integer(c_int), intent(in) :: trace

    valid = trace == int(zero_mode_trace_minus, c_int) .or. &
            trace == int(zero_mode_trace_principal_value, c_int) .or. &
            trace == int(zero_mode_trace_plus, c_int)
  end function trace_is_valid

end module bem_periodic_zero_mode_c