!> Zhao query の物理量への変換、残差式、Sagdeev 積分、接続プロファイルと流入応答。
!! 根の探索順・選択方針を持たず、与えられた状態の物理量と成立条件を評価する。
submodule(bem_matching_plane_zhao) bem_matching_plane_zhao_physics
  use bem_sheath_model_core, only: evaluate_zhao_rho_hat, &
                                   zhao_residuals_type_a, zhao_residuals_type_b, zhao_residuals_type_c
  implicit none

  integer, parameter :: rho_quadrature_panels = 256
  integer, parameter :: energy_quadrature_panels = 128
  integer, parameter :: profile_validation_samples = 32
  real(dp), parameter :: profile_negative_tolerance = 1.0e-7_dp
  real(dp), parameter :: profile_endpoint_tolerance = 1.0e-5_dp

contains

  module procedure prepare_matching_zhao_query

  real(dp) :: photoelectron_flux_m2_s, photoelectron_thermal_speed_mps

  params = zhao_params_type()
  photoelectron_temperature_ev = 0.0_dp
  photoelectron_source_density_m3 = 0.0_dp
  status = matching_plane_zhao_invalid_argument
  message = ''
  if (.not. self%initialized) then
    message = 'matching-plane Zhao model is not initialized.'
    return
  end if
  if (.not. all(ieee_is_finite(input))) then
    message = 'matching-plane Zhao query must be finite.'
    return
  end if
  if (input(matching_plane_input_photoelectron_outward_flux) < 0.0_dp .or. &
      input(matching_plane_input_photoelectron_mean_normal_energy) < 0.0_dp .or. &
      input(matching_plane_input_electron_outward_flux) < 0.0_dp .or. &
      input(matching_plane_input_ion_outward_flux) < 0.0_dp) then
    message = 'matching-plane Zhao fluxes and photoelectron energy must be nonnegative.'
    return
  end if

  photoelectron_flux_m2_s = input(matching_plane_input_photoelectron_outward_flux)
  if (photoelectron_flux_m2_s > 0.0_dp) then
    photoelectron_temperature_ev = input(matching_plane_input_photoelectron_mean_normal_energy)
    if (photoelectron_temperature_ev <= 0.0_dp) then
      message = 'positive matching-plane photoelectron flux requires positive mean normal energy.'
      return
    end if
  else
    photoelectron_temperature_ev = self%configured_photoelectron_temperature_ev
  end if

  photoelectron_thermal_speed_mps = sqrt( &
                                    2.0_dp*qe*photoelectron_temperature_ev/self%electron_mass_kg &
                                    )
  if (.not. ieee_is_finite(photoelectron_thermal_speed_mps) .or. &
      photoelectron_thermal_speed_mps <= 0.0_dp) then
    message = 'matching-plane Zhao photoelectron thermal speed is invalid.'
    return
  end if
  photoelectron_source_density_m3 = &
    2.0_dp*sqrt(pi)*photoelectron_flux_m2_s/photoelectron_thermal_speed_mps
  if (.not. ieee_is_finite(photoelectron_source_density_m3) .or. &
      photoelectron_source_density_m3 < 0.0_dp) then
    message = 'matching-plane Zhao photoelectron moment map is invalid.'
    return
  end if

  call prepare_matching_zhao_params( &
    self, photoelectron_temperature_ev, photoelectron_source_density_m3, params, status, message &
    )
  end procedure prepare_matching_zhao_query

  !> `build_zhao_params` の error-stop API を通さず online query を正規化する。
  subroutine prepare_matching_zhao_params( &
    self, photoelectron_temperature_ev, photoelectron_source_density_m3, &
    params, status, message &
    )
    class(matching_plane_zhao_model_type), intent(in) :: self
    real(dp), intent(in) :: photoelectron_temperature_ev, photoelectron_source_density_m3
    type(zhao_params_type), intent(out) :: params
    integer(i32), intent(out) :: status
    character(len=*), intent(out) :: message

    params = zhao_params_type()
    status = matching_plane_zhao_numerical_failure
    message = ''
    params%alpha_rad = 0.5_dp*pi
    params%n_swi_inf_m3 = self%ion_density_m3
    params%n_phe_ref_m3 = self%ion_density_m3
    params%n_phe0_m3 = photoelectron_source_density_m3
    params%photoelectron_population_fraction = 1.0_dp
    params%t_swe_ev = self%electron_temperature_ev
    params%t_phe_ev = photoelectron_temperature_ev
    params%v_d_electron_mps = self%electron_drift_mps
    params%v_d_ion_mps = self%ion_drift_mps
    params%m_i_kg = self%ion_mass_kg
    params%m_e_kg = self%electron_mass_kg
    params%v_swe_th_mps = sqrt(2.0_dp*qe*params%t_swe_ev/params%m_e_kg)
    params%v_phe_th_mps = sqrt(2.0_dp*qe*params%t_phe_ev/params%m_e_kg)
    params%cs_mps = sqrt(qe*params%t_swe_ev/params%m_i_kg)
    params%mach = params%v_d_ion_mps/params%cs_mps
    params%u = params%v_d_electron_mps/params%v_swe_th_mps
    params%tau = params%t_swe_ev/params%t_phe_ev
    params%lambda_d_phe_ref_m = sqrt( &
                                eps0*qe*params%t_phe_ev/(params%n_phe_ref_m3*qe*qe) &
                                )
    if (.not. all(ieee_is_finite([ &
                                 params%alpha_rad, params%n_swi_inf_m3, params%n_phe_ref_m3, params%n_phe0_m3, &
                                 params%photoelectron_population_fraction, params%t_swe_ev, params%t_phe_ev, &
                                 params%v_d_electron_mps, params%v_d_ion_mps, params%m_i_kg, params%m_e_kg, &
                                 params%v_swe_th_mps, params%v_phe_th_mps, params%cs_mps, params%mach, &
                                 params%u, params%tau, params%lambda_d_phe_ref_m &
                                 ])) .or. params%n_phe0_m3 < 0.0_dp .or. &
        min(params%v_swe_th_mps, params%v_phe_th_mps, params%cs_mps, &
            params%mach, params%tau, params%lambda_d_phe_ref_m) <= 0.0_dp) then
      params = zhao_params_type()
      message = 'matching-plane Zhao normalization produced invalid parameters.'
      return
    end if
    status = matching_plane_zhao_ok
  end subroutine prepare_matching_zhao_params

  module procedure encode_matching_unknowns

  y = 0.0_dp
  valid = density_m3 > 0.0_dp .and. params%n_phe_ref_m3 > 0.0_dp .and. params%t_phe_ev > 0.0_dp
  if (.not. valid) return
  select case (branch)
  case ('A')
    valid = phi0_v > 0.0_dp .and. phi_m_v < 0.0_dp
    if (.not. valid) return
    y(1) = log(phi0_v/params%t_phe_ev)
    y(2) = log(-phi_m_v/params%t_phe_ev)
    y(3) = log(density_m3/params%n_phe_ref_m3)
  case ('B')
    valid = phi0_v > 0.0_dp
    if (.not. valid) return
    y(1) = log(phi0_v/params%t_phe_ev)
    y(2) = log(density_m3/params%n_phe_ref_m3)
  case ('C')
    valid = phi0_v < 0.0_dp
    if (.not. valid) return
    y(1) = log(-phi0_v/params%t_phe_ev)
    y(2) = log(density_m3/params%n_phe_ref_m3)
  case default
    valid = .false.
  end select
  valid = valid .and. all(ieee_is_finite(y))
  end procedure encode_matching_unknowns

  module procedure decode_matching_unknowns

  phi0_v = 0.0_dp
  phi_m_v = 0.0_dp
  density_m3 = 0.0_dp
  valid = all(ieee_is_finite(y))
  if (.not. valid .or. y(1) < -50.0_dp .or. y(1) > log(200.0_dp)) then
    valid = .false.
    return
  end if
  select case (branch)
  case ('A')
    if (y(2) < -50.0_dp .or. y(2) > log(200.0_dp) .or. &
        y(3) < -30.0_dp .or. y(3) > log(1.0e6_dp)) then
      valid = .false.
      return
    end if
    phi0_v = params%t_phe_ev*exp(y(1))
    phi_m_v = -params%t_phe_ev*exp(y(2))
    density_m3 = params%n_phe_ref_m3*exp(y(3))
  case ('B')
    if (y(2) < -30.0_dp .or. y(2) > log(1.0e6_dp)) then
      valid = .false.
      return
    end if
    phi0_v = params%t_phe_ev*exp(y(1))
    phi_m_v = phi0_v
    density_m3 = params%n_phe_ref_m3*exp(y(2))
  case ('C')
    if (y(2) < -30.0_dp .or. y(2) > log(1.0e6_dp)) then
      valid = .false.
      return
    end if
    phi0_v = -params%t_phe_ev*exp(y(1))
    phi_m_v = phi0_v
    density_m3 = params%n_phe_ref_m3*exp(y(2))
  case default
    valid = .false.
  end select
  valid = valid .and. all(ieee_is_finite([phi0_v, phi_m_v, density_m3]))
  end procedure decode_matching_unknowns

  module procedure evaluate_charge_residual

  real(dp) :: phi0_v, phi_m_v, density_m3, phi0_hat, phi_m_hat, density_hat
  real(dp) :: raw(3), integral, field_squared, field_residual_scale
  real(dp) :: x3(3), x2(2)
  logical :: integral_ok

  residual = 0.0_dp
  call decode_matching_unknowns( &
    params, branch, y, phi0_v, phi_m_v, density_m3, valid &
    )
  if (.not. valid) return
  phi0_hat = phi0_v/params%t_phe_ev
  phi_m_hat = phi_m_v/params%t_phe_ev
  density_hat = density_m3/params%n_phe_ref_m3
  if (.not. ion_accessible(params, max(phi0_hat, 0.0_dp))) then
    valid = .false.
    return
  end if

  select case (branch)
  case ('A')
    x3 = [phi0_v, phi_m_v, density_m3]
    call zhao_residuals_type_a(params, x3, raw)
    call integrate_matching_rho_hat( &
      params, branch, 'lower', phi_m_hat, phi0_hat, phi0_hat, phi_m_hat, &
      density_hat, integral, integral_ok &
      )
    if (.not. integral_ok) then
      valid = .false.
      return
    end if
    field_squared = -2.0_dp*integral
    if (field_squared < -1.0e-10_dp) then
      valid = .false.
      return
    end if
    field_residual_scale = max(1.0_dp, target_field_hat*target_field_hat)
    residual(1) = raw(1)/params%n_phe_ref_m3
    residual(2) = (max(0.0_dp, field_squared) - target_field_hat*target_field_hat)/field_residual_scale
    residual(3) = raw(3)
  case ('B', 'C')
    x2 = [phi0_v, density_m3]
    if (branch == 'B') then
      call zhao_residuals_type_b(params, x2, raw(1:2))
    else
      call zhao_residuals_type_c(params, x2, raw(1:2))
    end if
    call integrate_matching_rho_hat( &
      params, branch, 'monotonic', phi0_hat, 0.0_dp, phi0_hat, phi_m_hat, &
      density_hat, integral, integral_ok &
      )
    if (.not. integral_ok) then
      valid = .false.
      return
    end if
    field_squared = 2.0_dp*integral
    if (field_squared < -1.0e-10_dp) then
      valid = .false.
      return
    end if
    field_residual_scale = max(1.0_dp, target_field_hat*target_field_hat)
    residual(1) = raw(1)/params%n_phe_ref_m3
    residual(2) = (max(0.0_dp, field_squared) - target_field_hat*target_field_hat)/field_residual_scale
  case default
    valid = .false.
    return
  end select
  valid = all(ieee_is_finite(residual))
  end procedure evaluate_charge_residual

  subroutine integrate_matching_rho_hat( &
    params, branch, side, lower_phi_hat, upper_phi_hat, phi0_hat, phi_m_hat, &
    density_hat, integral, success &
    )
    type(zhao_params_type), intent(in) :: params
    character(len=1), intent(in) :: branch
    character(len=*), intent(in) :: side
    real(dp), intent(in) :: lower_phi_hat, upper_phi_hat, phi0_hat, phi_m_hat, density_hat
    real(dp), intent(out) :: integral
    logical, intent(out) :: success

    real(dp) :: t, phi_hat, jacobian, rho_hat, summand, weight, h
    integer :: point

    integral = 0.0_dp
    success = .false.
    if (.not. all(ieee_is_finite([ &
                                 lower_phi_hat, upper_phi_hat, phi0_hat, phi_m_hat, density_hat &
                                 ])) .or. density_hat <= 0.0_dp) return
    h = 1.0_dp/real(rho_quadrature_panels, dp)
    do point = 0, rho_quadrature_panels
      t = real(point, dp)*h
      phi_hat = lower_phi_hat + (upper_phi_hat - lower_phi_hat)*sin(0.5_dp*pi*t)**2
      jacobian = (upper_phi_hat - lower_phi_hat)*0.5_dp*pi*sin(pi*t)
      if (.not. ion_accessible(params, phi_hat)) return
      call evaluate_zhao_rho_hat( &
        params, branch, side, phi_hat, phi0_hat, phi_m_hat, density_hat, rho_hat &
        )
      if (.not. ieee_is_finite(rho_hat)) return
      summand = rho_hat*jacobian
      if (point == 0 .or. point == rho_quadrature_panels) then
        weight = 1.0_dp
      else if (mod(point, 2) == 0) then
        weight = 2.0_dp
      else
        weight = 4.0_dp
      end if
      integral = integral + weight*summand
    end do
    integral = integral*h/3.0_dp
    success = ieee_is_finite(integral)
  end subroutine integrate_matching_rho_hat

  module procedure validate_matching_root_profile

  real(dp) :: phi0_hat, phi_m_hat, density_hat, phi_hat, fraction
  real(dp) :: integral, field_squared, interface_field_squared, upper_endpoint_field_squared
  real(dp) :: minimum_field_squared, field_squared_scale
  integer :: point
  logical :: integral_ok

  status = matching_plane_zhao_numerical_failure
  message = ''
  root%minimum_field_squared_hat = huge(1.0_dp)
  phi0_hat = root%phi0_v/params%t_phe_ev
  phi_m_hat = root%phi_m_v/params%t_phe_ev
  density_hat = root%ambient_electron_density_m3/params%n_phe_ref_m3
  field_squared_scale = max(1.0_dp, target_field_hat*target_field_hat)
  if (.not. all(ieee_is_finite([ &
                               phi0_hat, phi_m_hat, density_hat, field_squared_scale &
                               ])) .or. density_hat <= 0.0_dp) then
    message = 'matching-plane Zhao profile normalization is invalid.'
    return
  end if

  minimum_field_squared = huge(1.0_dp)
  interface_field_squared = huge(1.0_dp)
  upper_endpoint_field_squared = 0.0_dp
  select case (root%branch)
  case ('A')
    do point = 0, profile_validation_samples
      fraction = real(point, dp)/real(profile_validation_samples, dp)
      phi_hat = phi_m_hat + fraction*(phi0_hat - phi_m_hat)
      call integrate_matching_rho_hat( &
        params, 'A', 'lower', phi_m_hat, phi_hat, phi0_hat, phi_m_hat, &
        density_hat, integral, integral_ok &
        )
      if (.not. integral_ok) then
        message = 'matching-plane Zhao lower profile integration failed.'
        return
      end if
      field_squared = -2.0_dp*integral
      minimum_field_squared = min(minimum_field_squared, field_squared)
      if (point == profile_validation_samples) interface_field_squared = field_squared
    end do
    do point = 0, profile_validation_samples
      fraction = real(point, dp)/real(profile_validation_samples, dp)
      phi_hat = phi_m_hat + fraction*(0.0_dp - phi_m_hat)
      call integrate_matching_rho_hat( &
        params, 'A', 'upper', phi_m_hat, phi_hat, phi0_hat, phi_m_hat, &
        density_hat, integral, integral_ok &
        )
      if (.not. integral_ok) then
        message = 'matching-plane Zhao upper profile integration failed.'
        return
      end if
      field_squared = -2.0_dp*integral
      minimum_field_squared = min(minimum_field_squared, field_squared)
      if (point == profile_validation_samples) upper_endpoint_field_squared = field_squared
    end do
  case ('B', 'C')
    do point = 0, profile_validation_samples
      fraction = real(point, dp)/real(profile_validation_samples, dp)
      phi_hat = phi0_hat + fraction*(0.0_dp - phi0_hat)
      call integrate_matching_rho_hat( &
        params, root%branch, 'monotonic', phi_hat, 0.0_dp, phi0_hat, phi_m_hat, &
        density_hat, integral, integral_ok &
        )
      if (.not. integral_ok) then
        message = 'matching-plane Zhao monotonic profile integration failed.'
        return
      end if
      field_squared = 2.0_dp*integral
      minimum_field_squared = min(minimum_field_squared, field_squared)
      if (point == 0) interface_field_squared = field_squared
    end do
  case default
    message = 'matching-plane Zhao profile has an unknown branch.'
    return
  end select

  if (.not. all(ieee_is_finite([ &
                               minimum_field_squared, interface_field_squared, upper_endpoint_field_squared &
                               ]))) then
    message = 'matching-plane Zhao profile field is non-finite.'
    return
  end if
  root%minimum_field_squared_hat = minimum_field_squared
  if (minimum_field_squared < -profile_negative_tolerance*field_squared_scale) then
    status = matching_plane_zhao_no_physical_solution
    message = 'matching-plane Zhao profile requires an imaginary electric field.'
    return
  end if
  if (abs(interface_field_squared - target_field_hat*target_field_hat) > &
      profile_endpoint_tolerance*field_squared_scale) then
    message = 'matching-plane Zhao profile does not reproduce the interface field.'
    return
  end if
  if (root%branch == 'A' .and. &
      abs(upper_endpoint_field_squared) > profile_endpoint_tolerance*field_squared_scale) then
    message = 'matching-plane Zhao-A upper profile does not reach zero upstream field.'
    return
  end if
  status = matching_plane_zhao_ok
  end procedure validate_matching_root_profile

  module procedure evaluate_root_potential_energy

  real(dp) :: phi0_hat, phi_m_hat, density_hat, energy_hat, segment_energy_hat
  logical :: success

  status = matching_plane_zhao_numerical_failure
  message = ''
  root%potential_energy_j_m2 = huge(1.0_dp)
  phi0_hat = root%phi0_v/params%t_phe_ev
  phi_m_hat = root%phi_m_v/params%t_phe_ev
  density_hat = root%ambient_electron_density_m3/params%n_phe_ref_m3
  if (.not. all(ieee_is_finite([phi0_hat, phi_m_hat, density_hat])) .or. &
      density_hat <= 0.0_dp .or. params%lambda_d_phe_ref_m <= 0.0_dp) then
    message = 'matching-plane Zhao potential-energy normalization is invalid.'
    return
  end if

  energy_hat = 0.0_dp
  select case (root%branch)
  case ('A')
    call integrate_matching_field_energy_hat( &
      params, root%branch, 'lower', phi_m_hat, phi0_hat, phi0_hat, phi_m_hat, &
      density_hat, segment_energy_hat, success &
      )
    if (.not. success) then
      message = 'matching-plane Zhao-A lower potential-energy integral failed.'
      return
    end if
    energy_hat = energy_hat + segment_energy_hat
    call integrate_matching_field_energy_hat( &
      params, root%branch, 'upper', phi_m_hat, 0.0_dp, phi0_hat, phi_m_hat, &
      density_hat, segment_energy_hat, success &
      )
    if (.not. success) then
      message = 'matching-plane Zhao-A upper potential-energy integral failed.'
      return
    end if
    energy_hat = energy_hat + segment_energy_hat
  case ('B', 'C')
    call integrate_matching_field_energy_hat( &
      params, root%branch, 'monotonic', phi0_hat, 0.0_dp, phi0_hat, phi_m_hat, &
      density_hat, energy_hat, success &
      )
    if (.not. success) then
      message = 'matching-plane Zhao monotonic potential-energy integral failed.'
      return
    end if
  case default
    message = 'matching-plane Zhao potential-energy root has an unknown branch.'
    return
  end select
  root%potential_energy_j_m2 = -0.5_dp*eps0*params%t_phe_ev*params%t_phe_ev* &
                               energy_hat/params%lambda_d_phe_ref_m
  if (.not. ieee_is_finite(root%potential_energy_j_m2) .or. root%potential_energy_j_m2 > 0.0_dp) then
    root%potential_energy_j_m2 = huge(1.0_dp)
    message = 'matching-plane Zhao potential energy is invalid.'
    return
  end if
  status = matching_plane_zhao_ok
  end procedure evaluate_root_potential_energy

  subroutine integrate_matching_field_energy_hat( &
    params, branch, side, start_phi_hat, end_phi_hat, phi0_hat, phi_m_hat, &
    density_hat, energy_hat, success &
    )
    type(zhao_params_type), intent(in) :: params
    character(len=1), intent(in) :: branch
    character(len=*), intent(in) :: side
    real(dp), intent(in) :: start_phi_hat, end_phi_hat, phi0_hat, phi_m_hat, density_hat
    real(dp), intent(out) :: energy_hat
    logical, intent(out) :: success

    real(dp) :: t, phi_hat, jacobian, rho_integral, field_squared, summand, weight, h
    real(dp) :: field_squared_scale
    real(dp) :: energy_integrand(0:energy_quadrature_panels)
    integer :: point
    logical :: integral_ok, point_ok(0:energy_quadrature_panels)

    energy_hat = 0.0_dp
    success = .false.
    if (.not. all(ieee_is_finite([ &
                                 start_phi_hat, end_phi_hat, phi0_hat, phi_m_hat, density_hat &
                                 ])) .or. density_hat <= 0.0_dp) return
    h = 1.0_dp/real(energy_quadrature_panels, dp)
    field_squared_scale = max(1.0_dp, phi0_hat*phi0_hat, phi_m_hat*phi_m_hat)
    energy_integrand = 0.0_dp
    point_ok = .false.
    !$omp parallel do default(none) schedule(static) &
    !$omp shared(params,branch,side,start_phi_hat,end_phi_hat,phi0_hat,phi_m_hat,density_hat,h) &
    !$omp shared(field_squared_scale,energy_integrand,point_ok) &
    !$omp private(point,t,phi_hat,jacobian,rho_integral,field_squared,integral_ok)
    do point = 0, energy_quadrature_panels
      t = real(point, dp)*h
      phi_hat = start_phi_hat + (end_phi_hat - start_phi_hat)*sin(0.5_dp*pi*t)**2
      jacobian = (end_phi_hat - start_phi_hat)*0.5_dp*pi*sin(pi*t)
      if (branch == 'A') then
        call integrate_matching_rho_hat( &
          params, branch, side, phi_m_hat, phi_hat, phi0_hat, phi_m_hat, &
          density_hat, rho_integral, integral_ok &
          )
        field_squared = -2.0_dp*rho_integral
      else
        call integrate_matching_rho_hat( &
          params, branch, side, phi_hat, 0.0_dp, phi0_hat, phi_m_hat, &
          density_hat, rho_integral, integral_ok &
          )
        field_squared = 2.0_dp*rho_integral
      end if
      if (.not. integral_ok) cycle
      if (field_squared < -profile_negative_tolerance*field_squared_scale) cycle
      energy_integrand(point) = sqrt(max(0.0_dp, field_squared))*abs(jacobian)
      point_ok(point) = .true.
    end do
    !$omp end parallel do
    if (.not. all(point_ok)) return

    ! Preserve the original Simpson accumulation order for reproducible root ranking.
    do point = 0, energy_quadrature_panels
      summand = energy_integrand(point)
      if (point == 0 .or. point == energy_quadrature_panels) then
        weight = 1.0_dp
      else if (mod(point, 2) == 0) then
        weight = 2.0_dp
      else
        weight = 4.0_dp
      end if
      energy_hat = energy_hat + weight*summand
    end do
    energy_hat = energy_hat*h/3.0_dp
    success = ieee_is_finite(energy_hat) .and. energy_hat >= 0.0_dp
  end subroutine integrate_matching_field_energy_hat

  module procedure compose_matching_response

  real(dp) :: electron_cutoff, electron_term, number_flux_scale
  real(dp) :: electron_access_potential_v, photoelectron_barrier_potential_v

  output = 0.0_dp
  status = matching_plane_zhao_numerical_failure
  message = ''
  select case (root%branch)
  case ('A')
    electron_cutoff = sqrt(max(0.0_dp, -root%phi_m_v/params%t_swe_ev)) - params%u
    electron_access_potential_v = root%phi_m_v
    photoelectron_barrier_potential_v = root%phi_m_v
  case ('B')
    electron_cutoff = -params%u
    electron_access_potential_v = 0.0_dp
    photoelectron_barrier_potential_v = 0.0_dp
  case ('C')
    electron_cutoff = sqrt(max(0.0_dp, -root%phi0_v/params%t_swe_ev)) - params%u
    electron_access_potential_v = 0.0_dp
    photoelectron_barrier_potential_v = 0.0_dp
  case default
    message = 'matching-plane Zhao root has an unknown branch.'
    return
  end select

  electron_term = swe_free_current_term( &
                  params, root%ambient_electron_density_m3, electron_cutoff &
                  )
  number_flux_scale = params%v_phe_th_mps/(2.0_dp*sqrt(pi))
  output(matching_plane_output_matching_potential) = root%phi0_v
  output(matching_plane_output_electron_inward_flux) = number_flux_scale*electron_term
  output(matching_plane_output_ion_inward_flux) = params%n_swi_inf_m3*params%v_d_ion_mps
  output(matching_plane_output_electron_access_potential) = electron_access_potential_v
  output(matching_plane_output_ion_access_potential) = 0.0_dp
  output(matching_plane_output_photoelectron_barrier_potential) = &
    photoelectron_barrier_potential_v
  if (.not. all(ieee_is_finite(output)) .or. &
      any(output(matching_plane_output_electron_inward_flux:matching_plane_output_ion_inward_flux) < 0.0_dp)) then
    output = 0.0_dp
    message = 'matching-plane Zhao response is non-finite or has a negative inward flux.'
    return
  end if
  status = matching_plane_zhao_ok
  end procedure compose_matching_response

  pure logical function ion_accessible(params, phi_hat) result(accessible)
    type(zhao_params_type), intent(in) :: params
    real(dp), intent(in) :: phi_hat

    accessible = params%tau > 0.0_dp .and. params%mach > 0.0_dp .and. &
                 1.0_dp - 2.0_dp*phi_hat/(params%tau*params%mach*params%mach) > 0.0_dp
  end function ion_accessible

end submodule bem_matching_plane_zhao_physics