MODULE uncertains2
        USE nrtype
        USE make_model, ONLY: make_model_TGM,make_model_ABDerr
        USE share, ONLY: TGM,ABD,weights,sn,f_sp_norm,up_ABD,lo_ABD,sn_var
        USE error
        IMPLICIT NONE
        REAL(DP) :: chisq_lim,chisq_min
        INTEGER(I1B) :: posTGM,posABD
        REAL(DP), DIMENSION(3) :: up_TGM,lo_TGM
!        INTEGER(I2B) :: index


        CONTAINS

        SUBROUTINE TGM_errors(chisq_min)
        USE utils, ONLY: write_ABD_mask
        IMPLICIT NONE
        REAL(DP), INTENT(IN) :: chisq_min
        REAL(DP), DIMENSION(SIZE(f_sp_norm)) :: model
        REAL(DP), DIMENSION(SIZE(TGM)) :: parTGM
        REAL(DP), DIMENSION(SIZE(ABD)) :: parABD
        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(3),PARAMETER :: limit=(/10000.,6.,2.8/)
        REAL(DP), DIMENSION(4,3) :: incr
        REAL(DP), DIMENSION(4,3) :: p
        REAL(DP), DIMENSION(4) :: y
        REAL(DP) :: chi2,tol,minusone
        INTEGER(I1B) :: i,j,k
        INTEGER(I4B) :: iter

!        index=0
        minusone=-1

        !prob(nu) is chisq as function of the confidence level=68% 
        !where nu=degree of freedom
        chisq_lim=chisq_min+prob(1)

        tol=0.5_dp/chisq_lim

        incr(1,:)=(/0.,0.,0./)
        incr(2,:)=(/500.,0.,0./)
        incr(3,:)=(/0.,0.5,0./)
        incr(4,:)=(/0.,0.0,0.5/)


        DO k=1,3
        !compute the upper boundaries
          DO i=1,4
           DO j=1,3
            p(i,j)=TGM(j)+incr(i,j)
           END DO
          END DO
         parTGM=TGM
         posTGM=k
         DO i=1,4
          parTGM(1:3)=p(i,:)
          CALL make_model_TGM(model,parTGM)
          chi2=SUM(((f_sp_norm-model)*sn_var*weights)**2)
          IF(chi2>chisq_lim) THEN
           y(i)=1e9
          ELSE
           y(i)=-parTGM(posTGM)
          END IF
         END DO
         CALL amoeba(p,y,tol,func_sup_TGM,iter,k,limit(k))
         up_TGM(k)=MAXVAL(p(:,k))
        !compute the lower boundaries
          DO i=1,4
           DO j=1,3
            p(i,j)=TGM(j)-incr(i,j)
           END DO
          END DO
         parTGM=TGM
         posTGM=k
         DO i=1,4
          parTGM(1:3)=p(i,:)
          CALL make_model_TGM(model,parTGM)
          chi2=SUM(((f_sp_norm-model)*sn_var*weights)**2)
          IF(chi2>chisq_lim) THEN
           y(i)=1e9
          ELSE
           y(i)=parTGM(posTGM)
          END IF
         END DO
         CALL amoeba(p,y,tol,func_inf_TGM,iter,k,limit(k))
         lo_TGM(k)=MINVAL(p(:,k))
!        write(*,*) 'res TGM',k,lo_TGM(k)-TGM(k),up_TGM(k)-TGM(k)
!        write(*,*) 'res TGM',k,p(:,k),TGM(k)
        END DO

!####################################
!#### do the same for the abundances
!####################################

        DO k=1,INT(SIZE(ABD,1),I1B)
         IF(write_ABD_mask(k)) THEN
        !compute the upper boundaries
         posABD=k
          DO i=1,4
           DO j=1,2
            p(i,j)=TGM(j)+incr(i,j)
           END DO
           p(i,3)=ABD(k)+incr(i,3)
          END DO
         parTGM=TGM
         parABD=ABD
         posABD=k
         DO i=1,4
          parTGM(1:2)=p(i,1:2)
          parABD(posABD)=p(i,3)
          CALL make_model_ABDerr(model,parTGM,parABD)
          chi2=SUM(((f_sp_norm-model)*sn_var*weights)**2)
          IF(chi2>chisq_lim) THEN
           y(i)=1e9
          ELSE
           y(i)=-parABD(posABD)
          END IF
         END DO
         CALL amoeba(p,y,tol,func_sup_ABD,iter,3_I1B,1.0_dp)
         up_ABD(k)=MAXVAL(p(:,3))
        !compute the lower boundaries
          DO i=1,4
           DO j=1,2
            p(i,j)=TGM(j)-incr(i,j)
           END DO
           p(i,3)=ABD(k)-incr(i,3)
          END DO
         parTGM=TGM
         parABD=ABD
         posABD=k
         DO i=1,4
          parTGM(1:2)=p(i,1:2)
          parABD(posABD)=p(i,3)
          CALL make_model_ABDerr(model,parTGM,parABD)
          chi2=SUM(((f_sp_norm-model)*sn_var*weights)**2)
          IF(chi2>chisq_lim) THEN
           y(i)=1e9
          ELSE
           y(i)=parABD(posABD)
          END IF
         END DO
         CALL amoeba(p,y,tol,func_inf_ABD,iter,3_I1B,1.0_dp)
         lo_ABD(k)=MINVAL(p(:,3))
