!    This module contains a routine to compute the EWs of the lines from their
!    polynomial GCOG, and a function for the voigt profile.
!
!    It is part of the program SP_Ace, which derives stellar parameters, 
!    such as gravity, temperature, and element abundances from optical
!    stellar spectra, assuming Local Thermodynamic Equilibrium (LTE) 
!    and 1D stellar atmosphere models.
!
!    Copyright (C) 2016 Corrado Boeche
!    
!    This program is free software: you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation, either version 3 of the License, or
!    (at your option) any later version.

!
!    This program is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!    
!    You should have received a copy of the GNU General Public License
!    along with this program.  If not, see <http://www.gnu.org/licenses/>.


MODULE func_poly

        USE num_type
        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_lin_interp_a(pars,coeff)
!        USE data_lib, ONLY: temp_gridL,logg_gridL,met_gridL,aFe_gridL
!        USE read_GCOG, ONLY: load_GCOG_lin
!        IMPLICIT NONE
!        REAL(DP), DIMENSION(4), INTENT(IN) :: pars
!        REAL(DP), DIMENSION(48), INTENT(IN) :: coeff
!!        INTEGER(I1B) :: i1,i2,i3,i4
!        INTEGER(I1B) :: i,j,k,l,count,place
!        INTEGER(I1B) :: loc_a
!        REAL(DP) :: t_prox0,l_prox0,m_prox0,a_prox0
!        REAL(DP), SAVE :: t_prox_old,l_prox_old,m_prox_old
!        REAL(DP) :: t_step,l_step,m_step,a_step,xt,xl,xm,xa,a,b,c,d
!        LOGICAL :: flag_tm,flag_lm,flag_mm,flag_am,flag_tp,flag_lp,flag_mp,flag_ap
!!        REAL(DP) :: nt_steps,nl_steps,nm_steps,na_steps
!
!        t_prox0=par_prox(temp_gridL,pars(1),flag_tm,flag_tp)
!        l_prox0=par_prox(logg_gridL,pars(2),flag_lm,flag_lp)
!        m_prox0=par_prox(met_gridL,pars(3),flag_mm,flag_mp)
!        a_prox0=par_prox(aFe_gridL,pars(4),flag_am,flag_ap)
!        t_step=200_dp
!        l_step=0.4_dp
!        m_step=0.2_dp
!        a_step=0.2_dp
!        loc_a=INT(MINLOC(ABS(aFe_gridL-a_prox0),1),I1B)
!
!!   write(*,*) 'flags',pars,flag_tm,flag_tp,flag_lm,flag_lp,flag_mm,flag_mp,flag_am,flag_ap
!!   write(*,*) 'prox0',t_prox0,l_prox0,m_prox0,a_prox0,loc_a
!
!        IF(t_prox0/=t_prox_old.OR.l_prox0/=l_prox_old.OR.&
!        & m_prox0/=m_prox_old) THEN
!          CALL load_GCOG_lin(t_prox0,l_prox0,m_prox0)
!          t_prox_old=t_prox0
!          l_prox_old=l_prox0
!          m_prox_old=m_prox0
!        END IF
!
!
!          xt=(pars(1)-t_prox0)/t_step
!          xl=(pars(2)-l_prox0)/l_step
!          xm=(pars(3)-m_prox0)/m_step
!          xa=(pars(4)-a_prox0)/a_step
!
!!#### if pars<boundary define the number of steps between pars and the minimum boundary
!
!        ew_lin_interp_a=0._dp
!        DO l=0,1
!        count=-5_I1B
!        DO k=0,1
!        DO j=0,1
!        DO i=0,1
!          IF(i==0) a=1-xt
!          IF(i==1) a=xt
!          IF(j==0) b=1-xl
!          IF(j==1) b=xl
!          IF(k==0) c=1-xm
!          IF(k==1) c=xm
!          IF(l==0) d=1-xa
!          IF(l==1) d=xa
!          count=count+6_I1B
!          place=count+loc_a+l
!!        write(*,*) i,j,k,l,count,loc_a,l,place
!        ew_lin_interp_a=ew_lin_interp_a+coeff(place)*a*b*c*d
!
!        ENDDO
!        ENDDO
!        ENDDO
!        ENDDO
!
!!   IF(flag_ap.OR.flag_am) THEN
!!   IF(ABS(pars(4))>0.001) THEN
!!    write(*,*) 'pars',pars,t_prox0,l_prox0,m_prox0,a_prox0
!!    write(*,*) 'coeff',coeff
!!    write(*,*) 'ew_interp',ew_lin_interp
!!   END IF
!
!        ew_lin_interp_a=ew_lin_interp_a/1000._dp
!
!        END FUNCTION ew_lin_interp_a
!#######################################        
        REAL(DP) FUNCTION ew_lin_interp(pars,coeff)
        USE data_lib, ONLY: temp_gridL,logg_gridL,met_gridL,aFe_gridL
        USE read_GCOG, ONLY: load_GCOG_lin
        USE utils, ONLY : find_prox_L
        IMPLICIT NONE
        REAL(DP), DIMENSION(4), INTENT(IN) :: pars
        REAL(DP), DIMENSION(48), INTENT(IN) :: coeff
