Homec4science
Diffusion FENNECS lite (3D) Browse

fennecslite/src/distrib_mod.f90
6b37315661933D

Last Change
Raw File

distrib_mod.f90
View Options

!------------------------------------------------------------------------------
! 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_mod
!
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
c4science · Help