secondary_emission_yield Function

  pure real(dp) function secondary_emission_yield(model, energy_ev, cos_theta, yield_max, energy_max_ev, const_yield) &
    result(delta)
    character(len=*), intent(in) :: model
    real(dp), intent(in) :: energy_ev, cos_theta, yield_max, energy_max_ev, const_yield
    real(dp) :: cosine, ratio, exponent

    delta = 0.0_dp
    if (model == 'constant') then
      delta = const_yield
    else if (energy_ev > 0.0_dp) then
      cosine = max(0.1_dp, min(1.0_dp, cos_theta))
      ratio = energy_ev/energy_max_ev
      exponent = 2.28_dp*cosine*ratio**1.35_dp
      ! Avoid subtracting two nearly equal numbers at low incident energies.
      if (exponent < 1.0e-5_dp) then
        delta = exponent*(1.0_dp - 0.5_dp*exponent + exponent**2/6.0_dp)
      else
        delta = 1.0_dp - exp(-exponent)
      end if
      delta = 1.114_dp*yield_max/cosine*ratio**(-0.35_dp)*delta
    end if
  end function secondary_emission_yield