subroutine voigt (x,y,vgt)
C       -------------------------------------------------------
C           Drayson, S. R., JQSRT, v 16., pp 611-614, 1976
C       -------------------------------------------------------
        implicit double precision (a-h,o-z)
        double precision nby2
        logical tru

        common /vgt/ xn(15),h,yn(15),xx(3),hh(3),nby2(19),c(21)
     &   ,b(22),ri(15),d0(25),d1(25),d2(25),d3(25),d4(25),hn(25)
     &   ,co,tru

        if(tru) goto 104

C **       region I. Compute Dawson's Function at mesh points

        tru=.true.
        do 101 i = 1, 15
 101    ri(i) = -i/2.
        do 103 i = 1, 25
           hn(i) = h*(i-0.5)
           co = 4.*hn(i)*hn(i)/25.-2.
           do 102 j = 2, 21
 102       b(j+1) = co*b(j) -b(j-1) + c(j)
           d0(i) = hn(i)*(b(22)-b(21))/5.
           d1(i) = 1.-2.*hn(i)*d0(i)
           d2(i) = (hn(i)*d1(i)+d0(i))/ri(2)
           d3(i) = (hn(i)*d2(i)+d1(i))/ri(3)
 103     d4(i) = (hn(i)*d3(i)+d2(i))/ri(4)
 104     if (x-5.) 105,112,112
 105     if(y-1.) 110,110,106
 106     if(x.gt.1.85*(3.6-y)) goto 112

C **   Region II: Continued fraction. compute number of terms needed.

         if(y.lt.1.45) goto 107
         i = y + y
         goto 108
 107     i = 11.*y
 108     j = x+x+1.85
         imax = xn(j)*yn(i)+0.46         
         imin = min0(16,21-2*imax)

C **   Evaluate continued fraction

         uu = y
         vv = x
         do 109 j = imin, 19
            u = nby2(j)/(uu*uu+vv*vv)
            uu = y + u*uu
 109     vv = x-u*vv
         vgt = uu/(uu*uu+vv*vv)/1.772454
         return

 110     y2 = y*y
         if (x+y.ge.5.) goto 113

C **  Region I. Compute Dawson's function at x from Taylor series

         n = x/h
         dx = x - hn(n+1)
         u = (((d4(n+1)*dx+d3(n+1))*dx+d2(n+1))*dx+d1(n+1))*dx+d0(n+1)
         v = 1.-2.*x*u

C **  Tayloir series expansion about y = 0.0

         vv = exp(y2-x*x)*cos(2.*x*y)/1.128379-y*v
         uu = -y
         imax = 5.+(12.5-x)*0.8*y
         do 111 i = 2, imax, 2
            u = (x*v+u)/ri(i)
            v = (x*u+v)/ri(i+1)
            uu = -uu*y2
 111     vv=vv+v*uu
         vgt = 1.128379*vv
         return
  
 112     y2 = y*y
         if(y .lt. 11.-0.6875*x) goto 113

C **  Region IIIB. 2-point Gauss-Hermite Quadrature
         
         u = x-xx(3)
         v = x+xx(3)
         vgt = y*(hh(3)/(y2+u*u)+hh(3)/(y2+v*v))
         return

C **  Region IIIA. 4-point Gauss-Hermite Quadrature

 113     u = x-xx(1)
         v = x+xx(1)
         uu = x-xx(2)
         vv = x+xx(2)
         vgt = y*(hh(1)/(y2+u*u)+hh(1)/(y2+v*v)+hh(2)/(y2+uu*uu)
     &    +hh(2)/(y2+vv*vv))
         return

         end