!    
!    This module contains subroutines that estimate uncertains
!    in the stellar parameters and chemical abundances.
!
!    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 uncertains2
        USE num_type
        USE share, ONLY: up_ABD,lo_ABD
        USE utils, ONLY: ABD_mask,write_ABD_mask
        IMPLICIT NONE
        INTEGER (I1B) :: pos
        REAL (DP), DIMENSION(3) :: up_TGM,lo_TGM
        REAL (DP), DIMENSION(3), PARAMETER :: prob=(/3.53,6.25,8.02/)!for 3 deg of freedom, correspond to %68.3,%90,%95.4
        REAL(DP),DIMENSION(5) :: temp_gridS,logg_gridS,met_gridS
        LOGICAL :: flag_lim,flag_limS,flag_move,flag_TGM

        CONTAINS
!############################
       REAL(DP) FUNCTION chisq_tgm_e(x)
       USE make_model, ONLY: make_model_TGM
       USE share, ONLY: dimsp,TGM,f_sp_norm,sig_noise,TGM_prox
       USE utils, ONLY: find_proxS
       IMPLICIT NONE
       REAL (DP), DIMENSION(3), INTENT(IN) :: x 
       REAL (DP), DIMENSION(dimsp) :: model
       REAL (DP), DIMENSION(5) :: xTGM 

       xTGM(1:3)=x(1:3)
       xTGM(4:5)=TGM(4:5)
       !find the closest point in the sub-grid and verify if TGM is inside the grid
       CALL find_proxS(xTGM,temp_gridS,logg_gridS,met_gridS,TGM_prox,&
       &flag_lim,flag_limS,flag_move)

       IF(flag_limS) THEN
        CALL update_grid(xTGM)
       END IF

       CALL make_model_TGM(model,xTGM)
       chisq_tgm_e=SUM(((f_sp_norm-model)/sig_noise)**2)

       END FUNCTION chisq_tgm_e
!############################
       REAL(DP) FUNCTION chisq_abd_e(x)
       USE make_model, ONLY: make_model_ABDerr
       USE share, ONLY: dimsp,dim_ele,TGM,ABD,f_sp_norm,sig_noise
       USE share, ONLY: dimsp,TGM,f_sp_norm,sig_noise,TGM_prox
       USE utils, ONLY: find_proxS
       IMPLICIT NONE
       REAL (DP), DIMENSION(3), INTENT(IN) :: x 
       REAL (DP), DIMENSION(dimsp) :: model
       REAL (DP), DIMENSION(5) :: parTGM 
       REAL (DP), DIMENSION(dim_ele) :: parABD


       parTGM(1:2)=x(1:2)
       parTGM(3:5)=TGM(3:5)

       !find the closest point in the sub-grid and verify if TGM is inside the grid
       CALL find_proxS(parTGM,temp_gridS,logg_gridS,met_gridS,TGM_prox,&
       &flag_lim,flag_limS,flag_move)
       IF(flag_limS) THEN
        CALL update_grid(parTGM)
       END IF

       parABD=ABD
       parABD(pos)=x(3)

       CALL make_model_ABDerr(model,parTGM,parABD)
       chisq_abd_e=SUM(((f_sp_norm-model)/sig_noise)**2)

       END FUNCTION chisq_abd_e
!############################
SUBROUTINE update_grid(TGM_in)
       USE read_GCOG, ONLY: load_GCOG_lin
       USE share, ONLY: TGM_prox,sigma,coeff_ew_lin
       USE func_poly, ONLY: ew_lin_interp
       USE utils, ONLY: find_prox,find_proxS,make_gridS,select_lines
       IMPLICIT NONE
       REAL (DP), DIMENSION(:), INTENT(IN) :: TGM_in
       REAL(DP), DIMENSION(4) :: pars

        !make a new subgrids
        CALL find_prox(TGM_in,TGM_prox)
        CALL make_gridS(TGM_prox,temp_gridS,logg_gridS,met_gridS)
       !find the closest point in the sub-grid and verify if TGM is inside the grid
       CALL find_proxS(TGM_in,temp_gridS,logg_gridS,met_gridS,TGM_prox,&
       &flag_lim,flag_limS,flag_move)

       IF(flag_limS) THEN
        ! allocate and upload the table coeff_4deg which contains the coefficient
        ! to construct the polynomial GCOG
        CALL load_GCOG_lin(TGM_prox(1),TGM_prox(2),TGM_prox(3))
        pars=(/TGM_in(1),TGM_in(2),TGM_in(3),0._dp/)
        CALL select_lines(ew_lin_interp,coeff_ew_lin,pars,sigma,.FALSE.)
       END IF

