Homec4science
Diffusion FENNECS: Finite Element Non-Neutral Electron Cloud Simulator (stable) Browse

fennecs/dependencies/incomplete_gamma/incomplete_gamma_mod.f90
7d9796852d52stable

Last Change
Raw File

incomplete_gamma_mod.f90
View Options

!------------------------------------------------------------------------------
! EPFL/Swiss Plasma Center
!------------------------------------------------------------------------------
!
! MODULE: incomplete_gamma
!
!> @author
!> A. Macleod
!
!> Last modif.
!> 2008
!
! DESCRIPTION:
!
!> Module containing functions enabling to compute
!> incomplete gamma values to generate random values following gamma
!> distribution (see ion_induced_mod.f90 -> gen_E_gamma)
!------------------------------------------------------------------------------
MODULE incomplete_gamma
CONTAINS
function alngam ( xvalue, ifault )
!*****************************************************************************80
!
!! ALNGAM computes the logarithm of the gamma function.
!
! Modified:
!
! 13 January 2008
!
! Author:
!
! Allan Macleod
! FORTRAN90 version by John Burkardt
!
! Reference:
!
! Allan Macleod,
! Algorithm AS 245,
! A Robust and Reliable Algorithm for the Logarithm of the Gamma Function,
! Applied Statistics,
! Volume 38, Number 2, 1989, pages 397-402.
!
! Parameters:
!
! Input, real ( kind = 8 ) XVALUE, the argument of the Gamma function.
!
! Output, integer ( kind = 4 ) IFAULT, error flag.
! 0, no error occurred.
! 1, XVALUE is less than or equal to 0.
! 2, XVALUE is too big.
!
! Output, real ( kind = 8 ) ALNGAM, the logarithm of the gamma function of X.
!
implicit none
real ( kind = 8 ) alngam
real ( kind = 8 ), parameter :: alr2pi = 0.918938533204673D+00
integer ( kind = 4 ) ifault
real ( kind = 8 ), dimension ( 9 ) :: r1 = (/ &
-2.66685511495D+00, &
-24.4387534237D+00, &
-21.9698958928D+00, &
11.1667541262D+00, &
3.13060547623D+00, &
0.607771387771D+00, &
11.9400905721D+00, &
31.4690115749D+00, &
15.2346874070D+00 /)
real ( kind = 8 ), dimension ( 9 ) :: r2 = (/ &
-78.3359299449D+00, &
-142.046296688D+00, &
137.519416416D+00, &
78.6994924154D+00, &
4.16438922228D+00, &
47.0668766060D+00, &
313.399215894D+00, &
263.505074721D+00, &
43.3400022514D+00 /)
real ( kind = 8 ), dimension ( 9 ) :: r3 = (/ &
-2.12159572323D+05, &
2.30661510616D+05, &
2.74647644705D+04, &
-4.02621119975D+04, &
-2.29660729780D+03, &
-1.16328495004D+05, &
-1.46025937511D+05, &
-2.42357409629D+04, &
-5.70691009324D+02 /)
real ( kind = 8 ), dimension ( 5 ) :: r4 = (/ &
0.279195317918525D+00, &
0.4917317610505968D+00, &
0.0692910599291889D+00, &
3.350343815022304D+00, &
6.012459259764103D+00 /)
real ( kind = 8 ) x
real ( kind = 8 ) x1
real ( kind = 8 ) x2
real ( kind = 8 ), parameter :: xlge = 5.10D+05
real ( kind = 8 ), parameter :: xlgst = 1.0D+30
real ( kind = 8 ) xvalue
real ( kind = 8 ) y
x = xvalue
alngam = 0.0D+00
!
! Check the input.
!
if ( xlgst <= x ) then
ifault = 2
return
end if
if ( x <= 0.0D+00 ) then
ifault = 1
return
end if
ifault = 0
!
! Calculation for 0 < X < 0.5 and 0.5 <= X < 1.5 combined.
!
if ( x < 1.5D+00 ) then
if ( x < 0.5D+00 ) then
alngam = - log ( x )
y = x + 1.0D+00
!
! Test whether X < machine epsilon.
!
if ( y == 1.0D+00 ) then
return
end if
else
alngam = 0.0D+00
y = x
x = ( x - 0.5D+00 ) - 0.5D+00
end if
alngam = alngam + x * (((( &
r1(5) * y &
+ r1(4) ) * y &
+ r1(3) ) * y &
+ r1(2) ) * y &
+ r1(1) ) / (((( &
y &
+ r1(9) ) * y &
+ r1(8) ) * y &
+ r1(7) ) * y &
+ r1(6) )
return
end if
!
! Calculation for 1.5 <= X < 4.0.
!
if ( x < 4.0D+00 ) then
y = ( x - 1.0D+00 ) - 1.0D+00
alngam = y * (((( &
r2(5) * x &
+ r2(4) ) * x &
+ r2(3) ) * x &
+ r2(2) ) * x &
+ r2(1) ) / (((( &
x &
+ r2(9) ) * x &
+ r2(8) ) * x &
+ r2(7) ) * x &
+ r2(6) )
!
! Calculation for 4.0 <= X < 12.0.
!
else if ( x < 12.0D+00 ) then
alngam = (((( &
r3(5) * x &
+ r3(4) ) * x &
+ r3(3) ) * x &
+ r3(2) ) * x &
+ r3(1) ) / (((( &
x &
+ r3(9) ) * x &
+ r3(8) ) * x &
+ r3(7) ) * x &
+ r3(6) )
!
! Calculation for 12.0 <= X.
!
else
y = log ( x )
alngam = x * ( y - 1.0D+00 ) - 0.5D+00 * y + alr2pi
if ( x <= xlge ) then
x1 = 1.0D+00 / x
x2 = x1 * x1
alngam = alngam + x1 * ( ( &
r4(3) * &
x2 + r4(2) ) * &
x2 + r4(1) ) / ( ( &
x2 + r4(5) ) * &
x2 + r4(4) )
end if
end if
return
end function alngam
function gamain ( x, p, ifault )
!*****************************************************************************80
!
!! GAMAIN computes the incomplete gamma ratio.
!
! Discussion:
!
! A series expansion is used if P > X or X <= 1. Otherwise, a
! continued fraction approximation is used.
!
! Modified:
!
! 17 January 2008
!
! Author:
!
! G Bhattacharjee
! FORTRAN90 version by John Burkardt
!
! Reference:
!
! G Bhattacharjee,
! Algorithm AS 32:
! The Incomplete Gamma Integral,
! Applied Statistics,
! Volume 19, Number 3, 1970, pages 285-287.
!
! Parameters:
!
! Input, real ( kind = 8 ) X, P, the parameters of the incomplete
! gamma ratio. 0 <= X, and 0 < P.
!
! Output, integer ( kind = 4 ) IFAULT, error flag.
! 0, no errors.
! 1, P <= 0.
! 2, X < 0.
! 3, underflow.
! 4, error return from the Log Gamma routine.
!
! Output, real ( kind = 8 ) GAMAIN, the value of the incomplete
! gamma ratio.
!
implicit none
real ( kind = 8 ) a
real ( kind = 8 ), parameter :: acu = 1.0D-08
!real ( kind = 8 ) alngam
real ( kind = 8 ) an
real ( kind = 8 ) arg
real ( kind = 8 ) b
real ( kind = 8 ) dif
real ( kind = 8 ) factor
real ( kind = 8 ) g
real ( kind = 8 ) gamain
real ( kind = 8 ) gin
integer ( kind = 4 ) i
integer ( kind = 4 ) ifault
real ( kind = 8 ), parameter :: oflo = 1.0D+37
real ( kind = 8 ) p
real ( kind = 8 ) pn(6)
real ( kind = 8 ) rn
real ( kind = 8 ) term
real ( kind = 8 ), parameter :: uflo = 1.0D-37
real ( kind = 8 ) x
!
! Check the input.
!
if ( p <= 0.0D+00 ) then
ifault = 1
gamain = 0.0D+00
return
end if
if ( x < 0.0D+00 ) then
ifault = 2
gamain = 0.0D+00
return
end if
if ( x == 0.0D+00 ) then
ifault = 0
gamain = 0.0D+00
return
end if
g = alngam ( p, ifault )
if ( ifault /= 0 ) then
ifault = 4
gamain = 0.0D+00
return
end if
arg = p * log ( x ) - x - g
if ( arg < log ( uflo ) ) then
ifault = 3
gamain = 0.0D+00
return
end if
ifault = 0
factor = exp ( arg )
!
! Calculation by series expansion.
!
if ( x <= 1.0D+00 .or. x < p ) then
gin = 1.0D+00
term = 1.0D+00
rn = p
do
rn = rn + 1.0D+00
term = term * x / rn
gin = gin + term
if ( term <= acu ) then
exit
end if
end do
gamain = gin * factor / p
return
end if
!
! Calculation by continued fraction.
!
a = 1.0D+00 - p
b = a + x + 1.0D+00
term = 0.0D+00
pn(1) = 1.0D+00
pn(2) = x
pn(3) = x + 1.0D+00
pn(4) = x * b
gin = pn(3) / pn(4)
do
a = a + 1.0D+00
b = b + 2.0D+00
term = term + 1.0D+00
an = a * term
do i = 1, 2
pn(i+4) = b * pn(i+2) - an * pn(i)
end do
if ( pn(6) /= 0.0D+00 ) then
rn = pn(5) / pn(6)
dif = abs ( gin - rn )
!
! Absolute error tolerance satisfied?
!
if ( dif <= acu ) then
!
! Relative error tolerance satisfied?
!
if ( dif <= acu * rn ) then
gamain = 1.0D+00 - factor * gin
exit
end if
end if
gin = rn
end if
do i = 1, 4
pn(i) = pn(i+2)
end do
if ( oflo <= abs ( pn(5) ) ) then
do i = 1, 4
pn(i) = pn(i) / oflo
end do
end if
end do
return
end function gamain
subroutine gamma_inc_values ( n_data, a, x, fx )
!*****************************************************************************80
!
!! GAMMA_INC_VALUES returns some values of the incomplete Gamma function.
!
! Discussion:
!
! The (normalized) incomplete Gamma function P(A,X) is defined as:
!
! PN(A,X) = 1/Gamma(A) * Integral ( 0 <= T <= X ) T**(A-1) * exp(-T) dT.
!
! With this definition, for all A and X,
!
! 0 <= PN(A,X) <= 1
!
! and
!
! PN(A,INFINITY) = 1.0
!
! In Mathematica, the function can be evaluated by:
!
! 1 - GammaRegularized[A,X]
!
! Modified:
!
! 20 November 2004
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Milton Abramowitz, Irene Stegun,
! Handbook of Mathematical Functions,
! National Bureau of Standards, 1964,
! ISBN: 0-486-61272-4,
! LC: QA47.A34.
!
! Stephen Wolfram,
! The Mathematica Book,
! Fourth Edition,
! Cambridge University Press, 1999,
! ISBN: 0-521-64314-7,
! LC: QA76.95.W65.
!
! Parameters:
!
! Input/output, integer ( kind = 4 ) N_DATA. The user sets N_DATA to 0
! before the first call. On each call, the routine increments N_DATA by 1,
! and returns the corresponding data; when there is no more data, the
! output value of N_DATA will be 0 again.
!
! Output, real ( kind = 8 ) A, the parameter of the function.
!
! Output, real ( kind = 8 ) X, the argument of the function.
!
! Output, real ( kind = 8 ) FX, the value of the function.
!
implicit none
integer ( kind = 4 ), parameter :: n_max = 20
real ( kind = 8 ) a
real ( kind = 8 ), save, dimension ( n_max ) :: a_vec = (/ &
0.10D+00, &
0.10D+00, &
0.10D+00, &
0.50D+00, &
0.50D+00, &
0.50D+00, &
0.10D+01, &
0.10D+01, &
0.10D+01, &
0.11D+01, &
0.11D+01, &
0.11D+01, &
0.20D+01, &
0.20D+01, &
0.20D+01, &
0.60D+01, &
0.60D+01, &
0.11D+02, &
0.26D+02, &
0.41D+02 /)
real ( kind = 8 ) fx
real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ &
0.7382350532339351D+00, &
0.9083579897300343D+00, &
0.9886559833621947D+00, &
0.3014646416966613D+00, &
0.7793286380801532D+00, &
0.9918490284064973D+00, &
0.9516258196404043D-01, &
0.6321205588285577D+00, &
0.9932620530009145D+00, &
0.7205974576054322D-01, &
0.5891809618706485D+00, &
0.9915368159845525D+00, &
0.1018582711118352D-01, &
0.4421745996289254D+00, &
0.9927049442755639D+00, &
0.4202103819530612D-01, &
0.9796589705830716D+00, &
0.9226039842296429D+00, &
0.4470785799755852D+00, &
0.7444549220718699D+00 /)
integer ( kind = 4 ) n_data
real ( kind = 8 ) x
real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ &
0.30D-01, &
0.30D+00, &
0.15D+01, &
0.75D-01, &
0.75D+00, &
0.35D+01, &
0.10D+00, &
0.10D+01, &
0.50D+01, &
0.10D+00, &
0.10D+01, &
0.50D+01, &
0.15D+00, &
0.15D+01, &
0.70D+01, &
0.25D+01, &
0.12D+02, &
0.16D+02, &
0.25D+02, &
0.45D+02 /)
if ( n_data < 0 ) then
n_data = 0
end if
n_data = n_data + 1
if ( n_max < n_data ) then
n_data = 0
a = 0.0D+00
x = 0.0D+00
fx = 0.0D+00
else
a = a_vec(n_data)
x = x_vec(n_data)
fx = fx_vec(n_data)
end if
return
end subroutine gamma_inc_values
subroutine timestamp ( )
!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
! Example:
!
! 31 May 2001 9:45:54.872 AM
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 18 May 2013
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! None
!
implicit none
character ( len = 8 ) ampm
integer ( kind = 4 ) d
integer ( kind = 4 ) h
integer ( kind = 4 ) m
integer ( kind = 4 ) mm
character ( len = 9 ), parameter, dimension(12) :: month = (/ &
'January ', 'February ', 'March ', 'April ', &
'May ', 'June ', 'July ', 'August ', &
'September', 'October ', 'November ', 'December ' /)
integer ( kind = 4 ) n
integer ( kind = 4 ) s
integer ( kind = 4 ) values(8)
integer ( kind = 4 ) y
call date_and_time ( values = values )
y = values(1)
m = values(2)
d = values(3)
h = values(5)
n = values(6)
s = values(7)
mm = values(8)
if ( h < 12 ) then
ampm = 'AM'
else if ( h == 12 ) then
if ( n == 0 .and. s == 0 ) then
ampm = 'Noon'
else
ampm = 'PM'
end if
else
h = h - 12
if ( h < 12 ) then
ampm = 'PM'
else if ( h == 12 ) then
if ( n == 0 .and. s == 0 ) then
ampm = 'Midnight'
else
ampm = 'AM'
end if
end if
end if
write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )
return
end subroutine timestamp
END MODULE incomplete_gamma
c4science · Help