MODULE make_model
        USE nrtype
        CONTAINS
!#########################
        SUBROUTINE make_model_TGM_quick(model,params)
        USE share, ONLY: w_sp,wave_ll,select_ll_mask&
        &,damp_coeff,dim_ll,coeff_4deg_quick,dimsp,disp,ew
        USE func_poly, ONLY: ew_poly4_quick, voigt_logg
        IMPLICIT NONE
        
        REAL(DP) :: sig,c,Dc,dl,gamg,gaml,ew_dp
        INTEGER(I2B) :: i,j,pix_lo,pix_up,w_center_line,Dsig
        REAL(DP),DIMENSION(dimsp), INTENT(INOUT) :: model
        REAL(DP),DIMENSION(:) :: params
        REAL(DP),DIMENSION(4) :: st_pars
        !initialize model
        model=1._dp
        Dc=10.

        IF(params(1)<3600) params(1)=3600._dp
        IF(params(1)>7400) params(1)=7400._dp
        IF(params(2)<0.1) params(2)=0.1_dp
        IF(params(2)>5.5) params(2)=5.5_dp
        IF(params(3)<-2.6) params(3)=-2.6_dp
        IF(params(3)>0.5) params(3)=0.5_dp
        IF(ABS(params(4))>5.) params(4)=5.0_dp

        st_pars(1:3)=params(1:3)
        st_pars(4)=0.0001_dp        
        sig=params(4)
        c=1.0_dp+params(5)/299792.0_dp


        DO i=1,dim_ll

          IF(select_ll_mask(i)) THEN
            ew(i)=ew_poly4_quick(st_pars,coeff_4deg_quick(:,i))
!        write(*,*) 'TGM',ele_ll(i),ew(i)
            w_center_line=INT(MINLOC(ABS(w_sp-wave_ll(i)),1),I2B)
            ew_dp=ew(i)*damp_coeff(i)
            !compute the parameters used by the voigt function
            dl=1._dp+(5._dp-st_pars(2)*2_dp)*0.1_dp
!            gaml=0.8_dp*(ew_dp)*(1._dp-exp(-(ew_dp/(0.2_dp*dl))**2))
            gaml=gammaL(st_pars(2),ew_dp)
            gamg=2.0_dp*sig*1.17741002252_dp!### sqrt(2.*log(2.))=1.17741002252
            Dc=MAX(6.,gaml*6/sig)
            Dsig=INT(Dc*sig/disp(i),I2B)
            pix_up=w_center_line+Dsig
            pix_lo=w_center_line-Dsig
            pix_lo=MAX(1_I2B,pix_lo)
            pix_up=MIN(dimsp,pix_up)
            DO j=pix_lo,pix_up
             model(j)=model(j)-voigt_logg(w_sp(j),wave_ll(i)*c,gamg,ew(i),gaml)
            END DO
          END IF
        END DO
        END SUBROUTINE make_model_TGM_quick
!#########################
        SUBROUTINE make_model_TGM(model,params)
        USE share, ONLY: ele_ll,ABD,w_sp,wave_ll,select_ll_mask&
        &,damp_coeff,dim_ll,coeff_4deg,dimsp,disp,ew,ele2meas
        USE func_poly, ONLY: ew_poly4, voigt_logg
        IMPLICIT NONE
        
        REAL(DP) :: sig,c,Dc,dl,gamg,gaml,ew_dp
        INTEGER(I2B) :: i,j,pix_lo,pix_up,w_center_line,ele_pos,Dsig
        REAL(DP),DIMENSION(dimsp), INTENT(INOUT) :: model
        REAL(DP),DIMENSION(:) :: params
        REAL(DP),DIMENSION(4) :: st_pars
        !initialize model
        model=1._dp
        Dc=10.

        IF(params(1)<3600) params(1)=3600._dp
        IF(params(1)>7400) params(1)=7400._dp
        IF(params(2)<0.1) params(2)=0.1_dp
        IF(params(2)>5.5) params(2)=5.5_dp
        IF(params(3)<-2.6) params(3)=-2.6_dp
        IF(params(3)>0.5) params(3)=0.5_dp
        IF(ABS(params(4))>5.) params(4)=5.0_dp

        st_pars(1:3)=params(1:3)
        st_pars(4)=0.0001_dp        
        sig=params(4)
        c=1.0_dp+params(5)/299792.0_dp

        DO i=1,dim_ll
          IF(select_ll_mask(i)) THEN
            ele_pos=INT(MINLOC(ABS(ele_ll(i)-ele2meas),1),I2B)
            st_pars(4)=ABD(ele_pos)