!        write(*,*) 'res_ABD',k,lo_ABD(k)-ABD(k),up_ABD(k)-ABD(k)
         END IF
        END DO
!        write(*,*) lo_ABD(1:5)-ABD(1:5),up_ABD(1:5)-ABD(1:5)

        END SUBROUTINE TGM_errors
!###################################################
        SUBROUTINE TGM_errors_quick(covarTGM,covarABD)
        USE utils, ONLY: write_ABD_mask
        IMPLICIT NONE
        REAL(DP), DIMENSION(:,:) :: covarTGM,covarABD
        INTEGER(I1B) :: i

        DO i=1,3
         up_TGM(i)=sqrt(covarTGM(i,i))+TGM(i)
         lo_TGM(i)=-sqrt(covarTGM(i,i))+TGM(i)
        END DO
        DO i=1,INT(SIZE(ABD,1),I1B)
         IF(write_ABD_mask(i)) THEN
         up_ABD(i)=sqrt(covarABD(i,i))
         lo_ABD(i)=-up_ABD(i)
         END IF
        END DO

        END SUBROUTINE TGM_errors_quick
!###################################################
        SUBROUTINE TGM_errors_null
        USE utils, ONLY: write_ABD_mask
        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

!########################
!########################
        FUNCTION func_sup_TGM(x)
        USE nrtype
        IMPLICIT NONE
        REAL(DP), DIMENSION(:), INTENT(IN) :: x
        REAL(DP) :: func_sup_TGM,chi2
        REAL(DP), DIMENSION(SIZE(TGM)) :: parTGM
        REAL(DP), DIMENSION(SIZE(f_sp_norm)) :: model

!        index=index+1
                
        parTGM=TGM
        parTGM(1:3)=x
        CALL make_model_TGM(model,parTGM)
        chi2=SUM(((f_sp_norm-model)*sn_var*weights)**2)
        IF(chi2>chisq_lim) THEN
        func_sup_TGM=1e9
        ELSE
        func_sup_TGM=-x(posTGM)
        END IF
!        write(*,*) 'supTGM',posTGM,x,func_sup_TGM

!        OPEN(10,file='amoeba.dat',action='WRITE',position='APPEND')
!        write(10,fmt='(I4,1X,F10.1,1X,F8.3,1X,F8.3,1X,F8.3,1X)') index,chi2,x
!        close(10)

        END FUNCTION func_sup_TGM
!########################
        FUNCTION func_inf_TGM(x)
        USE nrtype
        IMPLICIT NONE
        REAL(DP), DIMENSION(:), INTENT(IN) :: x
        REAL(DP) :: func_inf_TGM,chi2
        REAL(DP), DIMENSION(SIZE(TGM)) :: parTGM
        REAL(DP), DIMENSION(SIZE(f_sp_norm)) :: model

        parTGM=TGM
        parTGM(1:3)=x
        CALL make_model_TGM(model,parTGM)
        chi2=SUM(((f_sp_norm-model)*sn_var*weights)**2)
        IF(chi2>chisq_lim) THEN
        func_inf_TGM=1e9
        ELSE
        func_inf_TGM=x(posTGM)
        END IF
!        write(*,*) 'infTGM',posTGM,x,func_inf_TGM

        END FUNCTION func_inf_TGM
!########################
        FUNCTION func_sup_ABD(x)
        USE nrtype
        IMPLICIT NONE
        REAL(DP), DIMENSION(:), INTENT(IN) :: x
        REAL(DP) :: func_sup_ABD,chi2
        REAL(DP), DIMENSION(SIZE(TGM)) :: parTGM
        REAL(DP), DIMENSION(SIZE(ABD)) :: parABD
        REAL(DP), DIMENSION(SIZE(f_sp_norm)) :: model

        parTGM=TGM
        parTGM(1:2)=x(1:2)
        parABD=ABD
        parABD(posABD)=x(3)
        CALL make_model_ABDerr(model,parTGM,parABD)
        chi2=SUM(((f_sp_norm-model)*sn_var*weights)**2)
        IF(chi2>chisq_lim) THEN
        func_sup_ABD=1e9
        ELSE
        func_sup_ABD=-x(3)
        END IF
!        write(*,*) 'supABD',posABD,x,func_sup_ABD

        END FUNCTION func_sup_ABD