!        INTEGER(I1B) :: i1,i2,i3,i4
        INTEGER(I1B) :: i,j,k,l,count,place
        INTEGER(I1B) :: loc_a
        REAL(DP),DIMENSION(4) :: tgm_prox
        REAL(DP), DIMENSION(3), SAVE :: tgm_prox_old
        REAL(DP) :: t_step,l_step,m_step,a_step,xt,xl,xm,xa,a,b,c,d
        LOGICAL :: flag_tm,flag_lm,flag_mm,flag_am,flag_tp,flag_lp,flag_mp,flag_ap
        REAL(DP) :: nt_steps,nl_steps,nm_steps,na_steps,ew_to_add
        REAL(DP) :: diff_t,diff_l,diff_m,diff_a

        CALL find_prox_L(pars,tgm_prox,flag_tm,flag_tp,flag_lm,flag_lp,flag_mm,flag_mp,flag_am,flag_ap)
        t_step=temp_gridL(2)-temp_gridL(1)
        l_step=logg_gridL(2)-logg_gridL(1)
        m_step=met_gridL(2)-met_gridL(1)
        a_step=aFe_gridL(2)-aFe_gridL(1)
        loc_a=INT(MINLOC(ABS(aFe_gridL-tgm_prox(4)),1),I1B)


!   write(*,*) 'flags',pars,flag_tm,flag_tp,flag_lm,flag_lp,flag_mm,flag_mp,flag_am,flag_ap


        IF(ANY(tgm_prox(1:3)/=tgm_prox_old)) THEN
          CALL load_GCOG_lin(tgm_prox(1),tgm_prox(2),tgm_prox(3))
          tgm_prox_old=tgm_prox(1:3)
        END IF


        IF(flag_tp.AND..NOT.flag_tm) THEN
          xt=1._dp
        ELSE IF (flag_tm.AND..NOT.flag_tp) THEN 
          xt=0._dp
        ELSE IF(.NOT.flag_tm.AND..NOT.flag_tp) THEN
          xt=(pars(1)-tgm_prox(1))/t_step
        END IF

        IF(flag_lp.AND..NOT.flag_lm) THEN
          xl=1._dp
        ELSE IF(flag_lm.AND..NOT.flag_lp) THEN
          xl=0._dp
        ELSE IF(.NOT.flag_lm.AND..NOT.flag_lp) THEN
          xl=(pars(2)-tgm_prox(2))/l_step
        END IF

        IF(flag_mp.AND..NOT.flag_mm) THEN
          xm=1._dp
        ELSE IF(flag_mm.AND..NOT.flag_mp) THEN
          xm=0._dp
        ELSE IF(.NOT.flag_mm.AND..NOT.flag_mp) THEN
          xm=(pars(3)-tgm_prox(3))/m_step
        END IF

        IF(flag_ap.AND..NOT.flag_am) THEN
          xa=1._dp
!     write(*,*) 'flag_ap',pars(4),t_prox0,l_prox0,m_prox0,a_prox0,xa
        ELSE IF(flag_am.AND..NOT.flag_ap) THEN
          xa=0._dp
!     write(*,*) 'flag_am',pars(4),t_prox0,l_prox0,m_prox0,a_prox0,xa
        ELSE IF(.NOT.flag_am.AND..NOT.flag_ap) THEN
          xa=(pars(4)-tgm_prox(4))/a_step
        END IF