!        write(*,*) 'TGM',ele_ll(i),ele_pos,st_pars(4),ele2meas(ele_pos)
            ew(i)=ew_poly4(st_pars,coeff_4deg(:,i))
            w_center_line=INT(MINLOC(ABS(w_sp-wave_ll(i)),1),I2B)
            ew_dp=ew(i)*damp_coeff(i)
            !compute the parameters used by the voigt function
            dl=1._dp+(5._dp-st_pars(2)*2_dp)*0.1_dp
!            gaml=0.8_dp*(ew_dp)*(1._dp-exp(-(ew_dp/(0.2_dp*dl))**2))
            gaml=gammaL(st_pars(2),ew_dp)
            gamg=2.0_dp*sig*1.17741002252_dp!### sqrt(2.*log(2.))=1.17741002252
            Dc=MAX(6.,gaml*6/sig)
            Dsig=INT(Dc*sig/disp(i),I2B)
            pix_up=w_center_line+Dsig
            pix_lo=w_center_line-Dsig
            pix_lo=MAX(1_I2B,pix_lo)
            pix_up=MIN(dimsp,pix_up)
            DO j=pix_lo,pix_up
             model(j)=model(j)-voigt_logg(w_sp(j),wave_ll(i)*c,gamg,ew(i),gaml)
            END DO
          END IF
        END DO
        END SUBROUTINE make_model_TGM
!#########################
        SUBROUTINE make_model_ABD(model,params)
        USE share, ONLY: ele_ll,TGM,w_sp,wave_ll,select_ll_mask&
        &,damp_coeff,dim_ll,coeff_4deg,dimsp,disp,ew,ele2meas
        USE func_poly, ONLY: ew_poly4, voigt_logg
        IMPLICIT NONE
        
        REAL(DP) :: sig,c,Dc,dl,gamg,gaml,ew_dp
        INTEGER(I2B) :: i,j,pix_lo,pix_up,w_center_line,ele_pos,Dsig
        REAL(DP),DIMENSION(dimsp), INTENT(INOUT) :: model
        REAL(DP),DIMENSION(:) :: params
        REAL(DP),DIMENSION(4) :: st_pars
        !initialize model
        model=1._dp
        Dc=10.

        st_pars(1:3)=TGM(1:3)
        st_pars(4)=0.0001_dp        
        sig=TGM(4)
        c=1.0_dp+TGM(5)/299792.0_dp

        WHERE(params>0.8) params=0.8_dp
        WHERE(params<-0.6) params=-0.6_dp

        DO i=1,dim_ll
          IF(select_ll_mask(i)) THEN
            ele_pos=INT(MINLOC(ABS(ele_ll(i)-ele2meas),1),I2B)
            st_pars(4)=params(ele_pos)
!        write(*,*) 'ABD',ele_ll(i),ele_pos(1),st_pars(4),ele2meas(ele_pos(1))
            ew(i)=ew_poly4(st_pars,coeff_4deg(:,i))
            w_center_line=INT(MINLOC(ABS(w_sp-wave_ll(i)),1),I2B)
            !compute the parameters used by the voigt function
            dl=1._dp+(5._dp-st_pars(2)*2_dp)*0.1_dp
            ew_dp=ew(i)*damp_coeff(i)
!            gaml=0.8_dp*(ew_dp)*(1._dp-exp(-(ew_dp/(0.2_dp*dl))**2))
            gaml=gammaL(st_pars(2),ew_dp)
            gamg=2.0_dp*sig*1.17741002252_dp!### sqrt(2.*log(2.))=1.17741002252
            Dc=MAX(6.,gaml*6/sig)
            Dsig=INT(Dc*sig/disp(i),I2B)
            pix_up=w_center_line+Dsig
            pix_lo=w_center_line-Dsig
            pix_lo=MAX(1_I2B,pix_lo)
            pix_up=MIN(dimsp,pix_up)
            DO j=pix_lo,pix_up
             model(j)=model(j)-voigt_logg(w_sp(j),wave_ll(i)*c,gamg,ew(i),gaml)
             END DO
          END IF
        END DO

        END SUBROUTINE make_model_ABD
