  module lineshape_module

    use vtype_module

    implicit none

    private

    public :: lineshape_voigt, lineshape_voigt_old, lineshape_lorentz, &
         lineshape_renorm

    !**************************************************************************

  contains

    !**************************************************************************

    !**************************************************************************
    ! function lineshape_voigt
    ! calculate line shape factor of voigt type
    !**************************************************************************
    ! lorwid        : lorentz halfwidth (m^-1)
    ! dopwid        : doppler halfwidth (m^-1)
    ! dv            : dv=(v-v0) : v  : wave number
    !               :           : v0 : wave number at line center
    !               : this value will be casted absolute value
    !**************************************************************************

    !**************************************************************************
    function lineshape_voigt( lorwid, dopwid, dv ) result( voigt )
      !************************************************************************
      !****
      !****  routine computes the voigt function:                     ***
      !****     y/pi*integral from - to + infinity of                 ***
      !                              exp(-t*t)/(y*y+(x-t)*(x-t)) dt   ***
      !****
      !****  Return value seems to be divided by dopwid * sqrt(pi).
      !
      ! This is written in Liou, Radiation and Cloud Processes in the 
      ! Atmosphere, p31. 

      use const_module, only : pi

      real(dp), intent(in) :: lorwid, dopwid, dv

      real(dp)             :: voigt


      !
      ! local variables
      !
      real(dp)            :: ln2, c1, voiwid, xi, eta, etasq, pivoiwid


      ln2      = log( 2.0d0 )
      c1       = sqrt( lorwid * lorwid + 4.0d0 * ln2 * dopwid * dopwid )
      voiwid   = 0.5d0 * ( lorwid + c1 ) &
           + 0.05d0 * lorwid * ( 1.0d0 - ( lorwid + lorwid ) / ( lorwid + c1 ) )
      xi       = lorwid / voiwid
      eta      = dv / voiwid
      etasq    = eta * eta
      pivoiwid = pi * voiwid

      voigt = sqrt( ln2 / pi ) / voiwid &
           * ( 1.0d0 - xi ) * exp( -ln2 * etasq ) &
           + xi / pivoiwid &
           * ( 1.0d0 / ( 1.0d0 + etasq ) &
           - ( 1.0d0 - xi ) * ( 1.5d0 / ln2 + 1.0d0 + xi ) &
           * ( 0.066d0 * exp( -0.4d0 * etasq ) &
           - 1.0d0 / ( 40.0d0 - 5.5d0 * etasq + etasq * etasq ) ) )


    end function lineshape_voigt

    !**************************************************************************
    ! real*8 function voigt
    ! calculate line shape factor of voigt type
    !**************************************************************************
    ! lorwid        : lorentz halfwidth (m^-1)
    ! dopwid        : doppler halfwidth (m^-1)
    ! dv            : dv=(v-v0) : v  : wave number
    !               :           : v0 : wave number at line center
    !               : this value will be casted absolute value
    !**************************************************************************

    !**************************************************************************
    FUNCTION lineshape_VOIGT_old( lorwid, dopwid, dv ) result( voigt )
      !************************************************************************
      !****
      !****  ROUTINE COMPUTES THE VOIGT FUNCTION: Y/PI*INTEGRAL FROM  ***
      !****   - TO + INFINITY OF EXP(-T*T)/(Y*Y+(X-T)*(X-T)) DT       ***
      !****

      !****  Return value seems to be divided by dopwid * sqrt(pi).

      use const_module, only : pi

      real(dp), intent(in) :: lorwid, dopwid, dv

      real(dp)             :: voigt


      !
      ! local variables
      !
      REAL(dp) :: B(22),RI(15),XN(15),YN(15)
      REAL(dp) :: D0(25),D1(25),D2(25),D3(25),D4(25)
      REAL(dp) :: HN(25),H,XX(3),HH(3),NBY2(19),C(21)
      LOGICAL  :: TRU
      DATA B(1)/0.0/,B(2)/0.7093682E-7/
      DATA XN/10.0,9.0,2*8.0,7.0,6.0,5.0,4.0,7*3.0/
      DATA YN/3*0.6,0.5,2*0.4,4*0.3,1.0,0.9,0.8,2*0.7/
      DATA H/0.201/
      DATA XX/0.5246476,1.65068,0.7071068/
      DATA HH/0.2562121,0.2588268E-1,0.2820948/
      DATA NBY2/9.5,9.0,8.5,8.0,7.5,7.0,6.5,6.0,5.5,5.0,4.5,4.0,3.5,3.0, &
           2.5,2.0,1.5,1.0,0.5/
      DATA C / 0.7093682E-7,-0.2518434E-6, 0.8566874E-6,-0.2787638E-5, &
               0.8668740E-5,-0.2565551E-4, 0.7228775E-4,-0.1933631E-3, &
               0.4899520E-3,-0.1173267E-2, 0.2648762E-2,-0.5623190E-2, &
               0.1119681E-1,-0.2084976E-1, 0.3621573E-1,-0.5851412E-1, &
               0.8770816E-1,-0.121664, 0.15584, -0.184, 0.2           /
      DATA TRU/.FALSE./

      real(dp)     :: x, y
      real(dp)     :: co, uu, vv, u, y2, d, dx, v
      integer(i4b) :: i, j, max, min, n

      x = abs( dv ) / dopwid
      y = lorwid / dopwid

      IF (TRU) GO TO 104
      !****
      !****  REGION I: COMPUTE DAWSON'S FUNCTION AT MESH POINTS. ***
      !****
      TRU=.TRUE.
      DO I=1,15
        RI(I)=-I/2.0
     END DO
      DO I=1,25
        HN(I)=H*(I-0.5)
        CO=4*HN(I)*HN(I)/25-2
        DO J=2,21
          B(J+1)=CO*B(J)-B(J-1)+C(J)
       END DO
        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)
        D4(I)=(HN(I)*D3(I)+D2(I))/RI(4)
     END DO
  104 CONTINUE
      IF (X.GE.5.0) GO TO 112
      IF (Y.LE.1.0) GO TO 110
      IF (X.GT.1.85*(3.6-Y)) GO TO 112
      !****
      !****  REGION II: CONTINUED FRACTION. COMPUTE NUMBER OF TERMS NEEDED. ***
      !****
      IF (Y.LT.1.45) GO TO 107
      I=Y+Y
      GO TO 108
  107 CONTINUE
      I=11*Y
  108 CONTINUE
      J=X+X+1.85
      MAX=XN(J)*YN(I)+0.46
      MIN=MIN0(16,21-2*MAX)
      !****
      !****  EVALUATE CONTINUED FRACTION. ***
      !****
      UU=Y
      VV=X
      DO J=MIN,19
        U=NBY2(J)/(UU*UU+VV*VV)
        UU=Y+U*UU
        VV=X-U*VV
     END DO
      VOIGT=UU/(UU*UU+VV*VV)/1.772454
      voigt=1.0e0/(dopwid*sqrt(pi))*voigt
      RETURN
      !****
  110 CONTINUE
      Y2=Y*Y
      IF (X+Y.GE.5.0) GO TO 113
      !****
      !****  REGION I: COMPUTE DAWSON'S FUNCTION AT X FROM TAYLOR SERIES. ***
      !****
      N=X/H
      D=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
      !****
      !****  TAYLOR SERIES EXPANSION ABOUT Y=0.0. ***
      !****
      VV=EXP(Y2-X*X)*COS(2*X*Y)/1.128379-Y*V
      UU=-Y
      MAX=5+(12.5-X)*0.8*Y
      DO I=2,MAX,2
        U=(X*V+U)/RI(I)
        V=(X*U+V)/RI(I+1)
        UU=-UU*Y2
        VV=VV+V*UU
     END DO
      VOIGT=1.128379*VV
      voigt=1.0e0/(dopwid*sqrt(pi))*voigt
      RETURN
      !****
  112 CONTINUE
      Y2=Y*Y
      IF (Y.LT.11-0.6875*X) GO TO 113
      !****
      !****  REGION IIIB: 2-POINT GAUSS-HERMITE QUADRATURE. ***
      !****
      U=X-XX(3)
      V=X+XX(3)
      VOIGT=Y*(HH(3)/(Y2+U*U)+HH(3)/(Y2+V*V))
      voigt=1.0e0/(dopwid*sqrt(pi))*voigt
      RETURN
      !****
      !****  REGION IIIA: 4-POINT GAUSS-HERMITE QUADRATURE. ***
      !****
  113 CONTINUE
      U=X-XX(1)
      V=X+XX(1)
      UU=X-XX(2)
      VV=X+XX(2)
      VOIGT=Y*(HH(1)/(Y2+U*U)+HH(1)/(Y2+V*V) &
           +HH(2)/(Y2+UU*UU)+HH(2)/(Y2+VV*VV))
      voigt=1.0e0/(dopwid*sqrt(pi))*voigt
      RETURN
    END FUNCTION LINESHAPE_VOIGT_OLD

    !**************************************************************************
    function lineshape_lorentz( lorwid, dv ) result( lorentz )
      !************************************************************************
      !****
      !****  routine computes the Lorentz function: 
      !****           1/pi*lorwid/(lorwid**2+dv**2)
      !****

      use const_module, only : pi

      real(dp), intent(in) :: lorwid, dv

      real(dp)             :: lorentz


      !
      ! local variables
      !

      lorentz = lorwid / ( pi * ( lorwid * lorwid + dv * dv ) )


    end function lineshape_lorentz

    !**************************************************************************
    ! real*8 function renorm
    ! Renormalize Line Shape
    !**************************************************************************

    function lineshape_renorm( lowid, dowid, wnmax ) result( renorm )

      use const_module

      real(dp), intent(in ) :: lowid,dowid
      real(dp), intent(in ) :: wnmax

      real(dp)              :: renorm


      !
      ! local variables
      !
      integer(i4b) :: l, lm
      real(dp)     :: dwn, wn
      real(dp)     :: val1, val2


      dwn    = dowid / 5.0d0
      lm     = int( wnmax / dwn ) + 1
      dwn    = wnmax / dble( lm )

      renorm = 0.0d0

      l = 1
      wn   = dwn * dble( l - 1 )
      val1 = lineshape_voigt( lowid, dowid, wn )
      val2 = val1
      do l = 1+1, lm
         wn     = dwn * dble( l - 1 )
         val1   = lineshape_voigt( lowid, dowid, wn )
         renorm = renorm + ( val1 + val2 ) * dwn * 0.5d0
         val2   = val1


!!$         val1   = lineshape_voigt( lowid, dowid, wn )
!!$         val2   = lineshape_voigt2( lowid, dowid, wn )
!!$         renorm = lineshape_lorentz( lowid, wn )
!!$         write( 6, * ) wn, val1, val2, renorm
!!$         if( l .eq. lm ) stop


      end do
      renorm = renorm + renorm

      renorm = 1.0d0 / renorm


    end function lineshape_renorm

    !**************************************************************************

  end module lineshape_module