!   write(*,*) 'x',xt,xl,xm,xa
!#### if pars<boundary define the number of steps between pars and the minimum boundary
        nt_steps=0.
        nl_steps=0
        nm_steps=0
        na_steps=0
        IF(flag_tm) nt_steps=(pars(1)-temp_gridL(1))/t_step
        IF(flag_lp) nl_steps=(pars(2)-logg_gridL(INT(SIZE(logg_gridL),I2B)))/l_step
        IF(flag_lm) nl_steps=(pars(2)-logg_gridL(1))/l_step
        IF(flag_mp) nm_steps=(pars(3)-met_gridL(INT(SIZE(met_gridL),I2B)))/m_step
        IF(flag_mm) nm_steps=(pars(3)-met_gridL(1))/m_step
        IF(flag_ap) na_steps=(pars(4)-aFe_gridL(INT(SIZE(aFe_gridL),I2B)))/a_step
        IF(flag_am) na_steps=(pars(4)-aFe_gridL(1))/a_step
! write(*,*) nt_steps,nl_steps,nm_steps,na_steps

        ew_lin_interp=0._dp
        DO l=0,1
        count=-6_I1B
        DO k=0,1
        DO j=0,1
        DO i=0,1
          IF(i==0) a=1-xt
          IF(i==1) a=xt
          IF(j==0) b=1-xl
          IF(j==1) b=xl
          IF(k==0) c=1-xm
          IF(k==1) c=xm
          IF(l==0) d=1-xa
          IF(l==1) d=xa
          count=count+6_I1B
          place=count+loc_a+l
!        write(*,*) i,j,k,l,count,loc_a,l,place,coeff(place)
!        write(*,*) a,b,c,d
        ew_lin_interp=ew_lin_interp+coeff(place)*a*b*c*d

        ENDDO
        ENDDO
        ENDDO
        ENDDO

        diff_t=(coeff(3+6)-coeff(3))/t_step
        diff_l=(coeff(3+12)-coeff(3))/l_step
        diff_m=(coeff(3+24)-coeff(3))/m_step
        diff_a=coeff(4)-coeff(3)

        ew_to_add=0.
        IF(flag_tp) ew_to_add=ew_to_add+diff_t*nt_steps
        IF(flag_lp) ew_to_add=ew_to_add+diff_t*nl_steps
        IF(flag_mp) ew_to_add=ew_to_add+diff_l*nm_steps
        IF(flag_ap) ew_to_add=ew_to_add+diff_l*na_steps
        IF(flag_tm) ew_to_add=ew_to_add+diff_m*nt_steps
        IF(flag_lm) ew_to_add=ew_to_add+diff_m*nl_steps
        IF(flag_mm) ew_to_add=ew_to_add+diff_a*nm_steps
        IF(flag_am) ew_to_add=ew_to_add+diff_a*na_steps

!   IF(flag_ap.OR.flag_am) THEN
!    write(*,*) 'ew_interp',ew_lin_interp
!   END IF
        ew_lin_interp=ew_lin_interp+MAX(0.,ew_to_add)
        ew_lin_interp=ew_lin_interp/1000._dp

        END FUNCTION ew_lin_interp
!#####################
!        REAL(DP) FUNCTION par_prox(array,par,flagm,flagp)
!        IMPLICIT NONE
!        REAL(DP), DIMENSION(:), INTENT(IN) :: array
!        REAL(DP), INTENT(IN) :: par
!        LOGICAL, INTENT(OUT) :: flagm,flagp
!        INTEGER(I1B) :: loc
!
!        flagm=.FALSE.
!        flagp=.FALSE.
!!### note: array(1:SIZE(array)-1) because at the last point the GCOG library has no results 
!        loc=INT(MINLOC(ABS(array(1:SIZE(array)-1)-par),1),I1B)
!
!        IF(par-array(loc)>=0) THEN
!          par_prox=array(loc)
!          IF(par>array(SIZE(array))) flagp=.TRUE.
!        ELSE
!         IF(par>array(2)) THEN
!          par_prox=array(loc-1_I1B)        
!         ELSE IF(par<array(2).AND.par>=array(1)) THEN
!          par_prox=array(1)
!         ELSE
!          par_prox=array(1)
!          flagm=.TRUE.
!         END IF
!        END IF
!
!        END FUNCTION par_prox
!#####################
        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,gamG,ew,gamL)
        USE num_type
        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
        REAL(DP), PARAMETER :: pi=3.14159265358979_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
!#####################

END MODULE func_poly