subroutine resolve_sheath_injection_context(cfg, ctx)
type(app_config), intent(in) :: cfg
type(sheath_injection_context), intent(out) :: ctx
integer(i32) :: electron_idx, ion_idx, photo_idx
type(particle_species_spec) :: spec_e_cfg, spec_i_cfg, spec_p_cfg
type(sheath_model_species) :: spec_e, spec_i, spec_p
type(zhao_params_type) :: params
type(zhao_local_state_type) :: local_state
real(dp) :: phi0_v, phi_m_v, n_swe_inf_m3
real(dp) :: n_e_inf_m3, n_i_inf_m3, n_phe_ref_m3
real(dp) :: t_swe_ev, t_phe_ev, v_d_electron_mps, v_d_ion_mps
real(dp) :: reference_coordinate, reference_inward_normal(3)
integer :: reference_axis
character(len=32) :: model
character(len=1) :: branch
ctx = sheath_injection_context()
model = trim(lower_ascii(cfg%sim%sheath_injection_model))
ctx%model = model
if (model == 'none') return
call detect_sheath_species(cfg, electron_idx, ion_idx, photo_idx)
spec_e_cfg = cfg%particle_species(electron_idx)
spec_i_cfg = cfg%particle_species(ion_idx)
spec_e = to_sheath_model_species(spec_e_cfg)
spec_i = to_sheath_model_species(spec_i_cfg)
ctx%enabled = .true.
ctx%electron_species = electron_idx
ctx%ion_species = ion_idx
ctx%reference_face = trim(lower_ascii(spec_e_cfg%inject_face))
call resolve_sheath_reference_plane( &
cfg%sim, ctx%reference_face, reference_axis, reference_coordinate, reference_inward_normal &
)
ctx%reference_axis = int(reference_axis, i32)
ctx%reference_coordinate = reference_coordinate
ctx%reference_inward_normal = reference_inward_normal
n_e_inf_m3 = spec_e%number_density_m3
n_i_inf_m3 = spec_i%number_density_m3
t_swe_ev = spec_e%temperature_k*k_boltzmann/qe
v_d_electron_mps = resolve_species_drift_speed( &
spec_e, trim(lower_ascii(cfg%sim%sheath_electron_drift_mode)), &
reference_inward_normal)
v_d_ion_mps = resolve_species_drift_speed( &
spec_i, trim(lower_ascii(cfg%sim%sheath_ion_drift_mode)), &
reference_inward_normal)
if (v_d_electron_mps <= 0.0d0) error stop 'sheath electron drift must point inward.'
if (v_d_ion_mps <= 0.0d0) error stop 'sheath ion drift must point inward.'
select case (model)
case ('floating_no_photo')
call solve_no_photo_floating_potential(spec_e, spec_i, reference_inward_normal, phi0_v)
ctx%branch = 'N'
ctx%phi0_v = phi0_v
ctx%phi_m_v = phi0_v
ctx%n_swe_inf_m3 = n_e_inf_m3
ctx%electron_number_density_m3 = n_e_inf_m3
ctx%electron_vmin_normal = sqrt(max(0.0d0, 2.0d0*qe*max(0.0d0, -phi0_v))/abs(spec_e%m_particle))
return
case ('zhao_auto', 'zhao_a', 'zhao_b', 'zhao_c')
if (photo_idx <= 0_i32) then
error stop 'Zhao sheath injection requires one enabled negative-q photo_raycast species.'
end if
spec_p_cfg = cfg%particle_species(photo_idx)
spec_p = to_sheath_model_species(spec_p_cfg)
if (abs(spec_p_cfg%q_particle) <= 0.0d0) error stop 'Zhao photoelectron species must have non-zero q_particle.'
ctx%photo_species = photo_idx
ctx%has_photo_species = .true.
t_phe_ev = spec_p%temperature_k*k_boltzmann/qe
n_phe_ref_m3 = cfg%sim%sheath_photoelectron_ref_density_cm3*1.0d6
call build_zhao_params( &
cfg%sim%sheath_alpha_deg, n_i_inf_m3, n_phe_ref_m3, t_swe_ev, t_phe_ev, v_d_electron_mps, v_d_ion_mps, &
spec_i%m_particle, abs(spec_e%m_particle), params &
)
call solve_zhao_unknowns(model, params, phi0_v, phi_m_v, n_swe_inf_m3, branch)
ctx%branch = branch
ctx%phi0_v = phi0_v
ctx%phi_m_v = phi_m_v
ctx%n_swe_inf_m3 = n_swe_inf_m3
ctx%electron_number_density_m3 = n_swe_inf_m3
ctx%electron_vmin_normal = zhao_electron_vmin_normal(branch, phi0_v, phi_m_v, abs(spec_e%m_particle))
ctx%photo_vmin_normal = zhao_photo_vmin_normal(branch, phi0_v, phi_m_v, abs(spec_p%m_particle))
ctx%photo_emit_current_density_a_m2 = zhao_photo_emit_current_density(branch, phi0_v, phi_m_v, spec_p_cfg%q_particle, params)
if (cfg%sim%has_sheath_reference_coordinate) then
call sample_zhao_reservoir_state(cfg%sim, spec_e_cfg, ctx, params, branch, phi0_v, phi_m_v, n_swe_inf_m3, local_state)
ctx%has_local_reservoir_profile = .true.
ctx%reservoir_phi_v = local_state%phi_v
ctx%electron_number_density_m3 = local_state%electron_source_density_m3
if (local_state%swe_reflected_active) then
ctx%electron_vmin_normal = 0.0d0
else
ctx%electron_vmin_normal = local_state%vcut_swe_mps
end if
ctx%ion_number_density_m3 = local_state%n_swi_m3
ctx%ion_normal_speed_mps = local_state%v_i_mps
end if
case default
error stop 'Unknown sim.sheath_injection_model in runtime.'
end select
end subroutine resolve_sheath_injection_context