function voigt (x, a)

! Kelly Chance      May 31, 1994

! Evaluation of the Voigt function using algorithm 363 of the Collected
! Algorithms from CACM. Also see "Efficient computation of the complex error
! function", W. Gautschi, SIAM J. Numer. Anal. 7, 187-198 (1970). The accuracy
! is stated as "10 decimal places after the decimal point, or better."
! Certification shows it to be 10 significant figures or better in most
! circumstances.

! The Voigt function is computed for an array of x values and any positive value
! of a, where x = (sigma - sigma0) / xg, xg = 4.30140E-7 * sigma0 * sqrt (T/m);
! a = (Lorentz HWHM)/xg. This routine gives a Voigt profile normalized to area
! xg (in the Doppler limit the line center is always 1/sqrt (pi). Divide the
! output by xg to obtain a Voigt profile normalized to 1.

! The present version closely follows that of D.G. Hummer, including the
! approximation for a near 0.

implicit real*8 (a - h, o - z)
integer capn
real*8 lamda
logical b
data sqrtpi /1.77245385090551d0/, twooverpi /0.63661977236758d0/, fouroverpi /1.27323954473516d0/

! a = 0.
if (a .lt. 1.0d-8) then
    voigt = dexp (-x**2) / sqrtpi
else
  sfac = 1.0d0 - a / 4.29d0
  absx = dabs (x)
  if ((a .lt. 4.29d0) .and. (absx .lt. 5.33d0)) then
    s = sfac * dsqrt (1.d0 - (x / 5.33d0)**2)
    h = 1.6d0 * s
    h2 = 2.0d0 * h
    capn = 6.d0 + 23.0d0 * s
    lamda = h2**capn
    nu = 9.0d0 + 21.0d0 * s
  else
    h = 0.0d0
    capn = 0
    nu = 8
  end if
  b = (h .eq. 0.0d0) .or. (lamda .eq. 0.0d0)
  r1 = 0.0d0
  r2 = 0.0d0
  s1 = 0.0d0
  s2 = 0.0d0
  n = nu
  do in = 1, nu + 1
    np1 = n + 1
    t1 = a + h + dfloat (np1) * r1
    t2 = absx - dfloat (np1) * r2
    c = .5d0 / (t1 * t1 + t2 * t2)
    r1 = c * t1
    r2 = c * t2
    if ((h .gt. 0.0d0) .and. (n .le. capn)) then
      t1 = lamda + s1
      s1 = r1 * t1 - r2 * s2
      s2 = r2 * t1 + r1 * s2
      lamda = lamda / h2
    end if
    n = n - 1
  end do
  if (b) then
    voigt = twooverpi * r1
  else
    voigt = twooverpi * s1
  end if
end if

return
end