FUNCTION VOIGT(V,AGAM)
C     ======================
C
C     Voigt function
C     Procedure after Matta and Reichel, 1971, Math.Comp. 25, 339.
C
      INCLUDE 'IMPLIC.FOR'
      DIMENSION HN(12),EN(12)
      PARAMETER (PI=3.141592653589793D0, M=12, HH=0.5D0, UN=1.D0)
      PARAMETER (PISQ=1.77245385090551D0,PISQ1=UN/PISQ)
      DATA ICOMP /0/
      SAVE EN,HN,PH,HP,ICOMP
C
C     Initialization of auxiliary quantities
C
      IF(ICOMP.EQ.0) THEN
         HP=HH*PISQ1
         PH=PI/HH
         DO I=1,M
            XI=I
            U=XI*XI*HH*HH
            EN(I)=EXP(-U)
            HN(I)=4.D0*U
         END DO
         ICOMP=1
      END IF
C
C     Main term
C
      AGAM1=UN/AGAM
      X=V*AGAM1
      T=0.25D0*AGAM1*AGAM1
      X2=X*X
      X4=4.D0*X2
      S1=UN+X2
      S2=UN-X2
      U0=0.
      DO I=1,M
         S0=HN(I)*T
         U=EN(I)/((S2+S0)*(S2+S0)+X4)
         U0=U0+U*(S1+S0)
      END DO
      S2=UN/S1
      U0=HP*(S2+2.D0*U0)
C
C     Correction term
C
      IF(T.GE.0.25/(PH*PH)) THEN
         U=X/2.D0/T
         A=COS(U)
         B=SIN(U)
         TSQ1=UN/SQRT(T)
         S1=PH*TSQ1
         S2=S1*X
         C=EXP(-S1)-COS(S2)
         D=SIN(S2)
         T4=0.25D0/T
         U=EXP(-X2*T4-S1+T4)*PISQ*TSQ1/(C*C+D*D)
         U0=U0+U*(A*C-B*D)
      END IF
C
      VOIGT=U0*AGAM1*PISQ1
      RETURN
      END