!  Kuntz, M.,1997 : A new implementation of the humlicek algorithm for the calculation of the voigt profile function, 
!  Journal of Quantitative Spectroscopy and Radiative Transfer, Volume 57, Issue 6,  June 1997, Pages 819-824 
!  URL : http://www.sciencedirect.com/science?_ob=MImg&_imagekey=B6TVR-4BJW9BP-1-W&_cdi=5541&_orig=search&
!_coverDate=06%2F15%2F2004&_sk=999139997&view=c&wchp=dGLbVzb-zSkzV&_acct=C000009418&_version=1&_userid=117185&md5=9c5621a94c2ab05c598310b143bade13&ie=f.pdf 
!
! Ruyten, W., 2004 : Comment on "A new implementation of the Humlicek algorithm for the calculation of the Voigt profile function" 
! by M. Kuntz [JQSRT 57(6) (1997) 819-824]
! Journal of Quantitative Spectroscopy and Radiative Transfer, Volume 86, Issue 2,  15 June 2004, Pages 231-233
! URL ; http://www.sciencedirect.com/science?_ob=MImg&_imagekey=B6TVR-3SPKTGS-8-3&_cdi=5541&_orig=search&
!_coverDate=06%2F30%2F1997&_sk=999429993&view=c&wchp=dGLbVzb-zSkzV&_acct=C000009418&_version=1&_userid=117185&md5=1468dff77dfd945efd68e0a336c997d7&ie=f.pdf
!

