MODULE func_poly

        USE nrtype
        CONTAINS
        
        REAL(DP) FUNCTION ew_poly4(pars,coeff)
        IMPLICIT NONE
        REAL(DP), DIMENSION(4), INTENT(IN) :: pars
        REAL(DP), DIMENSION(70), INTENT(IN) :: coeff
        INTEGER(I1B), PARAMETER :: grado=4
        INTEGER(I1B) :: i1,i2,i3,i4
        INTEGER(I1B) :: count
        REAL(DP) :: a,b,c

        !initialize variables
        ew_poly4=0.0_dp
        count=0

        DO i1=0,grado
         a=(pars(1)**i1)
         DO i2=0,grado-i1
          b=a*(pars(2)**i2)
          DO i3=0,grado-i2-i1
           c=b*(pars(3)**i3)
           DO i4=0,grado-i3-i2-i1
             count=count+1_I1B
             ew_poly4=ew_poly4+coeff(count)*c*(pars(4)**i4)
           END DO
          END DO
         END DO
        END DO                                 

        IF(ew_poly4<1e-6) ew_poly4=1e-6_dp

        END FUNCTION ew_poly4
!#####################
        REAL(DP) FUNCTION ew_poly4_quick(pars,coeff)
        IMPLICIT NONE
        REAL(DP), DIMENSION(3), INTENT(IN) :: pars
        REAL(DP), DIMENSION(84), INTENT(IN) :: coeff
        INTEGER(I1B), PARAMETER :: grado=6
        INTEGER(I1B) :: i1,i2,i3
        INTEGER(I1B) :: count
        REAL(DP) :: a,b

        !initialize variables
        ew_poly4_quick=0.0_dp
        count=0

        DO i1=0,grado
         a=(pars(1)**i1)
         DO i2=0,grado-i1
          b=a*(pars(2)**i2)
          DO i3=0,grado-i2-i1
             count=count+1_I1B
             ew_poly4_quick=ew_poly4_quick+coeff(count)*b*(pars(3)**i3)
          END DO
         END DO
        END DO                                 

        IF(ew_poly4_quick<1e-6) ew_poly4_quick=1e-6_dp

        END FUNCTION ew_poly4_quick
