evaluate_zhao_density_hat Subroutine

public subroutine evaluate_zhao_density_hat(p, branch, side, phi_hat, phi0_hat, phi_m_hat, n_swe_inf_hat, n_swi_hat, n_swe_f_hat, n_swe_r_hat, n_phe_f_hat, n_phe_c_hat)

Arguments

Type IntentOptional Attributes Name
type(zhao_params_type), intent(in) :: p
character(len=1), intent(in) :: branch
character(len=*), intent(in) :: side
real(kind=dp), intent(in) :: phi_hat
real(kind=dp), intent(in) :: phi0_hat
real(kind=dp), intent(in) :: phi_m_hat
real(kind=dp), intent(in) :: n_swe_inf_hat
real(kind=dp), intent(out) :: n_swi_hat
real(kind=dp), intent(out) :: n_swe_f_hat
real(kind=dp), intent(out) :: n_swe_r_hat
real(kind=dp), intent(out) :: n_phe_f_hat
real(kind=dp), intent(out) :: n_phe_c_hat

Called by

proc~~evaluate_zhao_density_hat~~CalledByGraph proc~evaluate_zhao_density_hat evaluate_zhao_density_hat proc~evaluate_zhao_rho_hat evaluate_zhao_rho_hat proc~evaluate_zhao_rho_hat->proc~evaluate_zhao_density_hat proc~evaluate_zhao_state_from_phi_hat evaluate_zhao_state_from_phi_hat proc~evaluate_zhao_state_from_phi_hat->proc~evaluate_zhao_density_hat proc~build_type_a_branch_from_minimum build_type_a_branch_from_minimum proc~build_type_a_branch_from_minimum->proc~evaluate_zhao_rho_hat proc~sample_monotonic_phi_hat_at_z sample_monotonic_phi_hat_at_z proc~sample_monotonic_phi_hat_at_z->proc~evaluate_zhao_rho_hat proc~sample_zhao_state_at_z sample_zhao_state_at_z proc~sample_zhao_state_at_z->proc~evaluate_zhao_state_from_phi_hat proc~sample_zhao_state_at_z->proc~sample_monotonic_phi_hat_at_z proc~sample_type_a_phi_hat_at_z sample_type_a_phi_hat_at_z proc~sample_zhao_state_at_z->proc~sample_type_a_phi_hat_at_z proc~sample_type_a_phi_hat_at_z->proc~build_type_a_branch_from_minimum proc~sample_zhao_reservoir_state sample_zhao_reservoir_state proc~sample_zhao_reservoir_state->proc~sample_zhao_state_at_z proc~resolve_sheath_injection_context resolve_sheath_injection_context proc~resolve_sheath_injection_context->proc~sample_zhao_reservoir_state proc~init_particle_batch_from_config init_particle_batch_from_config proc~init_particle_batch_from_config->proc~resolve_sheath_injection_context

Source Code

  subroutine evaluate_zhao_density_hat( &
    p, branch, side, phi_hat, phi0_hat, phi_m_hat, n_swe_inf_hat, &
    n_swi_hat, n_swe_f_hat, n_swe_r_hat, n_phe_f_hat, n_phe_c_hat &
    )
    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_hat
    real(dp), intent(out) :: n_swi_hat, n_swe_f_hat, n_swe_r_hat, n_phe_f_hat, n_phe_c_hat

    real(dp) :: arg_ion, s_swe, s_phe, sin_alpha

    sin_alpha = p%n_phe0_m3/p%n_phe_ref_m3
    arg_ion = 1.0d0 - 2.0d0*phi_hat/(p%tau*p%mach*p%mach)
    if (arg_ion <= 0.0d0) error stop 'Zhao ion density argument became non-positive.'
    n_swi_hat = (p%n_swi_inf_m3/p%n_phe_ref_m3)*arg_ion**(-0.5d0)

    select case (branch)
    case ('A')
      s_swe = sqrt(max(0.0d0, (phi_hat - phi_m_hat)/p%tau))
      s_phe = sqrt(max(0.0d0, phi_hat - phi_m_hat))
      n_swe_f_hat = 0.5d0*n_swe_inf_hat*exp(phi_hat/p%tau)*(1.0d0 - erf(s_swe - p%u))
      n_phe_f_hat = 0.5d0*sin_alpha*exp(phi_hat - phi0_hat)*(1.0d0 - erf(s_phe))
      if (trim(side) == 'lower') then
        n_swe_r_hat = 0.0d0
        n_phe_c_hat = sin_alpha*exp(phi_hat - phi0_hat)*erf(s_phe)
      else if (trim(side) == 'upper') then
        n_swe_r_hat = n_swe_inf_hat*exp(phi_hat/p%tau)*(erf(s_swe - p%u) + erf(p%u))
        n_phe_c_hat = 0.0d0
      else
        error stop 'Unknown Type-A Zhao side.'
      end if
    case ('B')
      s_phe = sqrt(max(0.0d0, phi_hat))
      n_swe_f_hat = 0.5d0*n_swe_inf_hat*exp(phi_hat/p%tau)*(1.0d0 + erf(p%u))
      n_swe_r_hat = 0.0d0
      n_phe_f_hat = 0.5d0*sin_alpha*exp(phi_hat - phi0_hat)*(1.0d0 - erf(s_phe))
      n_phe_c_hat = sin_alpha*exp(phi_hat - phi0_hat)*erf(s_phe)
    case ('C')
      s_swe = sqrt(max(0.0d0, (phi_hat - phi0_hat)/p%tau))
      s_phe = sqrt(max(0.0d0, phi_hat - phi0_hat))
      n_swe_f_hat = 0.5d0*n_swe_inf_hat*exp(phi_hat/p%tau)*(1.0d0 - erf(s_swe - p%u))
      n_swe_r_hat = n_swe_inf_hat*exp(phi_hat/p%tau)*(erf(s_swe - p%u) + erf(p%u))
      n_phe_f_hat = 0.5d0*sin_alpha*exp(phi_hat - phi0_hat)*erfc(s_phe)
      n_phe_c_hat = 0.0d0
    case default
      error stop 'Unknown Zhao branch in density evaluation.'
    end select
  end subroutine evaluate_zhao_density_hat