function voigt(x, y) result(out)
   implicit none

   INTRINSIC EXP
   real(8), intent(in) :: x, y
   real(8) :: out

   REAL(8) :: a1, b1, a2, b2
   REAL(8) :: a3, b3, c3, d3, a4, b4, c4, d4
   REAL(8) :: a5, b5, c5, d5, e5, a6, b6, c6, d6, e6
   REAL(8) :: a7, b7, c7, d7, e7, f7, g7, h7, o7, p7, q7, r7, s7, t7
   REAL(8) :: a8, b8, c8, d8, e8, f8, g8, h8, o8, p8, q8, r8, s8, t8


   if ( abs(x) + y  >= 15.0 ) then
      a1 =   0.2820948 * y  + 0.5641896 * y**3 
      b1 =   0.5641896 * y
      a2 =   0.5 + y**2 + y**4
      b2 =  -1.0 + 2.0*y**2
      out =   (a1 + b1*x**2       ) & 
           & / (a2 + b2*x**2 + x**4)

   elseif ( (abs(x) + y >= 5.5) .and. (abs(x) + y < 15.0) ) then
      a3 = 1.05786 * y + 4.65456 * y**3 + 3.10304 * y**5 + 0.56419 * y**7   
      b3 = 2.962   * y + 0.56419 * y**3 + 1.69257 * y**5
      c3 = 1.69257 * y**3 - 2.53885 * y
      d3 = 0.56419 * y
      a4 = 0.5625 + 4.5 * y**2 + 10.5 * y**4 + 6.0*y**6 + y**8
      b4 = -4.5 + 9.0 * y**2 + 6.0 * y**4 + 4.0 * y**6
      c4 = 10.5 - 6.0 * y**2 + 6.0 * y**4
      d4 = -6.0 + 4.0 * y**2
      out = (a3 + b3 * x**2 + c3 * x**4 + d3 * x**6)   &
           & /(a4 + b4 * x**2 + c4 * x**4 + d4 * x**6 + x**8)
   elseif ( ( abs(x) + y < 5.5) .and. (y >= 0.195*abs(x) - 0.176) ) then
      a5 =   272.102      + 973.778*y    + 1629.76*y**2 + 1678.33*y**3 &
           & +  1174.8*y**4  + 581.746*y**5 + 204.510*y**6 + 49.5213*y**7 &
           & + 7.55895*y**8 +  0.564224*y**9  
      b5 =  -60.5644 - 2.34403*y + 220.843*y**2 +  336.364*y**3 & 
           & + 247.198*y**4 +  100.705*y**5 & 
           & + 22.6778*y**6 + 2.25689*y**7
      c5 =  4.58029 +  18.546*y + 42.5683*y**2 + 52.8454*y**3  & 
           & + 22.6798*y**4 + 3.38534*y**5
      d5 =  -0.128922  +  1.66203*y +  7.56186*y**2 + 2.25689*y**3  
      e5 = 0.000971457 + 0.564224*y
      a6 = 272.102 +  1280.83*y + 2802.87*y**2 +  3764.97*y**3  & 
           & + 3447.63*y**4 + 2256.98*y**5 +  1074.41*y**6 + 369.199*y**7 &
           & + 88.2674*y**8  +  13.3988*y**9 + y**10
      b6 =  211.678  + 902.306*y +  1758.34*y**2 + 2037.31*y**3 &
           & + 1549.68*y**4 + 793.427*y**5 + 266.299*y**6 + 53.5952*y**7 & 
           & + 5.0*y**8 
      c6 = 78.866 +  308.186*y + 497.302*y**2 + 479.258*y**3 + 269.292*y**4 &
           & + 80.3928*y**5 +  10.0*y**6
      d6 = 22.0353 +  55.0293*y + 92.7568*y**2 +  53.5952*y**3 +  10.0*y**4
      e6 = 1.49645 +  13.3988*y +  5.0*y**2 
      out =  ( a5 +b5*x**2 + c5*x**4 + d5*x**6 + e5*x**8) &
           & /(a6 + b6*x**2 + c6*x**6 + d6*x**8 + e6*x**8 + x**10) 
   elseif ( (abs(x) + y < 5.5) .and.( y > 0.195*abs(x) - 0.176) ) then
      a7 =  1.16028E+9*y -  9.86604E+8*y**3 + 4.56662E+8*y**5 & 
           & -  1.53575E+8*y**7 + 4.08168E+7*y**9 - 9.69463E+6*y**11 & 
           & +  1.6841E+6*y**13 - 320772.0*y**15 + 40649.2*y**17 -  5860.68*y**19 &
           & +  571.687*y**21 - 72.9359*y**23 + 2.35944*y**25 - 0.56419*y**27  
      b7 = -5.60505E+8*y  -  9.85386E+8*y**3 + 8.06985e+8*y**5 & 
           & -  2.91876e+8*y**7 + 8.64829e+7*y**9 - 7.72359e+6*y**11 & 
           & + 3.59915e+6*y**13 - 234417.0*y**15 + 45251.3*y**17 &
           & -  2269.19*y**19 - 234.143*y**21 + 23.0312*y**23 - 7.33447*y**25
      c7 =  -6.51523e+8*y  + 2.47157e+8*y**3 + 2.94262e+8*y**5  & 
           & -  2.04467e+8*y**7 + 2.29302e+7*y**9 -  2.3818e+7*y**11 &
           & +  576054.0*y**13 +  98079.1*y**15 -  25338.3*y**17 &
           & +  1097.77*y**19 + 97.6203*y**21 -  44.0068*y**23
      d7 =  -2.63894e+8*y  + 2.70167e+8*y**3 -  9.96224e+7*y**5 & 
           & -  4.15013e+7*y**7 + 3.83112e+7*y**9 + 2.2404e+6*y**11 & 
           & -  303569.0*y**13 - 66431.2*y**15 + 8381.97*y**17 + 228.563*y**19  &
           & -  161.358*y**21
      e7 =  -6.31771e+7*y + 1.40677e+8*y**3 +  5.56965e+6*y**5 & 
           & + 2.46201e+7*y**7 + 468142.0*y**9 - 1.003e+6*y**11 - 66212.1*y**13 &
           & + 23507.6*y**15  +  296.38*y**17 - 403.396*y**19
      f7 = - 1.69846e+7*y + 4.07382e+6*y**3 - 3.32896e+7*y**5 - 1.93114e+6*y**7 &
           & -  934717.0*y**9 + 8820.94*y**11 + 37544.8*y**13 &
           &  +  125.591*y**15 -  726.113*y**17
      g7 = - 1.23165e+6*y + 7.52883e+6*y**3 - 900010.0*y**5 - 186682.0*y**7  &
           & + 79902.5*y**9 + 37371.9*y**11 -  260.198*y**13 -  968.15*y**15
      h7 =  -610622.0*y + 86407.6*y**3 + 153468.0*y**5 + 72520.9*y**7 & 
           & + 23137.1*y**9 - 571.645*y**11 -  968.15*y**13
      o7 =  -23586.5*y  + 49883.8*y**3 + 26538.5*y**5 +  8073.15*y**7 &
           & -  575.164*y**9 -  726.113*y**11
      p7 = -8009.1*y + 2198.86*y**3 + 953.655*y**5 - 352.467*y**7 - 403.396*y**9
      q7 =  -622.056*y  -  271.202*y**3 -  134.792*y**5 -  161.358*y**7
      r7 =  - 77.0535*y -  29.7896*y**3 -  44.0068*y**5 
      S7 =  -2.92264*y  -  7.33447*y**3
      t7 =  -0.56419*y
      a8 =  1.02827e+9 - 1.5599e+9*y**2  +  1.17022e+9*y**4 - 5.79099e+8*y**6 & 
           & + 2.11107e+8*y**8  -  6.11148e+7*y**10 +  1.44647e+7*y**12  &
           & -  2.85721e+6*y**14 + 483737.0*y**16 -  70946.1*y**18 + 9504.65*y**20 & 
           & -  955.194*y**22 +  126.532*y**24 -  3.68288*y**26 + 1.0*y**28
      b8 = 1.5599e+9 -  2.28855e+9*y**2 +  1.66421e+9*y**4 -  7.53828e+8*y**6 &
           &  + 2.89676e+8*y**8 - 7.01358e+7*y**10 +  1.39465e+7*y**12 &
           &  - 2.84954e+6*y**14 + 498334.0*y**16 - 55600.0*y**18 + 3058.26*y**20 & 
           &  + 533.254*y**22  -  40.5117*y**24 +  14.0*y**26  
      c8 =  1.17022e+9 - 1.66421e+9*y**2 + 1.06002e+9*y**4 - 6.60078e+8*y**6 & 
           & + 6.33496e+7*y**8 - 4.60396e+7*y**10 +  1.4841e+7*y**12 &
           & - 1.06352e+6*y**14 - 217801.0*y**16 + 48153.3*y**18 - 1500.17*y**20 & 
           & - 198.876*y**22 + 91.0*y**24
      d8 = 5.79099e+8 - 7.53828e+8*y**2 + 6.60078e+8*y**4 + 5.40367e+7*y**6 & 
           & + 1.99846e+8*y**8 -  6.87656e+6*y**10 -  6.89002e+6*y**12 &
           & + 280428.0*y**14 + 161461.0*y**16 - 16493.7*y**18 - 567.164*y**20 & 
           & + 364.0*y**22
      e8 = 2.11107e+8 -  2.89676e+8*y**2 + 6.33496e+7*y**4 -  1.99846e+8*y**6 &
           & -  5.01017e+7*y**8 - 5.25722e+6*y**10 +  1.9547e+6*y**12 &
           & + 240373.0*y**14  -  55582.0*y**16 -  1012.79*y**18 +  1001.0*y**20
      f8 = 6.11148e+7 -  7.01358e+7*y**2 + 4.60396e+7*y**4 -  6.87656e+6*y**6 &
           & +  5.25722e+6*y**8 + 3.04316e+6*y**10 +  123052.0*y**12  &
           &  - 106663.0*y**14  -  1093.82*y**16 + 2002.0*y**18
      g8 =  1.44647e+7 -  1.39465e+7*y**2 +  1.4841e+7*y**4 + 6.89002e+6*y**6 &
           &  +  1.9547e+6*y**8  -  123052.0*y**10 -  131337.0*y**12  &
           & -  486.14*y**14  +  3003.0*y**16
      h8 = 2.85721e+6 -  2.84954e+6*y**2 +  1.06352e+6*y**4 + 280428.0*y**6 & 
           & -  240373.0*y**8 - 106663.0*y**10 + 486.14*y**12 + 3432.0*y**14
      o8 = 483737.0 -  498334.0*y**2 -  217801.0*y**4 - 161461.0*y**6 & 
           & -  55582.0*y**8 +  1093.82*y**10 + 3003.0*y**12
      p8 =  70946.1 -  55600.0*y**2 -  48153.3*y**4 -  16493.7*y**6 & 
           & +  1012.79*y**8 + 2002.0*y**10
      q8 = 9504.65 - 3058.26*y**2 - 1500.17*y**4 + 567.164*y**6 + 1001.0*y**8
      r8 = 955.194 +  533.254*y**2 +  198.876*y**4 + 364.0*y**6
      s8 =  126.532 + 40.5117*y**2 +  91.0*y**4
      t8 =  3.68288 +  14.0*y**2  
      out = dexp(y**2-x**2) * cos(2*x*y)  &
           & - (a7 + b7*x**2 + c7*x**4 + d7*x**6 + e7*x**8 + f7*x**10 + g7*x**12 + h7*x**14  &
           &     + o7*x**16 + p7*x**18 + q7*x**20 + r7*x**22 + s7*x**24 + t7*x**26 )         &
           &  /(a8 + b8*x**2 + c8*x**4 + d8*x**6 + e8*x**8 + f8*x**10 + g8*x**12 + h8*x**14  &
           &     + o8*x**16 + p8*x**18 + q8*x**20 + r8*x**22 + s8*x**24 + t8*x**26 + x**28 )           
   endif
 end function voigt
