!> 同一時刻の粒子状態から、空間電場を中点評価した1ステップ候補を構築する。
module bem_particle_stepper
  use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
  use bem_kinds, only: dp, i32, i64
  use bem_types, only: mesh_type, sim_config, hit_info, bc_open, bc_reflect, bc_periodic, bc_redistributed_reflect
  use bem_electrostatic_snapshot, only: electrostatic_snapshot_type
  use bem_pusher, only: boris_push
  use bem_collision, only: collision_query_ok, find_first_hit
  use bem_boundary, only: boundary_event_type, boundary_event_ok, boundary_event_invalid_geometry, find_first_boundary_event, &
                          apply_escape_reflect_periodic_event
  use bem_external_boundary_contract, only: &
    external_boundary_contract_type, external_open_escape, external_open_potential_barrier
  implicit none
  private

  integer(i32), parameter, public :: particle_step_ok = collision_query_ok
  integer(i32), parameter, public :: particle_step_invalid_boundary = 1001_i32
  integer(i32), parameter, public :: particle_step_multiple_box_events = 1002_i32
  integer(i32), parameter, public :: particle_step_ambiguous_open_corner = 1003_i32
  integer(i32), parameter, public :: particle_step_unsupported_barrier_corner = particle_step_ambiguous_open_corner

  type, public :: particle_step_result
    real(dp) :: x(3) = 0.0_dp
    real(dp) :: v(3) = 0.0_dp
    logical :: absorbed = .false.
    logical :: escaped_boundary = .false.
    integer(i32) :: elem_idx = -1_i32
    integer(i32) :: status = particle_step_ok
    integer(i32) :: field_eval_count = 0_i32
    integer(i32) :: collision_query_count = 0_i32
    integer(i32) :: z_high_outward_event_count = 0_i32
    integer(i32) :: outer_barrier_return_count = 0_i32
    integer(i32) :: outer_barrier_escape_count = 0_i32
    real(dp) :: z_high_outward_normal_kinetic_energy_j_sum = 0.0_dp
  end type particle_step_result

  public :: build_particle_step_candidate
  public :: advance_particle_step
  public :: advance_particle_step_upper_panel_fourier
  public :: resolve_particle_boundary_candidate