END SUBROUTINE update_grid
!############################
       SUBROUTINE TGM_e(chisq_best)
       USE share, ONLY: TGM
       USE data_lib, ONLY: temp_gridL,logg_gridL,met_gridL
       IMPLICIT NONE
       REAL(DP), INTENT(IN) :: chisq_best
       REAL(DP), DIMENSION(3) :: TGM_best,TGM_test,step,max_step
       REAL(DP), DIMENSION(9,3) :: grid
       REAL(DP) :: chisq
       INTEGER (I2B) :: i,j
       REAL(DP) :: chisq_lim,zoom_par
       LOGICAL :: flag,flag_one



        chisq_lim=chisq_best+prob(1)
        CALL make_grid(pos,grid)
        !estimate the upper Teff, logg, [M/H] uncertains
           zoom_par=1.
           max_step=(/100.,0.2,0.1/);
           TGM_best=TGM(1:3)
           flag_one=.FALSE.
           DO j=1,200
            IF(zoom_par<0.05) EXIT
            flag=.FALSE.
           DO i=1,9
            step=grid(i,:)*max_step*zoom_par
            TGM_test=TGM_best(1:3)+step
            !do it only for the points inside the parameter grid
            IF(TGM_test(1)<MAXVAL(temp_gridL).AND.TGM_test(1)>MINVAL(temp_gridL)&
            &.AND.TGM_test(2)<MAXVAL(logg_gridL).AND.TGM_test(2)>MINVAL(logg_gridL)&
            &.AND.TGM_test(3)<MAXVAL(met_gridL).AND.TGM_test(3)>MINVAL(met_gridL)) THEN
             chisq=chisq_tgm_e(TGM_test)
             IF(chisq<chisq_lim.AND.TGM_test(pos)>TGM_best(pos)) THEN
              TGM_best=TGM_test
              flag=.TRUE.
              flag_one=.TRUE.
              EXIT
             ELSE
             END IF
            END IF
           END DO

           IF(.not.flag) zoom_par=zoom_par*0.5
           END DO

          IF(.NOT.flag_one) THEN
           up_TGM(pos)=100000
          ELSE
           up_TGM(pos)=TGM_best(pos)
          END IF
          
        CALL update_grid(TGM)

        !estimate the lower Teff, logg, [M/H] uncertains
           zoom_par=1. 
           TGM_best=TGM(1:3)
           flag_one=.FALSE. 
           DO j=1,200
            IF(zoom_par<0.01) EXIT
            flag=.FALSE.
           DO i=1,9
            step=grid(i,:)*max_step*zoom_par
            TGM_test=TGM_best(1:3)-step
            !do it only for the points inside the parameter grid
            IF(TGM_test(1)<MAXVAL(temp_gridL).AND.TGM_test(1)>MINVAL(temp_gridL)&
            &.AND.TGM_test(2)<MAXVAL(logg_gridL).AND.TGM_test(2)>MINVAL(logg_gridL)&
            &.AND.TGM_test(3)<MAXVAL(met_gridL).AND.TGM_test(3)>MINVAL(met_gridL)) THEN
             chisq=chisq_tgm_e(TGM_test)
             IF(chisq<chisq_lim.AND.TGM_test(pos)<TGM_best(pos)) THEN
              TGM_best=TGM_test
              flag=.TRUE.
              flag_one=.TRUE.
              EXIT
             ELSE
             END IF
            END IF
           END DO

           IF(.not.flag) zoom_par=zoom_par*0.5
           END DO


   CLOSE(9)

          IF(.NOT.flag_one) THEN
          lo_TGM(pos)=-100000
          ELSE
          lo_TGM(pos)=TGM_best(pos)
          END IF

       END SUBROUTINE TGM_e
