!> 1バッチ内で固定する静電場 snapshot。
module bem_electrostatic_snapshot
  use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
  use bem_kinds, only: dp, i32
  use bem_constants, only: eps0
  use bem_types, only: mesh_type, sim_config, bc_periodic
  use bem_field_solver, only: field_solver_type
  use bem_physics_config_types, only: field_physics_config, periodic2_physics_config, panel_kernel_config
  use bem_string_utils, only: lower_ascii
  use bem_periodic_zero_mode_plan, only: periodic_zero_mode_plan_type, periodic_zero_mode_state_type, &
                                         build_periodic_zero_mode_plan, refresh_periodic_zero_mode_state, &
                                         symmetric_vacuum_bottom_field, periodic_zero_mode_ok
  use bem_periodic_zero_mode_eval, only: eval_periodic_zero_mode, zero_mode_trace_plus
  use bem_coulomb_fmm_periodic_nonzero_reference, only: eval_periodic_nonzero_panel_reference
  use bem_coulomb_fmm_periodic_nonzero_upper_vacuum, only: &
    periodic_nonzero_upper_vacuum_plan_type
  use bem_panel_geometry, only: panel_geometry_type, init_panel_geometry, panel_geometry_ok
  use bem_panel_kernel, only: panel_potential_field, panel_side_principal_value
  implicit none
  private

  type, public :: electrostatic_diagnostics_type
    logical :: split_periodic_active = .false.
    logical :: top_reference_available = .false.
    integer(i32) :: top_reference_last_batch = -1_i32
    real(dp) :: top_reference_simulated_time = 0.0_dp
    real(dp) :: top_reference_z_high = 0.0_dp
    integer(i32) :: top_reference_sample_n = 0_i32
    real(dp) :: top_reference_potential_mean = 0.0_dp
    real(dp) :: top_reference_potential_std = 0.0_dp
    real(dp) :: top_reference_potential_min = 0.0_dp
    real(dp) :: top_reference_potential_max = 0.0_dp
    real(dp) :: gauss_residual = 0.0_dp
    character(len=64) :: status = 'not_applicable'
    logical :: periodic_cache_hit = .false.
    integer(i32) :: periodic_operator_build_count = 0_i32
    character(len=128) :: periodic_cache_fingerprint = ''
    character(len=512) :: periodic_cache_path = ''
  end type electrostatic_diagnostics_type

  type, public :: electrostatic_snapshot_type
    type(field_solver_type) :: nonzero_solver
    real(dp) :: prescribed_e(3) = [0.0_dp, 0.0_dp, 0.0_dp]
    real(dp) :: prescribed_phi_origin(3) = [0.0_dp, 0.0_dp, 0.0_dp]
    logical :: use_zero_mode = .false.
    logical :: use_panel_spectral_reference = .false.
    logical :: use_cached_kneq0 = .false.
    logical :: use_upper_panel_fourier_retry = .false.
    real(dp) :: periodic_length(2) = [0.0_dp, 0.0_dp]
    real(dp) :: periodic_origin(2) = [0.0_dp, 0.0_dp]
    integer(i32) :: reference_mode_layers = 0_i32
    integer(i32) :: panel_quadrature_order = 0_i32
    type(periodic_zero_mode_plan_type) :: zero_plan
    type(periodic_zero_mode_state_type) :: zero_state
    type(periodic_nonzero_upper_vacuum_plan_type) :: upper_fourier_plan
    type(periodic2_physics_config) :: periodic_options
    type(electrostatic_diagnostics_type) :: diagnostics
    real(dp) :: gauss_residual = 0.0_dp
    logical :: matching_plane_gauge_active = .false.
    real(dp) :: matching_plane_z = 0.0_dp
    real(dp) :: matching_plane_phi = 0.0_dp
  contains
    procedure :: init => init_electrostatic_snapshot
    procedure :: refresh => refresh_electrostatic_snapshot
    procedure :: eval_local_e => eval_snapshot_local_e
    procedure :: eval_upper_panel_fourier_e => eval_snapshot_upper_panel_fourier_e
    procedure :: eval_upper_panel_fourier_phi => eval_snapshot_upper_panel_fourier_phi
    procedure :: eval_local_phi => eval_snapshot_local_phi
    procedure :: eval_local_phi_without_primary_self => eval_snapshot_local_phi_without_primary_self
    procedure :: compute_mesh_potential => compute_snapshot_mesh_potential
    procedure :: measure_kneq0_potential_step => measure_snapshot_kneq0_potential_step
    procedure :: get_diagnostics => get_snapshot_diagnostics
    procedure :: set_matching_plane_gauge => set_snapshot_matching_plane_gauge
    procedure :: clear_matching_plane_gauge => clear_snapshot_matching_plane_gauge
    procedure :: get_matching_plane_displacement => get_snapshot_matching_plane_displacement
  end type electrostatic_snapshot_type

  interface
    module subroutine eval_snapshot_local_e(self, mesh, position, electric_field)
      class(electrostatic_snapshot_type), intent(inout) :: self
      type(mesh_type), intent(in) :: mesh
      real(dp), intent(in) :: position(3)
      real(dp), intent(out) :: electric_field(3)
    end subroutine eval_snapshot_local_e

    module subroutine eval_snapshot_upper_panel_fourier_e(self, mesh, position, electric_field, available)
      class(electrostatic_snapshot_type), intent(inout) :: self
      type(mesh_type), intent(in) :: mesh
      real(dp), intent(in) :: position(3)
      real(dp), intent(out) :: electric_field(3)
      logical, intent(out) :: available
    end subroutine eval_snapshot_upper_panel_fourier_e

    module subroutine eval_snapshot_upper_panel_fourier_phi(self, mesh, position, potential, available)
      class(electrostatic_snapshot_type), intent(inout) :: self
      type(mesh_type), intent(in) :: mesh
      real(dp), intent(in) :: position(3)
      real(dp), intent(out) :: potential
      logical, intent(out) :: available
    end subroutine eval_snapshot_upper_panel_fourier_phi

    module subroutine eval_snapshot_local_phi(self, mesh, sim, position, potential)
      class(electrostatic_snapshot_type), intent(inout) :: self
      type(mesh_type), intent(in) :: mesh
      type(sim_config), intent(in) :: sim
      real(dp), intent(in) :: position(3)
      real(dp), intent(out) :: potential
    end subroutine eval_snapshot_local_phi

    module subroutine eval_snapshot_local_phi_without_primary_self(self, mesh, sim, element, potential)
      class(electrostatic_snapshot_type), intent(inout) :: self
      type(mesh_type), intent(in) :: mesh
      type(sim_config), intent(in) :: sim
      integer(i32), intent(in) :: element
      real(dp), intent(out) :: potential
    end subroutine eval_snapshot_local_phi_without_primary_self

    module subroutine compute_snapshot_mesh_potential(self, mesh, sim, potential)
      class(electrostatic_snapshot_type), intent(inout) :: self
      type(mesh_type), intent(in) :: mesh
      type(sim_config), intent(in) :: sim
      real(dp), intent(out) :: potential(:)
    end subroutine compute_snapshot_mesh_potential

    module subroutine measure_snapshot_kneq0_potential_step( &
      self, mesh, candidate_charge, max_abs_delta_phi_v, delta_phi_v &
      )
      class(electrostatic_snapshot_type), intent(inout) :: self
      type(mesh_type), intent(in) :: mesh
      real(dp), intent(in) :: candidate_charge(:)
      real(dp), intent(out) :: max_abs_delta_phi_v
      real(dp), intent(out), optional :: delta_phi_v(:)
    end subroutine measure_snapshot_kneq0_potential_step
  end interface