!#########################
        SUBROUTINE make_model_ABDerr(model,parTGM,parABD)
        USE share, ONLY: ele_ll,w_sp,wave_ll,select_ll_mask&
        &,damp_coeff,dim_ll,coeff_4deg,dimsp,disp,ew,ele2meas
        USE func_poly, ONLY: ew_poly4, voigt_logg
        IMPLICIT NONE
        
        REAL(DP) :: sig,c,Dc,gamg,gaml,ew_dp,dl
        INTEGER(I2B) :: i,j,pix_lo,pix_up,w_center_line,Dsig
        INTEGER(I1B) :: ele_pos
        REAL(DP),DIMENSION(dimsp), INTENT(INOUT) :: model
        REAL(DP),DIMENSION(:) :: parTGM,parABD
        REAL(DP),DIMENSION(4) :: st_pars
        !initialize model
        model=1._dp
        Dc=10.

        st_pars(1:3)=parTGM(1:3)
        st_pars(4)=0.0001_dp        
        sig=parTGM(4)
        c=1._dp+parTGM(5)/299792._dp

        WHERE(parABD>0.8) parABD=0.8_dp
        WHERE(parABD<-0.6) parABD=-0.6_dp

        DO i=1,dim_ll
          IF(select_ll_mask(i)) THEN
            ele_pos=INT(MINLOC(ABS(ele_ll(i)-ele2meas),1),I1B)
            st_pars(4)=parABD(ele_pos)
            ew(i)=ew_poly4(st_pars,coeff_4deg(:,i))
            w_center_line=INT(MINLOC(ABS(w_sp-wave_ll(i)),1),I2B)
            ew_dp=ew(i)*damp_coeff(i)
            !compute the parameters used by the voigt function
            dl=1._dp+(5._dp-st_pars(2)*2_dp)*0.1_dp
            ew_dp=ew(i)*damp_coeff(i)
!            gaml=0.8_dp*(ew_dp)*(1._dp-exp(-(ew_dp/(0.2_dp*dl))**2))
            gaml=gammaL(st_pars(2),ew_dp)
            gamg=2.0_dp*sig*1.17741002252_dp!### sqrt(2.*log(2.))=1.17741002252
            Dc=MAX(6.,gaml*6/sig)
            Dsig=INT(Dc*sig/disp(i),I2B)
            pix_up=w_center_line+Dsig
            pix_lo=w_center_line-Dsig
            pix_lo=MAX(1_I2B,pix_lo)
            pix_up=MIN(dimsp,pix_up)
            DO j=pix_lo,pix_up
             model(j)=model(j)-voigt_logg(w_sp(j),wave_ll(i)*c,gamg,ew(i),gaml)
            END DO
          END IF
        END DO

        END SUBROUTINE make_model_ABDerr
!#########################
        REAL(DP) FUNCTION gammaL1(logg,ewdp)
        IMPLICIT NONE
        REAL(DP), INTENT(IN) :: logg, ewdp
        REAL(DP) :: dl
        dl=1._dp+(5._dp-logg*2_dp)*0.1_dp
        gammaL1=0.8_dp*(ewdp)*(1._dp-exp(-(ewdp/(0.2_dp*dl))**2))
        
        END FUNCTION gammaL1
!#########################
        REAL(DP) FUNCTION gammaL(logg,ewdp)
        IMPLICIT NONE
        REAL(DP), INTENT(IN) :: logg, ewdp

        IF(logg>4.5) THEN
         gammaL=0.8*(ewdp)*(1.-exp(-((ewdp)/(0.14))**2))
        ELSE IF(logg<=4.5.AND.logg>3.5) THEN
         gammaL=(0.7+(logg-3.5)*0.1)*(ewdp)*(1.-exp(-((ewdp)/(0.16))**2))
        ELSE IF (logg<=3.5.AND.logg>2.5) THEN
         gammaL=(0.6+(logg-2.5)*0.1)*(ewdp)*(1.-exp(-((ewdp)/(0.2))**2))  
        ELSE IF (logg<=2.5.AND.logg>1.5) THEN
         gammaL=0.6*(ewdp)*(1.-exp(-((ewdp)/(0.2))**2))  
        ELSE IF (logg<=1.5) THEN
         gammaL=(0.6+(1.5-logg)*0.1)*(ewdp)*(1.-exp(-((ewdp)/(0.2))**2))  
       END IF
        
        END FUNCTION gammaL
!########################
END MODULE make_model