!------------------------------------------------------------------------------
! EPFL/Swiss Plasma Center
!------------------------------------------------------------------------------
!
! MODULE: distrib
!
!> @author
!> Guillaume Le Bars EPFL/SPC
!> Trach Minh Tran   EPFL/SPC
!
! DESCRIPTION:
!> Contains the routines used to load several type of distribution functions
!------------------------------------------------------------------------------

MODULE distrib
    USE constants
    USE omp_lib
    USE random
    !
    IMPLICIT NONE
contains

!---------------------------------------------------------------------------
!> @author
!> Trach Minh Tran EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Load a uniform distribution using the Hammersley's sequence (0<=y<=1).
!> (nbase = 0 => Random sampling)
!
!> @param [in] nbase Base of the Hammersley's sequence. (0 => Random)
!> @param [out] y array containing the random variables.
!---------------------------------------------------------------------------
SUBROUTINE loduni(nbase, y)
    !
    !   Load a uniform distribution using the Hammersley's sequence.
    !   (nbase = 0 => Random sampling)
    !
        INTEGER, INTENT(in)        :: nbase
        REAL(kind=db), INTENT(inout) :: y(:)
        INTEGER       :: n, i, ompthread
    !
        n = SIZE(y)
    !
        SELECT CASE (nbase)
        CASE(0)
            ompthread=omp_get_thread_num()+1
            CALL random_array(y,n,ran_index(ompthread),ran_array(:,ompthread))
           !CALL RANDOM_NUMBER(y)
        CASE(1)
           DO i=1,n
              y(i) = (i-0.5_db)/n
           END DO
        CASE(2:)
           DO i=1,n
              y(i) = rev(nbase,i)
           END DO
        CASE default
           WRITE(*,'(a,i5)') 'Invalid value of NBASE =', nbase
        END SELECT
    !
      END SUBROUTINE loduni
    
    !---------------------------------------------------------------------------
    !> @author
    !> Trach Minh Tran EPFL/SPC
    !
    ! DESCRIPTION:
    !>
    !> @brief Load a gaussian distribution using a uniform distribution according to the Hammersley's sequence.
    !>   (nbase = 0 => Random sampling)
    !
    !> @param [in] nbase Base of the Hammersley's sequence. (0 => Random)
    !> @param [out] y array containing the random variables.
    !---------------------------------------------------------------------------
      SUBROUTINE lodgaus(nbase, y)
    !
    !   Sample y from the Gauss distributrion
    !
        REAL(kind=db), INTENT(inout) :: y(:)
        INTEGER, INTENT(in)        :: nbase
    !
        INTEGER       :: n, i, iflag
        REAL(kind=db) :: r(SIZE(y))
        REAL(kind=db) :: t
    !
        REAL(kind=db) :: c0, c1, c2
        REAL(kind=db) :: d1, d2, d3
        DATA c0, c1, c2/2.515517_db, 0.802853_db, 0.010328_db/
        DATA d1, d2, d3/1.432788_db, 0.189269_db, 0.001308_db/
    !
        n = SIZE(y)
        CALL loduni(nbase, r)
        r = 1.0E-5 + 0.99998*r
    !
        DO i=1,n
           iflag = 1
           IF (r(i) .GT. 0.5) THEN
              r(i) = 1.0 - r(i)
              iflag = -1
           END IF
           t = SQRT(LOG(1.0/(r(i)*r(i))))
           y(i) = t - (c0+c1*t+c2*t**2) / (1.0+d1*t+d2*t**2+d3*t**3)
           y(i) = y(i) * iflag
        END DO
        y = y - SUM(y)/REAL(n)
        
      END SUBROUTINE lodgaus
      !---------------------------------------------------------------------------
      !> @author
      !> Guillaume Le Bars EPFL/SPC
      !
      ! DESCRIPTION:
      !>
      !> @brief Load the array y according to a probability distribution function proportional to 1/r
      !> between ra and rb
      !
      !> @param [in] nbase Base of the Hammersley's sequence. (0 => Random)
      !> @param [in] ra Lower limit for the values in y
      !> @param [in] rb Upper limit for the values in y
      !> @param [out] y array containing the random variables.
      !---------------------------------------------------------------------------
      SUBROUTINE lodinvr(nbase, y, ra, rb)
        !
        !
            REAL(kind=db), INTENT(out) :: y(:)
            INTEGER, INTENT(in)        :: nbase
            REAL(kind=db), INTENT(in)  :: rb, ra
        !
            REAL(kind=db) :: r(SIZE(y))
        !
            CALL loduni(nbase, r)
            r=r*log(rb/ra)
            y=ra*exp(r)
      END SUBROUTINE lodinvr
      !---------------------------------------------------------------------------
      !> @author
      !> Guillaume Le Bars EPFL/SPC
      !
      ! DESCRIPTION:
      !>
      !> @brief Load the array y according to a probability distribution function proportional to r
      !> between ra and rb
      !
      !> @param [in] nbase Base of the Hammersley's sequence. (0 => Random)
      !> @param [in] ra Lower limit for the values in y
      !> @param [in] rb Upper limit for the values in y
      !> @param [out] y array containing the random variables.
      !---------------------------------------------------------------------------
      SUBROUTINE lodlinr(nbase, y, ra, rb)
            REAL(kind=db), INTENT(out) :: y(:)
            INTEGER, INTENT(in)        :: nbase
            REAL(kind=db), INTENT(in)  :: rb, ra
        !
            REAL(kind=db) :: r(SIZE(y))
        !
            CALL loduni(nbase, r)
            y=sqrt(r*(rb**2-ra**2)+ra**2)
      END SUBROUTINE lodlinr
      !---------------------------------------------------------------------------
      !> @author
      !> Guillaume Le Bars EPFL/SPC
      !
      ! DESCRIPTION:
      !>
      !> @brief Load the array y according to a uniform probability distribution function
      !> between ra and rb
      !
      !> @param [in] nbase Base of the Hammersley's sequence. (0 => Random)
      !> @param [in] ra Lower limit for the values in y
      !> @param [in] rb Upper limit for the values in y
      !> @param [out] y array containing the random variables.
      !---------------------------------------------------------------------------
      SUBROUTINE lodunir(nbase, y, ra, rb)
            REAL(kind=db), INTENT(out) :: y(:)
            INTEGER, INTENT(in)        :: nbase
            REAL(kind=db), INTENT(in)  :: rb, ra
        !
            REAL(kind=db) :: r(SIZE(y))
        !
            CALL loduni(nbase, r)
            y=ra+(rb-ra)*r
      END SUBROUTINE lodunir
      !---------------------------------------------------------------------------
      !> @author
      !> Guillaume Le Bars EPFL/SPC
      !
      ! DESCRIPTION:
      !>
      !> @brief Load the array y according to a gaussian probability distribution function
      !> around the mean of rb and ra and a sigma of (rb-ra)/10
      !
      !> @param [in] nbase Base of the Hammersley's sequence. (0 => Random)
      !> @param [in] ra Lower limit for the values in y
      !> @param [in] rb Upper limit for the values in y
      !> @param [out] y array containing the random variables.
      !---------------------------------------------------------------------------
      SUBROUTINE lodgausr(nbase, y, ra, rb)
            REAL(kind=db), INTENT(out) :: y(:)
            INTEGER, INTENT(in)        :: nbase
            REAL(kind=db), INTENT(in)  :: rb, ra
        !
            REAL(kind=db) :: r(SIZE(y))
        !
            CALL lodgaus(nbase, r)
            y=0.5*(rb+ra)+ 0.1*(rb-ra)*r
      END SUBROUTINE lodgausr
    !---------------------------------------------------------------------------
    !> @author
    !> Trach Minh Tran EPFL/SPC
    !
    ! DESCRIPTION:
    !>
    !> @brief Return an element of the Hammersley's sequence
    !
    !> @param [in] nbase Base of the Hammersley's sequence. (0 => Random)
    !> @param [in] i Index of the sequence.
    !> @param [out] rev Returned element of the Hammersley's sequence 
    !---------------------------------------------------------------------------
       REAL(kind=db) FUNCTION rev(nbase,i)
    !
        INTEGER, INTENT(IN) :: nbase   ! Base of the Hammersley's sequence (IN)
        INTEGER, INTENT(IN) :: i       ! Index of the sequence (IN)
    !
    !   Local vars
        INTEGER :: j1, j2
        REAL(kind=db) :: xs, xsi
    !-----------------------------------------------------------------------
        xs = 0.0
        xsi = 1.0
        j2 = i
        DO
           xsi = xsi/nbase
           j1 = j2/nbase
           xs = xs + (j2-nbase*j1)*xsi
           j2 = j1
           IF( j2.LE.0 )  EXIT
        END DO
        rev = xs
      END FUNCTION rev
    !---------------------------------------------------------------------------
    !> @author
    !> Trach Minh Tran EPFL/SPC
    !
    ! DESCRIPTION:
    !>
    !> @brief Modified Bessel functions of the first kind of the first order 1
    !
    !---------------------------------------------------------------------------
      
      FUNCTION bessi1(x)
        REAL(kind=db) ::  bessi1,x
        REAL(kind=db) ::  ax
        REAL(kind=db) p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y
        SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
        DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/
        DATA q1,q2,q3,q4,q5,q6,q7,q8,q9 /0.39894228d0, -0.3988024d-1, -0.362018d-2,  0.163801d-2,-0.1031555d-1, &
             &                           0.2282967d-1, -0.2895312d-1,  0.1787654d-1,-0.420059d-2/
        if (abs(x).lt.3.75) then
           y=(x/3.75)**2
           bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
        else
           ax=abs(x)
           y=3.75/ax
           bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))))
           if(x.lt.0.)bessi1=-bessi1
        endif
        return
      END FUNCTION bessi1


      ! find the index for the value in array closest but inferior to val
    function closest(array, val, siz)
        real(kind=db):: array(1:), val
        integer:: closest
        integer:: siz
        integer:: i, j, mid
  
        ! Corner cases
        if (val.lt.array(1)) then
          closest=-1
          RETURN
        end if
        if (val .ge. array(siz)) then
          closest=siz
          return 
        end if
   
      !Start binary search
      i=1
      j=siz
      mid=1
      do while (i .lt. j) 
          mid = (i + j) / 2
   
          if (val .ge. array(mid) .and. val .lt. array(mid+1)) then
            closest = mid
            return
          end if
          
          ! if val > array(mid) search in right 
          if (val .gt. array(mid)) then
            ! update i
              i = mid+1
          else  ! if val < array(mid) search in left
            ! update j
              j = mid 
          end if
      end do
   
      closest=-1 ! not found
    end function closest

END MODULE distrib