contains

  subroutine init_electrostatic_snapshot(self, mesh, sim, field_config, periodic_config, panel_config)
    class(electrostatic_snapshot_type), intent(inout) :: self
    type(mesh_type), intent(in) :: mesh
    type(sim_config), intent(in) :: sim
    type(field_physics_config), intent(in), optional :: field_config
    type(periodic2_physics_config), intent(in), optional :: periodic_config
    type(panel_kernel_config), intent(in), optional :: panel_config

    self%prescribed_e = sim%e0
    self%prescribed_phi_origin = 0.0_dp
    self%use_zero_mode = .false.
    self%use_panel_spectral_reference = .false.
    self%use_cached_kneq0 = .false.
    self%use_upper_panel_fourier_retry = .false.
    self%gauss_residual = 0.0_dp
    self%matching_plane_gauge_active = .false.
    self%matching_plane_z = 0.0_dp
    self%matching_plane_phi = 0.0_dp
    self%diagnostics = electrostatic_diagnostics_type()

    if (present(periodic_config)) then
      if (trim(lower_ascii(periodic_config%nonzero_mode_backend)) == 'panel_spectral_reference') then
        call init_panel_spectral_snapshot(self, mesh, sim, periodic_config)
        return
      end if
      if (trim(lower_ascii(periodic_config%nonzero_mode_backend)) == 'cached_kneq0') then
        if (.not. present(field_config) .or. .not. present(panel_config)) then
          error stop 'cached_kneq0 requires typed field and panel configuration.'
        end if
        call init_cached_periodic_snapshot(self, mesh, sim, field_config, periodic_config, panel_config)
        return
      end if
    end if

    if (present(field_config) .and. present(periodic_config) .and. present(panel_config)) then
      call self%nonzero_solver%init(mesh, sim, field_config, periodic_config, panel_config)
    else
      call self%nonzero_solver%init(mesh, sim)
    end if
  end subroutine init_electrostatic_snapshot

  subroutine refresh_electrostatic_snapshot(self, mesh)
    class(electrostatic_snapshot_type), intent(inout) :: self
    type(mesh_type), intent(in) :: mesh

    if (.not. self%use_panel_spectral_reference) call self%nonzero_solver%refresh(mesh)
    if (self%use_zero_mode) then
      if (self%matching_plane_gauge_active) then
        call refresh_periodic_zero_mode_state( &
          self%zero_plan, mesh%q_elem, zero_mode_bottom_field(self, mesh%q_elem), &
          self%matching_plane_z, self%matching_plane_phi, self%zero_state &
          )
      else
        call refresh_periodic_zero_mode_state( &
          self%zero_plan, mesh%q_elem, zero_mode_bottom_field(self, mesh%q_elem), &
          self%zero_plan%break_z(1), 0.0_dp, self%zero_state &
          )
      end if
    end if
  end subroutine refresh_electrostatic_snapshot

  !> 周期k=0成分の電位ゲージをmatching planeへ移す。
  !! 非ゼロモードのsolver/cacheは変更しない。次回refreshでもこのゲージを保持する。
  subroutine set_snapshot_matching_plane_gauge(self, mesh, z_high, phi_high)
    class(electrostatic_snapshot_type), intent(inout) :: self
    type(mesh_type), intent(in) :: mesh
    real(dp), intent(in) :: z_high, phi_high

    if (.not. self%use_zero_mode) then
      error stop 'matching-plane gauge requires an active periodic zero mode.'
    end if
    if (.not. ieee_is_finite(z_high) .or. .not. ieee_is_finite(phi_high)) then
      error stop 'matching-plane gauge values must be finite.'
    end if
    if (z_high <= self%zero_plan%break_z(self%zero_plan%nbreak)) then
      error stop 'matching-plane gauge must lie strictly above the charged geometry.'
    end if
    self%matching_plane_gauge_active = .true.
    self%matching_plane_z = z_high
    self%matching_plane_phi = phi_high
    call refresh_periodic_zero_mode_state( &
      self%zero_plan, mesh%q_elem, zero_mode_bottom_field(self, mesh%q_elem), &
      z_high, phi_high, self%zero_state &
      )
  end subroutine set_snapshot_matching_plane_gauge

  subroutine clear_snapshot_matching_plane_gauge(self, mesh)
    class(electrostatic_snapshot_type), intent(inout) :: self
    type(mesh_type), intent(in) :: mesh

    self%matching_plane_gauge_active = .false.
    self%matching_plane_z = 0.0_dp
    self%matching_plane_phi = 0.0_dp
    call self%refresh(mesh)
  end subroutine clear_snapshot_matching_plane_gauge

  !> matching plane直下の法線電束密度 D_H = eps0 E_bottom + Q/A を返す。
  function get_snapshot_matching_plane_displacement(self) result(displacement_c_m2)
    class(electrostatic_snapshot_type), intent(in) :: self
    real(dp) :: displacement_c_m2

    if (.not. self%use_zero_mode) then
      error stop 'matching-plane displacement requires an active periodic zero mode.'
    end if
    if (self%zero_plan%area_xy <= 0.0_dp) then
      error stop 'matching-plane displacement requires positive periodic area.'
    end if
    displacement_c_m2 = eps0*self%zero_state%e_bottom + &
                        self%zero_state%total_charge/self%zero_plan%area_xy
  end function get_snapshot_matching_plane_displacement

  subroutine get_snapshot_diagnostics(self, diagnostics)
    class(electrostatic_snapshot_type), intent(in) :: self
    type(electrostatic_diagnostics_type), intent(out) :: diagnostics

    diagnostics = self%diagnostics
    diagnostics%gauss_residual = self%gauss_residual
  end subroutine get_snapshot_diagnostics

  subroutine init_panel_spectral_snapshot(self, mesh, sim, periodic_config)
    class(electrostatic_snapshot_type), intent(inout) :: self
    type(mesh_type), intent(in) :: mesh
    type(sim_config), intent(in) :: sim
    type(periodic2_physics_config), intent(in) :: periodic_config
    integer(i32) :: status
    character(len=256) :: message

    if (trim(lower_ascii(periodic_config%zero_mode_policy)) /= 'exclude_k0' .or. &
        .not. supported_lower_boundary(periodic_config%lower_boundary_model)) then
      error stop 'panel_spectral_reference requires exclude_k0 and a supported lower boundary model.'
    end if
    if (.not. sim%use_box .or. any(sim%bc_low(1:2) /= bc_periodic) .or. &
        any(sim%bc_high(1:2) /= bc_periodic) .or. sim%bc_low(3) == bc_periodic .or. &
        sim%bc_high(3) == bc_periodic) then
      error stop 'panel_spectral_reference requires x/y periodic and z open.'
    end if

    self%periodic_length = sim%box_max(1:2) - sim%box_min(1:2)
    self%periodic_origin = sim%box_min(1:2)
    if (any(self%periodic_length <= 0.0_dp)) then
      error stop 'panel_spectral_reference requires positive periods.'
    end if
    call build_periodic_zero_mode_plan( &
      mesh, product(self%periodic_length), self%zero_plan, status, message &
      )
    if (status /= periodic_zero_mode_ok) then
      error stop 'panel_spectral_reference zero-mode plan: '//trim(message)
    end if
    self%periodic_options = periodic_config
    call refresh_periodic_zero_mode_state( &
      self%zero_plan, mesh%q_elem, zero_mode_bottom_field(self, mesh%q_elem), &
      self%zero_plan%break_z(1), 0.0_dp, self%zero_state &
      )
    self%reference_mode_layers = periodic_config%reference_mode_layers
    self%panel_quadrature_order = periodic_config%panel_quadrature_order
    self%use_panel_spectral_reference = .true.
    self%use_zero_mode = .true.
    self%diagnostics%split_periodic_active = .true.
    self%diagnostics%status = 'panel_spectral_reference'
  end subroutine init_panel_spectral_snapshot

  subroutine init_cached_periodic_snapshot(self, mesh, sim, field_config, periodic_config, panel_config)
    class(electrostatic_snapshot_type), intent(inout) :: self
    type(mesh_type), intent(in) :: mesh
    type(sim_config), intent(in) :: sim
    type(field_physics_config), intent(in) :: field_config
    type(periodic2_physics_config), intent(in) :: periodic_config
    type(panel_kernel_config), intent(in) :: panel_config
    integer(i32) :: status
    character(len=256) :: message

    if (trim(lower_ascii(periodic_config%zero_mode_policy)) /= 'exclude_k0' .or. &
        .not. supported_lower_boundary(periodic_config%lower_boundary_model)) then
      error stop 'cached_kneq0 requires exclude_k0 and a supported lower boundary model.'
    end if
    if (.not. sim%use_box .or. any(sim%bc_low(1:2) /= bc_periodic) .or. &
        any(sim%bc_high(1:2) /= bc_periodic) .or. sim%bc_low(3) == bc_periodic .or. &
        sim%bc_high(3) == bc_periodic) then
      error stop 'cached_kneq0 snapshot currently requires x/y periodic and z open.'
    end if

    call build_periodic_zero_mode_plan( &
      mesh, product(sim%box_max(1:2) - sim%box_min(1:2)), self%zero_plan, status, message &
      )
    if (status /= periodic_zero_mode_ok) error stop 'cached_kneq0 zero-mode plan: '//trim(message)
    self%periodic_options = periodic_config
    call refresh_periodic_zero_mode_state( &
      self%zero_plan, mesh%q_elem, zero_mode_bottom_field(self, mesh%q_elem), &
      self%zero_plan%break_z(1), 0.0_dp, self%zero_state &
      )
    call self%nonzero_solver%init(mesh, sim, field_config, periodic_config, panel_config)
    self%use_cached_kneq0 = .true.
    self%use_zero_mode = .true.
    self%periodic_length = sim%box_max(1:2) - sim%box_min(1:2)
    self%periodic_origin = sim%box_min(1:2)
    self%reference_mode_layers = periodic_config%reference_mode_layers
    self%panel_quadrature_order = periodic_config%panel_quadrature_order
    self%diagnostics%split_periodic_active = .true.
    self%diagnostics%status = 'cached_kneq0'
    self%diagnostics%periodic_cache_hit = self%nonzero_solver%fmm_core_plan%periodic_cache_hit
    self%diagnostics%periodic_operator_build_count = &
      self%nonzero_solver%fmm_core_plan%periodic_operator_build_count
    self%diagnostics%periodic_cache_fingerprint = &
      self%nonzero_solver%fmm_core_plan%periodic_cache_fingerprint
    self%diagnostics%periodic_cache_path = self%nonzero_solver%fmm_core_plan%periodic_cache_path
    if (trim(lower_ascii(sim%multiple_box_events_retry_backend)) == 'upper_panel_fourier') then
      call self%upper_fourier_plan%build( &
        mesh, self%periodic_length(1), self%periodic_length(2), self%reference_mode_layers &
        )
      self%use_upper_panel_fourier_retry = .true.
    end if
  end subroutine init_cached_periodic_snapshot

  pure real(dp) function zero_mode_bottom_field(self, charge) result(field)
    class(electrostatic_snapshot_type), intent(in) :: self
    real(dp), intent(in) :: charge(:)

    select case (trim(lower_ascii(self%periodic_options%lower_boundary_model)))
    case ('e_bottom_zero')
      field = 0.0_dp
    case ('symmetric_vacuum')
      field = symmetric_vacuum_bottom_field(self%zero_plan, charge)
    case default
      error stop 'unsupported periodic2 lower boundary model.'
    end select
  end function zero_mode_bottom_field

  pure logical function supported_lower_boundary(model) result(supported)
    character(len=*), intent(in) :: model

    select case (trim(lower_ascii(model)))
    case ('e_bottom_zero', 'symmetric_vacuum')
      supported = .true.
    case default
      supported = .false.
    end select
  end function supported_lower_boundary

end module bem_electrostatic_snapshot