!#####################
!        REAL(DP) FUNCTION voigt_logg(w,mu,sigma,ew,dd)
!        USE nrtype
!        IMPLICIT NONE
!        REAL(DP),INTENT(IN) :: w,mu,sigma,ew,dd
!        REAL(DP),PARAMETER :: sqrtln2=0.832554611_dp,sqrtpi=1.772453851_dp
!        REAL(DP), DIMENSION(4), PARAMETER :: A=(/-1.2150,-1.3509,-1.2150,-1.3509/)
!        REAL(DP), DIMENSION(4), PARAMETER :: B=(/1.2359,0.3786,-1.2359,-0.3786/)
!        REAL(DP), DIMENSION(4), PARAMETER :: C=(/-0.3085,0.5906,-0.3085,0.5906/)
!        REAL(DP), DIMENSION(4), PARAMETER :: D=(/0.0210,-1.1858,-0.0210,1.1858/)
!        REAL(DP), DIMENSION(4) :: V
!        REAL(DP) :: gamL,gamG,sigmaL,aL,X,Y
!
!
!        gamG=2.0_dp*sigma*1.17741002252_dp!### sqrt(2.*log(2.))=1.17741002252
!
!        gamL=1.0_dp-exp(-(ew/dd)**1.5_dp)
!
!        sigmaL=(gamL*0.5_dp)/sqrtpi
!        aL=ew/(sigmaL*pi*sqrtpi)
!        X=(w-mu)*2.0_dp*sqrtln2/gamG
!        Y=gamL*sqrtln2/gamG
!
!        V=(C*(Y-A)+D*(X-B))/((Y-A)**2+(X-B)**2)
!        voigt_logg=SUM(V)*(gamL*aL*sqrtpi*sqrtln2/gamG)
!        voigt_logg=MAX(voigt_logg,1e-6_dp)
!        
!        END FUNCTION voigt_logg
!#####################
!        REAL(DP) FUNCTION voigt_logg1(w,w1,mu,sigma,ew,dd)
!        USE nrtype
!        IMPLICIT NONE
!        REAL(DP),INTENT(IN) :: w,w1,mu,sigma,ew,dd
!        REAL(DP),PARAMETER :: sqrtln2=0.832554611_dp,sqrtpi=1.772453851_dp
!        REAL(DP), DIMENSION(4), PARAMETER :: A=(/-1.2150,-1.3509,-1.2150,-1.3509/)
!        REAL(DP), DIMENSION(4), PARAMETER :: B=(/1.2359,0.3786,-1.2359,-0.3786/)
!        REAL(DP), DIMENSION(4), PARAMETER :: C=(/-0.3085,0.5906,-0.3085,0.5906/)
!        REAL(DP), DIMENSION(4), PARAMETER :: D=(/0.0210,-1.1858,-0.0210,1.1858/)
!        REAL(DP), DIMENSION(4) :: V
!        REAL(DP) :: gamL,gamG,sigmaL,aL,X,Y,del,w_lo,w_med,vv
!        INTEGER(I1B) :: i
!
!        gamG=2.0_dp*sigma*1.17741002252_dp!### sqrt(2.*log(2.))=1.17741002252
!        gamL=1.0_dp-exp(-(ew/dd)**1.5_dp)
!
!        sigmaL=(gamL*0.5_dp)/sqrtpi
!        aL=ew/(sigmaL*pi*sqrtpi)
!        Y=gamL*sqrtln2/gamG
!
!        vv=0.
!        
!        del=(w1-w)/6._dp
!        w_lo=w-(w1-w)*0.5_dp
!        DO i=2,6,2
!        w_med=w_lo+del*(i-1)
!        X=(w_med-mu)*2.0_dp*sqrtln2/gamG
!        V=(C*(Y-A)+D*(X-B))/((Y-A)**2+(X-B)**2)
!        vv=vv+(SUM(V)*(gamL*aL*sqrtpi*sqrtln2/gamG))
!        END DO
!        voigt_logg1=MAX(vv/3._dp,1e-6_dp)
!
!        END FUNCTION voigt_logg1
!############################
        REAL(DP) FUNCTION voigt_logg(w,mu,gamG,ew,gamL)
        USE nrtype
        IMPLICIT NONE
        REAL(DP),INTENT(IN) :: w,mu,gamG,ew,gamL
        REAL(DP),PARAMETER :: sqrtln2=0.832554611_dp,sqrtpi=1.772453851_dp
        REAL(DP), DIMENSION(4), PARAMETER :: A=(/-1.2150,-1.3509,-1.2150,-1.3509/)
        REAL(DP), DIMENSION(4), PARAMETER :: B=(/1.2359,0.3786,-1.2359,-0.3786/)
        REAL(DP), DIMENSION(4), PARAMETER :: C=(/-0.3085,0.5906,-0.3085,0.5906/)
        REAL(DP), DIMENSION(4), PARAMETER :: D=(/0.0210,-1.1858,-0.0210,1.1858/)
        REAL(DP), DIMENSION(4) :: V
        REAL(DP) :: sigmaL,aL,X,Y


        sigmaL=(gamL*0.5_dp)/sqrtpi
        aL=ew/(sigmaL*pi*sqrtpi)
        X=(w-mu)*2.0_dp*sqrtln2/gamG
        Y=gamL*sqrtln2/gamG

        V=(C*(Y-A)+D*(X-B))/((Y-A)**2+(X-B)**2)
        voigt_logg=SUM(V)*(gamL*aL*sqrtpi*sqrtln2/gamG)
        voigt_logg=MAX(voigt_logg,1e-6_dp)
        
        END FUNCTION voigt_logg
!#####################

END MODULE func_poly