evaluate_zhao_state_from_phi_hat Subroutine

  subroutine evaluate_zhao_state_from_phi_hat(p, branch, side, phi_hat, phi0_hat, phi_m_hat, n_swe_inf_m3, z_m_m, state)
    type(zhao_params_type), intent(in) :: p
    character(len=1), intent(in) :: branch
    character(len=*), intent(in) :: side
    real(dp), intent(in) :: phi_hat, phi0_hat, phi_m_hat, n_swe_inf_m3, z_m_m
    type(zhao_local_state_type), intent(out) :: state

    real(dp) :: n_swi_hat, n_swe_f_hat, n_swe_r_hat, n_phe_f_hat, n_phe_c_hat
    real(dp) :: arg_ion, a_swe, a_phe

    call evaluate_zhao_density_hat( &
      p, branch, side, phi_hat, phi0_hat, phi_m_hat, n_swe_inf_m3/p%n_phe_ref_m3, &
      n_swi_hat, n_swe_f_hat, n_swe_r_hat, n_phe_f_hat, n_phe_c_hat &
      )

    arg_ion = 1.0d0 - 2.0d0*phi_hat/(p%tau*p%mach*p%mach)
    if (arg_ion <= 0.0d0) error stop 'Zhao local ion energy argument became non-positive.'

    select case (branch)
    case ('A')
      a_swe = sqrt(max(0.0d0, (phi_hat - phi_m_hat)/p%tau))
      a_phe = sqrt(max(0.0d0, phi_hat - phi_m_hat))
      state%swe_reflected_active = trim(side) == 'upper'
      state%phe_captured_active = trim(side) == 'lower'
    case ('B')
      a_swe = 0.0d0
      a_phe = sqrt(max(0.0d0, phi_hat))
      state%swe_reflected_active = .false.
      state%phe_captured_active = .true.
    case ('C')
      a_swe = sqrt(max(0.0d0, (phi_hat - phi0_hat)/p%tau))
      a_phe = sqrt(max(0.0d0, phi_hat - phi0_hat))
      state%swe_reflected_active = .true.
      state%phe_captured_active = .false.
    case default
      error stop 'Unknown Zhao branch in local state evaluation.'
    end select

    state%branch = branch
    state%side = side
    state%z_m_m = z_m_m
    state%phi_hat = phi_hat
    state%phi_v = phi_hat*p%t_phe_ev
    state%n_swi_m3 = n_swi_hat*p%n_phe_ref_m3
    state%n_swe_f_m3 = n_swe_f_hat*p%n_phe_ref_m3
    state%n_swe_r_m3 = n_swe_r_hat*p%n_phe_ref_m3
    state%n_phe_f_m3 = n_phe_f_hat*p%n_phe_ref_m3
    state%n_phe_c_m3 = n_phe_c_hat*p%n_phe_ref_m3
    state%electron_source_density_m3 = n_swe_inf_m3*exp(phi_hat/p%tau)
    state%vcut_swe_mps = a_swe*p%v_swe_th_mps
    state%vcut_phe_mps = a_phe*p%v_phe_th_mps
    state%v_i_mps = p%v_d_ion_mps*sqrt(arg_ion)
  end subroutine evaluate_zhao_state_from_phi_hat