!########################
        FUNCTION func_inf_ABD(x)
        USE nrtype
        IMPLICIT NONE
        REAL(DP), DIMENSION(:), INTENT(IN) :: x
        REAL(DP) :: func_inf_ABD,chi2
        REAL(DP), DIMENSION(SIZE(TGM)) :: parTGM
        REAL(DP), DIMENSION(SIZE(ABD)) :: parABD
        REAL(DP), DIMENSION(SIZE(f_sp_norm)) :: model

        parTGM=TGM
        parABD=ABD
        parTGM(1:2)=x(1:2)
        parABD(posABD)=x(3)
        CALL make_model_ABDerr(model,parTGM,parABD)
        chi2=SUM(((f_sp_norm-model)*sn_var*weights)**2)
        IF(chi2>chisq_lim) THEN
        func_inf_ABD=1e9
        ELSE
        func_inf_ABD=x(3)
        END IF
!        write(*,*) 'infABD',posABD,x,func_inf_ABD

        END FUNCTION func_inf_ABD
!########################        
        SUBROUTINE amoeba(p,y,ftol,func,iter,param,limit)
        USE util_minim, ONLY : assert_eq,imaxloc,iminloc,swap
        IMPLICIT NONE
        INTEGER(I4B), INTENT(OUT) :: iter
        REAL(DP), INTENT(IN) :: ftol,limit
        INTEGER(I1B), INTENT(IN) :: param
        REAL(DP), DIMENSION(:), INTENT(INOUT) :: y
        REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: p
        INTERFACE
                FUNCTION func(x)
                USE nrtype
                IMPLICIT NONE
                REAL(DP), DIMENSION(:), INTENT(IN) :: x
                REAL(DP) :: func
                END FUNCTION func
        END INTERFACE
        INTEGER(I4B), PARAMETER :: ITMAX=500
        REAL(DP), PARAMETER :: TINY=1.0e-10_dp
        INTEGER(I4B) :: ihi,ndim
        REAL(DP), DIMENSION(size(p,2)) :: psum
        call amoeba_private
        CONTAINS
!BL
        SUBROUTINE amoeba_private
        IMPLICIT NONE
        INTEGER(I4B) :: i,ilo,inhi
        REAL(DP) :: rtol,ysave,ytry,ytmp
        ndim=assert_eq(size(p,2),size(p,1)-1,size(y)-1,'amoeba')
        iter=0
        psum(:)=sum(p(:,:),dim=1)
        do
                ilo=iminloc(y(:))
                ihi=imaxloc(y(:))
                ytmp=y(ihi)
                y(ihi)=y(ilo)
                inhi=imaxloc(y(:))
                y(ihi)=ytmp
                rtol=2.0_dp*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo))+TINY)
                !if all the values in p(:,param) are beyond the grid, means there are no boundaries
                !then exit
                if (all(abs(p(:,param))>limit)) then
                        call swap(y(1),y(ilo))
                        call swap(p(1,:),p(ilo,:))
                        RETURN
                end if
                if (rtol < ftol) then
                        call swap(y(1),y(ilo))
                        call swap(p(1,:),p(ilo,:))
                        RETURN
                end if
!                if (iter >= ITMAX) call nrerror('ITMAX exceeded in amoeba')
                if (iter >= ITMAX) then
                        call swap(y(1),y(ilo))
                        call swap(p(1,:),p(ilo,:))
!                write(*,*) posTGM,posABD,'ITMAX exceeded in amoeba'
                return
                end if
                ytry=amotry(-1.0_dp)
                iter=iter+1
                if (ytry <= y(ilo)) then
                        ytry=amotry(2.0_dp)
                        iter=iter+1
                else if (ytry >= y(inhi)) then
                        ysave=y(ihi)
                        ytry=amotry(0.5_dp)
                        iter=iter+1
                        if (ytry >= ysave) then
                                p(:,:)=0.5_dp*(p(:,:)+spread(p(ilo,:),1,size(p,1)))
                                do i=1,ndim+1
                                        if (i /= ilo) y(i)=func(p(i,:))
                                end do
                                iter=iter+ndim
                                psum(:)=sum(p(:,:),dim=1)
                        end if
                end if
        end do
        END SUBROUTINE amoeba_private
!BL
        FUNCTION amotry(fac)
        IMPLICIT NONE
        REAL(DP), INTENT(IN) :: fac
        REAL(DP) :: amotry
        REAL(DP) :: fac1,fac2,ytry
        REAL(DP), DIMENSION(size(p,2)) :: ptry
        fac1=(1.0_dp-fac)/ndim
        fac2=fac1-fac
        ptry(:)=psum(:)*fac1-p(ihi,:)*fac2
        ytry=func(ptry)
        if (ytry < y(ihi)) then
                y(ihi)=ytry
                psum(:)=psum(:)-p(ihi,:)+ptry(:)
                p(ihi,:)=ptry(:)
        end if
        amotry=ytry
        END FUNCTION amotry
        END SUBROUTINE amoeba

END MODULE uncertains2