NB.  Voigt Function, normalized to unity.  
NB.  (Multiply by pi to get H)             
NB.  The Method is that of A. Riechel      
NB.  1968 J.Q.S.R.T., Vol 8, P 1601. 
NB.  The frequency from line center in Doppler widths
NB.  is denoted by x = (nu-nu_0)%Delta nu_D , where
NB.  Delta nu_D = (lambda_0/c) sqrt(2kT/M). The parameter
NB.  a = gamma/(4*pi*Delta nu_D), the ratio of the damping
NB.  parameter gamma to 4 pi times the Doppler width.        
NB.  Usage: "a voigt x"  (x may be a vector)

voigt=: 4 : 'x vvvv"0 y' NB. y item by item
vvvv=: dyad define
if. 7 >: (u=. |y) do.       
  f=. (1p1*x)%~ ^-u*u
  f1=. +:x*x
  u1=. 1p0.5* x
  sum=. u1+ u2=. -f1
  n=. 2
  while. n<: 40 do.
    u3=. f1*u1%n
    u1=. u2
    sum=. sum+ u2=. u3
    if. (|u2%sum)< 1e_11 do. break. end.
    n=. >:n
  end.
  fct=. n=. 1
  usq=. *:u
  s=. sum
  while. n<:150 do.
    fct=. usq*fct%n
    tm=. fct* s=. f1*(1-s)%_1++:n
    sum=. sum+ tm
    if. (|tm%sum)< 1e_11 do. break. end.
    n=. >:n
  end.
f*sum
else.             
NB. ------------- far wings: u > 7 -------
  v1=. - sign=. 1
  v2=. - sum=. f1=. x%u
  f2=. -: %u*u
  n=. 3
  while. n<:55 do.
    v3=. (f1*v2)-(n-2)*f2*v1
    v4=. (f1*v3)-(n-1)*f2*v2
    n=. n+2
    if. ((|v4)>|v2) do. break. end.
    sum=. sum+ sign*v4
    if. (|v4%sum) < 1e_10 do. break. end.
    'v1 v2'=. v3;v4
    sign=. -sign
  end. 
  sum%1p1*u 
end.
)