!############################
       SUBROUTINE ABD_e(chisq_best)
       USE share, ONLY: TGM,ABD
       USE data_lib, ONLY: temp_gridL,logg_gridL,met_gridL
       IMPLICIT NONE
       REAL(DP), INTENT(IN) :: chisq_best
       REAL(DP), DIMENSION(3) :: TGM_best,TGM_test,step,max_step
       REAL(DP), DIMENSION(9,3) :: grid
       REAL(DP) :: chisq
       INTEGER (I2B) :: i,j
       REAL(DP) :: chisq_lim,zoom_par
       LOGICAL :: flag,flag_one

        CALL make_grid(INT(3,I1B),grid)
        chisq_lim=chisq_best+prob(1)

        !estimate the upper [El/m] uncertains
           zoom_par=1.
           max_step=(/50.,0.2,0.05/)
           TGM_best=TGM(1:3)
           TGM_best(3)=ABD(pos)
           flag_one=.FALSE.
           DO j=1,200

            IF(zoom_par<0.05) EXIT
            flag=.FALSE.
            DO i=1,9
             step=grid(i,:)*max_step*zoom_par
             TGM_test=TGM_best(1:3)+step
            IF(TGM_test(1)<MAXVAL(temp_gridL).AND.TGM_test(1)>MINVAL(temp_gridL)&
            &.AND.TGM_test(2)<MAXVAL(logg_gridL).AND.TGM_test(2)>MINVAL(logg_gridL)) THEN
             chisq=chisq_abd_e(TGM_test)

            IF(chisq<chisq_lim.AND.TGM_test(3)>TGM_best(3)) THEN
             TGM_best=TGM_test
             flag=.TRUE.
             flag_one=.TRUE.
              EXIT
            END IF
           END IF
           END DO
           IF(.not.flag) zoom_par=zoom_par*0.5
            IF(TGM_best(3)>0.8) THEN
             TGM_best(3)=10.
             EXIT
            END IF
          END DO

          IF(.NOT.flag_one) THEN
          up_ABD(pos)=9.
          ELSE
          up_ABD(pos)=TGM_best(3)
          END IF


        CALL update_grid(TGM)

        !estimate the lower [El/m] uncertains
           zoom_par=1.
           TGM_best=TGM(1:3)
           TGM_best(3)=ABD(pos)
           flag_one=.FALSE.
           DO j=1,200
            IF(zoom_par<0.01) EXIT
            flag=.FALSE.
            DO i=1,9
             step=grid(i,:)*max_step*zoom_par
             TGM_test=TGM_best(1:3)-step
            IF(TGM_test(1)<MAXVAL(temp_gridL).AND.TGM_test(1)>MINVAL(temp_gridL)&
            &.AND.TGM_test(2)<MAXVAL(logg_gridL).AND.TGM_test(2)>MINVAL(logg_gridL)) THEN
             chisq=chisq_abd_e(TGM_test)
            IF(chisq<chisq_lim.AND.TGM_test(3)<TGM_best(3)) THEN
             TGM_best=TGM_test
             flag=.TRUE.
             flag_one=.TRUE.
              EXIT
            END IF
           END IF
           END DO
           IF(.not.flag) zoom_par=zoom_par*0.5
            IF(TGM_best(3)<-0.6) THEN
             TGM_best(3)=-10.
             EXIT
            END IF
          END DO

          IF(.NOT.flag_one) THEN
          lo_ABD(pos)=-9.
          ELSE
          lo_ABD(pos)=TGM_best(3)
          END IF

       END SUBROUTINE ABD_e
!############################
SUBROUTINE TGM_errors(chisq_best)
        USE share, ONLY: ABD,TGM
        IMPLICIT NONE
        REAL(DP), INTENT(IN) :: chisq_best
        INTEGER (I1B) :: i


        !estimate the Teff errors
        pos=1
        CALL update_grid(TGM)
        CALL TGM_e(chisq_best)

        !estimate the logg errors
        pos=2
        CALL update_grid(TGM)
        CALL TGM_e(chisq_best)
        !estimate the [M/H] errors
        pos=3
        CALL update_grid(TGM)
        CALL TGM_e(chisq_best)

        !estimate the [El/m] errors
        DO i=1,INT(SIZE(ABD,1),I1B)
         IF(write_ABD_mask(i).AND.ABD_mask(i)) THEN
          pos=i
        CALL update_grid(TGM)
          CALL ABD_e(chisq_best)
         END IF
        END DO
END SUBROUTINE TGM_errors
!############################
SUBROUTINE make_grid(param,grid)
        IMPLICIT NONE
        INTEGER (I1B),INTENT(IN) :: param
        REAL (DP), DIMENSION(9,3), INTENT(OUT) :: grid
        REAL (DP), DIMENSION(9,3) :: grid_swap

      grid(1,:)=(/1.,0.,0./)
      grid(2,:)=(/0.5,0.5,-0.5/)
      grid(3,:)=(/0.5,0.5,0.5/)
      grid(4,:)=(/0.5,-0.5,0.5/)
      grid(5,:)=(/0.5,-0.5,-0.5/)
      grid(6,:)=(/0.01,1.,-1./)
      grid(7,:)=(/0.01,1.,1./)
      grid(8,:)=(/0.01,-1.,1./)
      grid(9,:)=(/0.01,-1.,-1./)

      IF(param==2) THEN
       grid_swap(:,1)=grid(:,3)
       grid_swap(:,2)=grid(:,1)
       grid_swap(:,3)=grid(:,2)
       grid=grid_swap
      ELSE IF(param==3) THEN
       grid_swap(:,1)=grid(:,2)
       grid_swap(:,2)=grid(:,3)
       grid_swap(:,3)=grid(:,1)
       grid=grid_swap
      END IF


END SUBROUTINE make_grid
!###################################################
        SUBROUTINE TGM_errors_null
        USE utils, ONLY: write_ABD_mask
        USE share, ONLY: ABD
        IMPLICIT NONE
        INTEGER(I1B) :: i

        DO i=1,3
         up_TGM(i)=9999._dp
         lo_TGM(i)=-9999._dp
        END DO
        DO i=1,INT(SIZE(ABD,1),I1B)
         IF(write_ABD_mask(i)) THEN
         up_ABD(i)=9.99_dp
         lo_ABD(i)=-9.99_dp
         END IF
        END DO

        END SUBROUTINE TGM_errors_null
!############################


END MODULE uncertains2