sample_flux_weighted_normal_component Subroutine

  subroutine sample_flux_weighted_normal_component(mu, sigma, vn, vmin_normal, sigma_cutoff)
    real(dp), intent(in) :: mu, sigma
    real(dp), intent(in), optional :: vmin_normal
    real(dp), intent(in), optional :: sigma_cutoff
    real(dp), intent(out) :: vn(:)
    integer :: i
    real(dp) :: target, low, high, mid, vmin, cutoff, tail_min, target_tail, tail_high, width, tolerance

    if (size(vn) == 0) return
    vmin = 0.0_dp
    if (present(vmin_normal)) vmin = max(0.0_dp, vmin_normal)
    cutoff = default_velocity_sigma_cutoff
    if (present(sigma_cutoff)) cutoff = sigma_cutoff
    if (cutoff <= 0.0_dp) error stop "sigma_cutoff must be > 0"
    if (sigma <= 0.0_dp) then
      vn = max(mu, vmin)
      return
    end if
    tail_min = flux_weighted_normal_tail(vmin, mu, sigma)
    if (tail_min <= 0.0_dp) then
      vn = vmin
      return
    end if

    do i = 1, size(vn)
      call random_number(target)
      target_tail = tail_min*(1.0_dp - target)
      low = vmin
      width = max(sigma, spacing(max(abs(vmin), sigma)))
      high = vmin + width
      tail_high = flux_weighted_normal_tail(high, mu, sigma)
      do while (tail_high > target_tail)
        width = 2.0_dp*width
        high = vmin + width
        if (.not. ieee_is_finite(high)) error stop "flux-weighted normal inverse survival bracket overflow"
        tail_high = flux_weighted_normal_tail(high, mu, sigma)
      end do
      tolerance = sqrt(epsilon(1.0_dp))*max(sigma, abs(vmin), abs(high), tiny(1.0_dp))
      do while ((high - low) > tolerance)
        mid = 0.5_dp*(low + high)
        if (flux_weighted_normal_tail(mid, mu, sigma) > target_tail) then
          low = mid
        else
          high = mid
        end if
      end do
      vn(i) = 0.5_dp*(low + high)
    end do
  end subroutine sample_flux_weighted_normal_component