contains

  !> 予測中点の電場と一様磁場を使い、次時刻の位置・速度候補を返す。
  subroutine build_particle_step_candidate( &
    mesh, sim, snapshot, bfield, x0, v0, q, m, dt, x1, v1, sampled_electric_field, &
    use_upper_panel_fourier, field_available &
    )
    type(mesh_type), intent(in) :: mesh
    type(sim_config), intent(in) :: sim
    type(electrostatic_snapshot_type), intent(inout) :: snapshot
    real(dp), intent(in) :: bfield(3), x0(3), v0(3), q, m, dt
    real(dp), intent(out) :: x1(3), v1(3)
    real(dp), intent(out), optional :: sampled_electric_field(3)
    logical, intent(in), optional :: use_upper_panel_fourier
    logical, intent(out), optional :: field_available
    real(dp) :: x_mid(3), e_mid(3)
    logical :: available, use_upper

    x_mid = x0 + 0.5d0*v0*dt
    call project_field_sample_to_box(sim, x_mid)
    use_upper = .false.
    if (present(use_upper_panel_fourier)) use_upper = use_upper_panel_fourier
    available = .true.
    if (use_upper) then
      call snapshot%eval_upper_panel_fourier_e(mesh, x_mid, e_mid, available)
    else
      call snapshot%eval_local_e(mesh, x_mid, e_mid)
    end if
    if (present(field_available)) field_available = available
    if (.not. available) then
      x1 = x0
      v1 = v0
      if (present(sampled_electric_field)) sampled_electric_field = 0.0_dp
      return
    end if
    call boris_push(x0, v0, q, m, dt, e_mid, bfield, x1, v1)
    if (present(sampled_electric_field)) sampled_electric_field = e_mid
  end subroutine build_particle_step_candidate

  !> 境界を越えるcandidateでも、場評価点はsolverのprimitive target box内に保つ。
  pure subroutine project_field_sample_to_box(sim, position)
    type(sim_config), intent(in) :: sim
    real(dp), intent(inout) :: position(3)
    integer(i32) :: axis
    real(dp) :: span

    if (.not. sim%use_box) return
    do axis = 1_i32, 3_i32
      if (.not. ieee_is_finite(sim%box_min(axis)) .or. .not. ieee_is_finite(sim%box_max(axis))) cycle
      span = sim%box_max(axis) - sim%box_min(axis)
      if (.not. ieee_is_finite(span) .or. span <= 0.0_dp) cycle
      if (sim%bc_low(axis) == bc_periodic .and. sim%bc_high(axis) == bc_periodic) then
        position(axis) = sim%box_min(axis) + modulo(position(axis) - sim%box_min(axis), span)
      else
        position(axis) = min(max(position(axis), sim%box_min(axis)), sim%box_max(axis))
      end if
    end do
  end subroutine project_field_sample_to_box

  !> 一つの粒子stepについてmesh/boxの最早eventを順序付ける。
  subroutine advance_particle_step( &
    mesh, sim, snapshot, bfield, x0, v0, q, m, dt, result, boundary_contract, boundary_rng_counter &
    )
    type(mesh_type), intent(in) :: mesh
    type(sim_config), intent(in) :: sim
    type(electrostatic_snapshot_type), intent(inout) :: snapshot
    real(dp), intent(in) :: bfield(3), x0(3), v0(3), q, m, dt
    type(particle_step_result), intent(out) :: result
    type(external_boundary_contract_type), intent(in), optional :: boundary_contract
    integer(i64), intent(in), optional :: boundary_rng_counter(4)
    logical :: field_available

    field_available = .true.
    call advance_particle_step_impl( &
      mesh, sim, snapshot, bfield, x0, v0, q, m, dt, result, boundary_contract, boundary_rng_counter, 0_i32, &
      .false., field_available &
      )
  end subroutine advance_particle_step

  !> 上部真空域のfactorized P0-panel Fourier場だけで、失敗した1 stepを最初から再試行する。
  subroutine advance_particle_step_upper_panel_fourier( &
    mesh, sim, snapshot, bfield, x0, v0, q, m, dt, result, field_available, &
    boundary_contract, boundary_rng_counter &
    )
    type(mesh_type), intent(in) :: mesh
    type(sim_config), intent(in) :: sim
    type(electrostatic_snapshot_type), intent(inout) :: snapshot
    real(dp), intent(in) :: bfield(3), x0(3), v0(3), q, m, dt
    type(particle_step_result), intent(out) :: result
    logical, intent(out) :: field_available
    type(external_boundary_contract_type), intent(in), optional :: boundary_contract
    integer(i64), intent(in), optional :: boundary_rng_counter(4)

    field_available = .true.
    call advance_particle_step_impl( &
      mesh, sim, snapshot, bfield, x0, v0, q, m, dt, result, boundary_contract, boundary_rng_counter, 0_i32, &
      .true., field_available &
      )
  end subroutine advance_particle_step_upper_panel_fourier

  !> periodic event の適応分割深さを内部だけで引き回す1 step実装。
  recursive subroutine advance_particle_step_impl( &
    mesh, sim, snapshot, bfield, x0, v0, q, m, dt, result, boundary_contract, boundary_rng_counter, adaptive_depth, &
    use_upper_panel_fourier, field_available &
    )
    type(mesh_type), intent(in) :: mesh
    type(sim_config), intent(in) :: sim
    type(electrostatic_snapshot_type), intent(inout) :: snapshot
    real(dp), intent(in) :: bfield(3), x0(3), v0(3), q, m, dt
    type(particle_step_result), intent(out) :: result
    type(external_boundary_contract_type), intent(in), optional :: boundary_contract
    integer(i64), intent(in), optional :: boundary_rng_counter(4)
    integer(i32), intent(in) :: adaptive_depth
    logical, intent(in) :: use_upper_panel_fourier
    logical, intent(inout) :: field_available

    type(hit_info) :: hit
    real(dp) :: x_candidate(3), v_candidate(3), candidate_electric_field(3)
    integer(i32) :: query_status
    integer(i64) :: active_boundary_rng_counter(4)
    logical :: has_boundary_rng_counter

    result = particle_step_result()
    result%x = x0
    result%v = v0
    active_boundary_rng_counter = 0_i64
    has_boundary_rng_counter = present(boundary_rng_counter)
    if (has_boundary_rng_counter) active_boundary_rng_counter = boundary_rng_counter
    if (.not. valid_particle_step_input(x0, v0, bfield, q, m, dt)) then
      result%status = particle_step_invalid_boundary
      return
    end if

    call build_particle_step_candidate( &
      mesh, sim, snapshot, bfield, x0, v0, q, m, dt, x_candidate, v_candidate, candidate_electric_field, &
      use_upper_panel_fourier, field_available &
      )
    result%field_eval_count = 1_i32
    if (.not. field_available) then
      result%status = particle_step_invalid_boundary
      return
    end if
    if (.not. all(ieee_is_finite(x_candidate)) .or. .not. all(ieee_is_finite(v_candidate))) then
      result%status = particle_step_invalid_boundary
      return
    end if

    call find_first_hit(mesh, x0, x_candidate, hit, sim=sim, status=query_status)
    result%collision_query_count = 1_i32
    if (query_status /= collision_query_ok) then
      if (sim%use_box .and. .not. point_strictly_inside_box(sim, x_candidate)) then
        call advance_particle_boundary_crossing( &
          mesh, sim, snapshot, bfield, x0, v0, q, m, dt, x_candidate, v_candidate, candidate_electric_field, result=result, &
          boundary_contract=boundary_contract, boundary_rng_counter=active_boundary_rng_counter, &
          has_boundary_rng_counter=has_boundary_rng_counter, adaptive_depth=adaptive_depth, &
          use_upper_panel_fourier=use_upper_panel_fourier, field_available=field_available &
          )
        return
      end if
      result%status = query_status
      return
    end if
    if (.not. sim%use_box .or. point_strictly_inside_box(sim, x_candidate)) then
      if (hit%has_hit) then
        call accept_particle_hit(v0, v_candidate, hit, result)
        return
      end if
      result%x = x_candidate
      result%v = v_candidate
      return
    end if

    call advance_particle_boundary_crossing( &
      mesh, sim, snapshot, bfield, x0, v0, q, m, dt, x_candidate, v_candidate, candidate_electric_field, hit, result, &
      boundary_contract=boundary_contract, boundary_rng_counter=active_boundary_rng_counter, &
      has_boundary_rng_counter=has_boundary_rng_counter, adaptive_depth=adaptive_depth, &
      use_upper_panel_fourier=use_upper_panel_fourier, field_available=field_available &
      )
  end subroutine advance_particle_step_impl

  !> 構築済みcandidateとその場評価値からbox eventを解決する。
  subroutine resolve_particle_boundary_candidate( &
    mesh, sim, snapshot, bfield, x0, v0, q, m, dt, x_candidate, v_candidate, hit, result, boundary_contract, &
    boundary_rng_counter, sampled_electric_field &
    )
    type(mesh_type), intent(in) :: mesh
    type(sim_config), intent(in) :: sim
    type(electrostatic_snapshot_type), intent(inout) :: snapshot
    real(dp), intent(in) :: bfield(3), x0(3), v0(3), q, m, dt, x_candidate(3), v_candidate(3)
    type(hit_info), intent(in), optional :: hit
    type(particle_step_result), intent(out) :: result
    type(external_boundary_contract_type), intent(in), optional :: boundary_contract
    integer(i64), intent(in), optional :: boundary_rng_counter(4)
    real(dp), intent(in), optional :: sampled_electric_field(3)
    integer(i64) :: active_boundary_rng_counter(4)
    real(dp) :: active_electric_field(3), x_mid(3)
    logical :: has_boundary_rng_counter
    logical :: field_available

    result = particle_step_result()
    result%x = x0
    result%v = v0
    active_boundary_rng_counter = 0_i64
    has_boundary_rng_counter = present(boundary_rng_counter)
    if (has_boundary_rng_counter) active_boundary_rng_counter = boundary_rng_counter
    result%field_eval_count = 1_i32
    result%collision_query_count = merge(1_i32, 0_i32, present(hit))
    if (.not. valid_particle_step_input(x0, v0, bfield, q, m, dt) .or. &
        .not. all(ieee_is_finite(x_candidate)) .or. .not. all(ieee_is_finite(v_candidate))) then
      result%status = particle_step_invalid_boundary
      return
    end if
    if (present(sampled_electric_field)) then
      active_electric_field = sampled_electric_field
    else
      x_mid = x0 + 0.5_dp*v0*dt
      call project_field_sample_to_box(sim, x_mid)
      call snapshot%eval_local_e(mesh, x_mid, active_electric_field)
      result%field_eval_count = result%field_eval_count + 1_i32
    end if
    if (.not. all(ieee_is_finite(active_electric_field))) then
      result%status = particle_step_invalid_boundary
      return
    end if
    field_available = .true.
    call advance_particle_boundary_crossing( &
      mesh, sim, snapshot, bfield, x0, v0, q, m, dt, x_candidate, v_candidate, active_electric_field, hit, result, &
      boundary_contract, active_boundary_rng_counter, has_boundary_rng_counter, 0_i32, .false., field_available &
      )
  end subroutine resolve_particle_boundary_candidate

  !> box crossing時だけevent用stateを確保し、periodic過多時はremainderを適応分割する。
  subroutine advance_particle_boundary_crossing( &
    mesh, sim, snapshot, bfield, x0, v0, q, m, dt, x_candidate, v_candidate, candidate_electric_field, hit, result, &
    boundary_contract, boundary_rng_counter, has_boundary_rng_counter, adaptive_depth, &
    use_upper_panel_fourier, field_available &
    )
    type(mesh_type), intent(in) :: mesh
    type(sim_config), intent(in) :: sim
    type(electrostatic_snapshot_type), intent(inout) :: snapshot
    real(dp), intent(in) :: bfield(3), x0(3), v0(3), q, m, dt, x_candidate(3), v_candidate(3)
    real(dp), intent(in) :: candidate_electric_field(3)
    type(hit_info), intent(in), optional :: hit
    type(particle_step_result), intent(inout) :: result
    type(external_boundary_contract_type), intent(in), optional :: boundary_contract
    integer(i64), intent(in) :: boundary_rng_counter(4)
    logical, intent(in) :: has_boundary_rng_counter
    integer(i32), intent(in) :: adaptive_depth
    logical, intent(in) :: use_upper_panel_fourier
    logical, intent(inout) :: field_available

    integer(i32), parameter :: max_boundary_events = 8_i32
    type(boundary_event_type) :: event
    type(hit_info) :: remainder_hit
    real(dp) :: x_start(3), v_start(3), x_trial(3), v_trial(3), x_event(3), v_event(3), segment_electric_field(3)
    real(dp) :: dt_segment, dt_remaining
    real(dp) :: redistribution_uniform(3)
    integer(i32) :: query_status, boundary_status, event_count
    logical :: alive, escaped, z_high_outward_event, z_high_barrier_event
    logical :: first_segment, requires_redistribution
    type(external_boundary_contract_type) :: active_boundary_contract

    active_boundary_contract = external_boundary_contract_type()
    select case (trim(sim%open_boundary_model))
    case ('escape')
      active_boundary_contract%ordinary_open_model = external_open_escape
    case ('potential_barrier')
      active_boundary_contract%ordinary_open_model = external_open_potential_barrier
    case default
      result%status = particle_step_invalid_boundary
      return
    end select
    if (present(boundary_contract)) active_boundary_contract = boundary_contract
    x_start = x0
    v_start = v0
    x_trial = x_candidate
    v_trial = v_candidate
    segment_electric_field = candidate_electric_field
    dt_segment = dt
    event_count = 0_i32
    first_segment = .true.

    do
      event = boundary_event_type()
      if (.not. point_strictly_inside_box(sim, x_trial)) then
        call find_first_boundary_event(sim, x_start, x_trial, event, boundary_status)
        if (boundary_status /= boundary_event_ok) then
          result%status = particle_step_invalid_boundary
          return
        end if
      end if
      if (.not. event%has_event) then
        if (first_segment .and. present(hit)) then
          if (hit%has_hit) then
            call accept_particle_hit(v_start, v_trial, hit, result)
            return
          end if
        else
          call query_particle_chord(mesh, sim, x_start, v_start, x_trial, v_trial, result, remainder_hit, query_status)
          if (query_status /= collision_query_ok .or. result%absorbed) return
        end if
        result%x = x_trial
        result%v = v_trial
        return
      end if

      if (first_segment .and. present(hit)) then
        if (hit%has_hit) then
          if (hit%t <= event%fraction + 64.0_dp*epsilon(1.0_dp)*max(1.0_dp, abs(event%fraction))) then
            call accept_particle_hit(v_start, v_trial, hit, result)
            return
          end if
        end if
      end if

      call evaluate_boundary_state( &
        sim, event, x_start, v_start, x_trial, v_trial, bfield, q, m, dt_segment, segment_electric_field, x_event, v_event, &
        boundary_status &
        )
      if (boundary_status /= boundary_event_ok) then
        result%status = particle_step_invalid_boundary
        return
      end if
      if (.not. (first_segment .and. present(hit))) then
        call query_particle_chord(mesh, sim, x_start, v_start, x_event, v_event, result, remainder_hit, query_status)
        if (query_status /= collision_query_ok .or. result%absorbed) return
      end if

      if (event_count >= max_boundary_events) then
        if (periodic_event_can_subdivide(sim, event) .and. adaptive_depth < 12_i32) then
          call advance_periodic_substeps( &
            mesh, sim, snapshot, bfield, x_start, v_start, q, m, dt_segment, result, active_boundary_contract, &
            boundary_rng_counter, has_boundary_rng_counter, adaptive_depth, use_upper_panel_fourier, field_available &
            )
          return
        end if
        result%status = particle_step_multiple_box_events
        return
      end if
      event_count = event_count + 1_i32
      z_high_outward_event = btest(event%face_mask, 5_i32) .and. &
                             event%face_bc(6_i32) == bc_open .and. v_event(3) > 0.0_dp
      z_high_barrier_event = z_high_outward_event .and. &
                             open_face_uses_potential_barrier(6_i32, active_boundary_contract)
      if (z_high_outward_event) then
        result%z_high_outward_event_count = result%z_high_outward_event_count + 1_i32
        result%z_high_outward_normal_kinetic_energy_j_sum = &
          result%z_high_outward_normal_kinetic_energy_j_sum + 0.5_dp*m*v_event(3)*v_event(3)
      end if
      requires_redistribution = event_requires_redistribution_counter(event, active_boundary_contract)
      redistribution_uniform = 0.0_dp
      if (requires_redistribution) then
        if (.not. has_boundary_rng_counter) then
          result%status = particle_step_invalid_boundary
          return
        end if
        call generate_redistribution_uniform( &
          sim%rng_seed, boundary_rng_counter, event_count, redistribution_uniform &
          )
      end if
      dt_remaining = (1.0_dp - event%fraction)*dt_segment
      alive = .true.
      escaped = .false.
      if (event_uses_potential_barrier(event, active_boundary_contract)) then
        call apply_potential_barrier_event( &
          mesh, sim, snapshot, event, active_boundary_contract, q, m, x_event, v_event, alive, escaped, boundary_status, &
          redistribution_uniform, use_upper_panel_fourier, field_available &
          )
      else
        call apply_escape_reflect_periodic_event( &
          sim, event, x_event, v_event, alive, escaped, boundary_status, redistribution_uniform &
          )
      end if
      if (boundary_status /= boundary_event_ok) then
        if (boundary_status == boundary_event_invalid_geometry) then
          result%status = particle_step_invalid_boundary
        else
          result%status = boundary_status
        end if
        return
      end if
      if (z_high_barrier_event) then
        if (alive) then
          result%outer_barrier_return_count = result%outer_barrier_return_count + 1_i32
        else if (escaped) then
          result%outer_barrier_escape_count = result%outer_barrier_escape_count + 1_i32
        end if
      end if
      if (.not. alive) then
        result%x = x_event
        result%v = v_event
        result%escaped_boundary = escaped
        return
      end if
      if (dt_remaining <= 0.0_dp) then
        result%x = x_event
        result%v = v_event
        return
      end if

      x_start = x_event
      v_start = v_event
      dt_segment = dt_remaining
      call build_particle_step_candidate( &
        mesh, sim, snapshot, bfield, x_start, v_start, q, m, dt_segment, x_trial, v_trial, segment_electric_field, &
        use_upper_panel_fourier, field_available &
        )
      result%field_eval_count = result%field_eval_count + 1_i32
      if (.not. field_available) then
        result%status = particle_step_invalid_boundary
        return
      end if
      if (.not. all(ieee_is_finite(x_trial)) .or. .not. all(ieee_is_finite(v_trial))) then
        result%status = particle_step_invalid_boundary
        return
      end if
      first_segment = .false.
    end do
  end subroutine advance_particle_boundary_crossing

  !> 一つのperiodic remainderを二分し、各半stepのbox/collision eventを順番に解く。
  subroutine advance_periodic_substeps( &
    mesh, sim, snapshot, bfield, x0, v0, q, m, dt, result, boundary_contract, boundary_rng_counter, &
    has_boundary_rng_counter, adaptive_depth, use_upper_panel_fourier, field_available &
    )
    type(mesh_type), intent(in) :: mesh
    type(sim_config), intent(in) :: sim
    type(electrostatic_snapshot_type), intent(inout) :: snapshot
    real(dp), intent(in) :: bfield(3), x0(3), v0(3), q, m, dt
    type(particle_step_result), intent(inout) :: result
    type(external_boundary_contract_type), intent(in) :: boundary_contract
    integer(i64), intent(in) :: boundary_rng_counter(4)
    logical, intent(in) :: has_boundary_rng_counter
    integer(i32), intent(in) :: adaptive_depth
    logical, intent(in) :: use_upper_panel_fourier
    logical, intent(inout) :: field_available

    type(particle_step_result) :: first_half, second_half
    integer(i32) :: prior_field_evals, prior_collision_queries
    integer(i32) :: prior_z_high_events, prior_outer_returns, prior_outer_escapes
    real(dp) :: prior_z_high_energy

    prior_field_evals = result%field_eval_count
    prior_collision_queries = result%collision_query_count
    prior_z_high_events = result%z_high_outward_event_count
    prior_outer_returns = result%outer_barrier_return_count
    prior_outer_escapes = result%outer_barrier_escape_count
    prior_z_high_energy = result%z_high_outward_normal_kinetic_energy_j_sum
    if (has_boundary_rng_counter) then
      call advance_particle_step_impl( &
        mesh, sim, snapshot, bfield, x0, v0, q, m, 0.5_dp*dt, first_half, boundary_contract, &
        boundary_rng_counter, adaptive_depth + 1_i32, use_upper_panel_fourier, field_available &
        )
    else
      call advance_particle_step_impl( &
        mesh, sim, snapshot, bfield, x0, v0, q, m, 0.5_dp*dt, first_half, boundary_contract=boundary_contract, &
        adaptive_depth=adaptive_depth + 1_i32, use_upper_panel_fourier=use_upper_panel_fourier, &
        field_available=field_available &
        )
    end if
    if (.not. field_available) then
      result = first_half
      return
    end if
    if (first_half%status /= particle_step_ok .or. first_half%absorbed .or. first_half%escaped_boundary) then
      result = first_half
      result%field_eval_count = prior_field_evals + first_half%field_eval_count
      result%collision_query_count = prior_collision_queries + first_half%collision_query_count
      result%z_high_outward_event_count = prior_z_high_events + first_half%z_high_outward_event_count
      result%outer_barrier_return_count = prior_outer_returns + first_half%outer_barrier_return_count
      result%outer_barrier_escape_count = prior_outer_escapes + first_half%outer_barrier_escape_count
      result%z_high_outward_normal_kinetic_energy_j_sum = prior_z_high_energy + &
                                                          first_half%z_high_outward_normal_kinetic_energy_j_sum
      return
    end if

    if (has_boundary_rng_counter) then
      call advance_particle_step_impl( &
        mesh, sim, snapshot, bfield, first_half%x, first_half%v, q, m, 0.5_dp*dt, second_half, boundary_contract, &
        boundary_rng_counter, adaptive_depth + 1_i32, use_upper_panel_fourier, field_available &
        )
    else
      call advance_particle_step_impl( &
        mesh, sim, snapshot, bfield, first_half%x, first_half%v, q, m, 0.5_dp*dt, second_half, &
        boundary_contract=boundary_contract, adaptive_depth=adaptive_depth + 1_i32, &
        use_upper_panel_fourier=use_upper_panel_fourier, field_available=field_available &
        )
    end if
    result = second_half
    result%field_eval_count = prior_field_evals + first_half%field_eval_count + second_half%field_eval_count
    result%collision_query_count = prior_collision_queries + first_half%collision_query_count + &
                                   second_half%collision_query_count
    result%z_high_outward_event_count = prior_z_high_events + first_half%z_high_outward_event_count + &
                                        second_half%z_high_outward_event_count
    result%outer_barrier_return_count = prior_outer_returns + first_half%outer_barrier_return_count + &
                                        second_half%outer_barrier_return_count
    result%outer_barrier_escape_count = prior_outer_escapes + first_half%outer_barrier_escape_count + &
                                        second_half%outer_barrier_escape_count
    result%z_high_outward_normal_kinetic_energy_j_sum = prior_z_high_energy + &
                                                        first_half%z_high_outward_normal_kinetic_energy_j_sum + &
                                                        second_half%z_high_outward_normal_kinetic_energy_j_sum
  end subroutine advance_periodic_substeps

  !> eventの全faceがperiodicで、適応分割中に乱数境界へ入らない場合だけ再試行する。
  pure logical function periodic_event_can_subdivide(sim, event) result(can_subdivide)
    type(sim_config), intent(in) :: sim
    type(boundary_event_type), intent(in) :: event
    integer(i32) :: axis

    can_subdivide = event%has_event
    if (any(sim%bc_low == bc_redistributed_reflect) .or. any(sim%bc_high == bc_redistributed_reflect)) then
      can_subdivide = .false.
      return
    end if
    do axis = 1_i32, 3_i32
      if (btest(event%face_mask, 2_i32*axis - 2_i32)) then
        if (sim%bc_low(axis) /= bc_periodic) can_subdivide = .false.
      end if
      if (btest(event%face_mask, 2_i32*axis - 1_i32)) then
        if (sim%bc_high(axis) /= bc_periodic) can_subdivide = .false.
      end if
    end do
  end function periodic_event_can_subdivide

  !> 既に選択済みのmesh hitをresultへ反映する。
  subroutine accept_particle_hit(va, vb, hit, result)
    real(dp), intent(in) :: va(3), vb(3)
    type(hit_info), intent(in) :: hit
    type(particle_step_result), intent(inout) :: result

    result%x = hit%pos
    result%v = va + hit%t*(vb - va)
    result%absorbed = .true.
    result%elem_idx = hit%elem_idx
  end subroutine accept_particle_hit

  !> 線分衝突を一度照会し、hit/statusを共通resultへ反映する。
  subroutine query_particle_chord(mesh, sim, xa, va, xb, vb, result, hit, query_status)
    type(mesh_type), intent(in) :: mesh
    type(sim_config), intent(in) :: sim
    real(dp), intent(in) :: xa(3), va(3), xb(3), vb(3)
    type(particle_step_result), intent(inout) :: result
    type(hit_info), intent(out) :: hit
    integer(i32), intent(out) :: query_status

    call find_first_hit(mesh, xa, xb, hit, sim=sim, status=query_status)
    result%collision_query_count = result%collision_query_count + 1_i32
    if (query_status /= collision_query_ok) then
      result%status = query_status
      return
    end if
    if (.not. hit%has_hit) return

    call accept_particle_hit(va, vb, hit, result)
  end subroutine query_particle_chord

  !> Chord上のevent位置と、その接線方向・離散workに整合する速度を返す。
  subroutine evaluate_boundary_state( &
    sim, event, x0, v0, x1, v1, bfield, q, m, dt, electric_field, x_event, v_event, status &
    )
    type(sim_config), intent(in) :: sim
    type(boundary_event_type), intent(in) :: event
    real(dp), intent(in) :: x0(3), v0(3), x1(3), v1(3), bfield(3), q, m, dt, electric_field(3)
    real(dp), intent(out) :: x_event(3), v_event(3)
    integer(i32), intent(out) :: status
    integer(i32) :: axis
    real(dp) :: chord_velocity(3), chord_speed2, event_speed2, speed_scale, speed_tolerance
    logical :: force_free

    status = boundary_event_ok
    x_event = x0 + event%fraction*(x1 - x0)
    do axis = 1_i32, 3_i32
      if (btest(event%face_mask, 2_i32*axis - 2_i32)) x_event(axis) = sim%box_min(axis)
      if (btest(event%face_mask, 2_i32*axis - 1_i32)) x_event(axis) = sim%box_max(axis)
    end do
    chord_velocity = (x1 - x0)/dt
    chord_speed2 = sum(chord_velocity*chord_velocity)
    event_speed2 = sum(v0*v0) + 2.0_dp*(q/m)*dot_product(electric_field, x_event - x0)
    speed_tolerance = 256.0_dp*epsilon(1.0_dp)*max(tiny(1.0_dp), sum(v0*v0), chord_speed2, abs(event_speed2))
    if (.not. ieee_is_finite(chord_speed2) .or. .not. ieee_is_finite(event_speed2) .or. &
        chord_speed2 <= tiny(1.0_dp) .or. event_speed2 <= speed_tolerance) then
      v_event = 0.0_dp
      status = boundary_event_invalid_geometry
      return
    end if
    force_free = q == 0.0_dp .or. (all(electric_field == 0.0_dp) .and. all(bfield == 0.0_dp))
    if (force_free .and. all(v1 == v0)) then
      v_event = v0
      if (.not. all(ieee_is_finite(x_event))) status = boundary_event_invalid_geometry
      return
    end if
    speed_scale = sqrt(event_speed2/chord_speed2)
    v_event = speed_scale*chord_velocity
    if (.not. all(ieee_is_finite(x_event)) .or. .not. all(ieee_is_finite(v_event))) then
      status = boundary_event_invalid_geometry
    end if
  end subroutine evaluate_boundary_state

  !> 単一barrier-open面のpotential-barrier式をevent位置とevent時速度で評価する。
  subroutine apply_potential_barrier_event( &
    mesh, sim, snapshot, event, boundary_contract, q, m, x, v, alive, escaped, status, redistribution_uniform, &
    use_upper_panel_fourier, field_available &
    )
    type(mesh_type), intent(in) :: mesh
    type(sim_config), intent(in) :: sim
    type(electrostatic_snapshot_type), intent(inout) :: snapshot
    type(boundary_event_type), intent(in) :: event
    type(external_boundary_contract_type), intent(in) :: boundary_contract
    real(dp), intent(in) :: q, m
    real(dp), intent(inout) :: x(3), v(3)
    logical, intent(inout) :: alive
    logical, intent(out) :: escaped
    integer(i32), intent(out) :: status
    real(dp), intent(in) :: redistribution_uniform(3)
    logical, intent(in) :: use_upper_panel_fourier
    logical, intent(inout) :: field_available

    type(sim_config) :: action_sim
    type(boundary_event_type) :: action_event
    real(dp) :: phi_boundary, barrier_potential_v, outward_v, kinetic_normal, potential_barrier
    integer(i32) :: axis, face_index, barrier_open_count, ordinary_open_count
    logical :: high_side, reflect_open

    status = boundary_event_ok
    escaped = .false.
    barrier_open_count = 0_i32
    ordinary_open_count = 0_i32
    face_index = 0_i32
    high_side = .false.
    do axis = 1_i32, 3_i32
      if (btest(event%face_mask, 2_i32*axis - 2_i32) .and. event%face_bc(2_i32*axis - 1_i32) == bc_open) then
        if (open_face_uses_potential_barrier(2_i32*axis - 1_i32, boundary_contract)) then
          barrier_open_count = barrier_open_count + 1_i32
          face_index = 2_i32*axis - 1_i32
          high_side = .false.
        else
          ordinary_open_count = ordinary_open_count + 1_i32
        end if
      end if
      if (btest(event%face_mask, 2_i32*axis - 1_i32) .and. event%face_bc(2_i32*axis) == bc_open) then
        if (open_face_uses_potential_barrier(2_i32*axis, boundary_contract)) then
          barrier_open_count = barrier_open_count + 1_i32
          face_index = 2_i32*axis
          high_side = .true.
        else
          ordinary_open_count = ordinary_open_count + 1_i32
        end if
      end if
    end do
    if (ordinary_open_count > 0_i32) then
      alive = .false.
      escaped = .true.
      return
    end if
    if (barrier_open_count > 1_i32) then
      status = particle_step_ambiguous_open_corner
      return
    end if

    action_sim = sim
    action_sim%open_boundary_model = 'escape'
    action_event = event
    if (barrier_open_count == 0_i32) then
      call apply_escape_reflect_periodic_event( &
        action_sim, action_event, x, v, alive, escaped, status, redistribution_uniform &
        )
      return
    end if
    if (.not. ieee_is_finite(q) .or. .not. ieee_is_finite(m) .or. m <= 0.0_dp) then
      status = particle_step_invalid_boundary
      return
    end if

    axis = (face_index + 1_i32)/2_i32
    barrier_potential_v = sim%phi_infty
    if (high_side) then
      if (boundary_contract%barrier_override_high(axis)) then
        barrier_potential_v = boundary_contract%barrier_potential_high_v(axis)
      end if
    else
      if (boundary_contract%barrier_override_low(axis)) then
        barrier_potential_v = boundary_contract%barrier_potential_low_v(axis)
      end if
    end if
    if (use_upper_panel_fourier) then
      call snapshot%eval_upper_panel_fourier_phi(mesh, x, phi_boundary, field_available)
      if (.not. field_available) then
        status = particle_step_invalid_boundary
        return
      end if
    else
      call snapshot%eval_local_phi(mesh, sim, x, phi_boundary)
    end if
    if (.not. ieee_is_finite(phi_boundary)) then
      status = particle_step_invalid_boundary
      return
    end if
    if (high_side) then
      outward_v = v(axis)
    else
      outward_v = -v(axis)
    end if
    kinetic_normal = 0.5_dp*m*outward_v*outward_v
    potential_barrier = q*(barrier_potential_v - phi_boundary)
    if (.not. ieee_is_finite(kinetic_normal) .or. .not. ieee_is_finite(potential_barrier)) then
      status = particle_step_invalid_boundary
      return
    end if
    reflect_open = outward_v > 0.0_dp .and. potential_barrier > 0.0_dp .and. kinetic_normal < potential_barrier
    if (.not. reflect_open) then
      alive = .false.
      escaped = .true.
      return
    end if

    action_event%face_bc(face_index) = bc_reflect
    if (high_side) then
      action_sim%bc_high(axis) = bc_reflect
    else
      action_sim%bc_low(axis) = bc_reflect
    end if
    call apply_escape_reflect_periodic_event( &
      action_sim, action_event, x, v, alive, escaped, status, redistribution_uniform &
      )
  end subroutine apply_potential_barrier_event

  !> eventに含まれるopen面へglobalまたはspecies別のpotential barrierを適用するか返す。
  pure logical function event_uses_potential_barrier(event, boundary_contract) result(uses_barrier)
    type(boundary_event_type), intent(in) :: event
    type(external_boundary_contract_type), intent(in) :: boundary_contract
    integer(i32) :: face

    uses_barrier = .false.
    do face = 1_i32, 6_i32
      if (.not. btest(event%face_mask, face - 1_i32) .or. event%face_bc(face) /= bc_open) cycle
      uses_barrier = open_face_uses_potential_barrier(face, boundary_contract)
      if (uses_barrier) return
    end do
  end function event_uses_potential_barrier

  !> face bit順のopen面がglobalまたはspecies別barrierの対象かを返す。
  pure logical function open_face_uses_potential_barrier(face, boundary_contract) result(uses_barrier)
    integer(i32), intent(in) :: face
    type(external_boundary_contract_type), intent(in) :: boundary_contract
    integer(i32) :: axis

    axis = (face + 1_i32)/2_i32
    uses_barrier = boundary_contract%ordinary_open_model == external_open_potential_barrier
    if (mod(face, 2_i32) == 1_i32) then
      uses_barrier = uses_barrier .or. boundary_contract%barrier_override_low(axis)
    else
      uses_barrier = uses_barrier .or. boundary_contract%barrier_override_high(axis)
    end if
  end function open_face_uses_potential_barrier

  !> event作用が面内再配置用の一意なcounterを必要とするかを返す。
  pure logical function event_requires_redistribution_counter(event, boundary_contract) result(requires)
    type(boundary_event_type), intent(in) :: event
    type(external_boundary_contract_type), intent(in) :: boundary_contract
    integer(i32) :: face
    logical :: has_open, has_ordinary_open, has_redistributed_reflect

    has_open = .false.
    has_ordinary_open = .false.
    has_redistributed_reflect = .false.
    do face = 1_i32, 6_i32
      if (.not. btest(event%face_mask, face - 1_i32)) cycle
      has_open = has_open .or. event%face_bc(face) == bc_open
      if (event%face_bc(face) == bc_open) then
        has_ordinary_open = has_ordinary_open .or. .not. open_face_uses_potential_barrier(face, boundary_contract)
      end if
      has_redistributed_reflect = has_redistributed_reflect .or. &
                                  event%face_bc(face) == bc_redistributed_reflect
    end do
    requires = has_redistributed_reflect .and. .not. has_ordinary_open .and. &
               (.not. has_open .or. event_uses_potential_barrier(event, boundary_contract))
  end function event_requires_redistribution_counter

  !> 粒子event識別子からOpenMP実行順序に依存しない一様乱数を構築する。
  pure subroutine generate_redistribution_uniform(seed, counter, event_index, uniform)
    integer(i32), intent(in) :: seed, event_index
    integer(i64), intent(in) :: counter(4)
    real(dp), intent(out) :: uniform(3)
    integer(i32) :: axis, item
    integer(i64) :: state
    integer(i64), parameter :: hash_modulus = 2147483647_i64

    do axis = 1_i32, 3_i32
      state = counter_hash_mix(104729_i64, int(seed, i64), 37_i64*int(axis, i64))
      state = counter_hash_mix(state, int(axis, i64), 1009_i64)
      do item = 1_i32, 4_i32
        state = counter_hash_mix(state, counter(item), 7919_i64*int(item, i64))
      end do
      state = counter_hash_mix(state, int(event_index, i64), 104729_i64)
      state = counter_hash_mix(state, counter(3), 13007_i64*int(axis, i64))
      state = counter_hash_mix(state, counter(1), 433494437_i64)
      uniform(axis) = (real(state, dp) - 0.5_dp)/real(hash_modulus, dp)
    end do
  end subroutine generate_redistribution_uniform

  !> 31-bit素数体上の非線形写像でtuple成分をcounter hashへ混合する。
  pure integer(i64) function counter_hash_mix(state, value, salt) result(mixed)
    integer(i64), intent(in) :: state, value, salt
    integer(i64), parameter :: hash_modulus = 2147483647_i64
    integer(i64), parameter :: hash_multiplier = 48271_i64

    mixed = modulo(state*state, hash_modulus)
    mixed = modulo(mixed + modulo(hash_multiplier*state, hash_modulus), hash_modulus)
    mixed = modulo(mixed + modulo(value, hash_modulus), hash_modulus)
    mixed = modulo(mixed + modulo(salt, hash_modulus) + 1_i64, hash_modulus)
    if (mixed == 0_i64) mixed = 1_i64
  end function counter_hash_mix

  !> Candidate endpointがboxの全faceからstrictly interiorかを返すfast-path判定。
  pure logical function point_strictly_inside_box(sim, x) result(inside)
    type(sim_config), intent(in) :: sim
    real(dp), intent(in) :: x(3)

    inside = all(x > sim%box_min) .and. all(x < sim%box_max)
  end function point_strictly_inside_box

  !> 粒子stepの有限値・質量・時間刻み入力を検証する。
  pure logical function valid_particle_step_input(x, v, b, q, m, dt) result(valid)
    real(dp), intent(in) :: x(3), v(3), b(3), q, m, dt

    valid = all(ieee_is_finite(x)) .and. all(ieee_is_finite(v)) .and. all(ieee_is_finite(b)) .and. &
            ieee_is_finite(q) .and. ieee_is_finite(m) .and. m > 0.0_dp .and. &
            ieee_is_finite(dt) .and. dt >= 0.0_dp
  end function valid_particle_step_input

end module bem_particle_stepper