Merge branch '12-fixOrientationSampling' into 'development'
Resolve "Fix and improve random orientation generation" Closes #12 See merge request damask/DAMASK!21
This commit is contained in:
Franz Roters
commit
0e4812014a
|
|
@ -324,6 +324,13 @@ HybridIA:
|
|||
- master
|
||||
- release
|
||||
|
||||
TextureComponents:
|
||||
stage: spectral
|
||||
script: TextureComponents/test.py
|
||||
except:
|
||||
- master
|
||||
- release
|
||||
|
||||
###################################################################################################
|
||||
Marc_compileIfort2014:
|
||||
stage: compileMarc2014
|
||||
|
|
|
|||
|
|
@ -23,7 +23,7 @@ if which $1 &> /dev/null; then
|
|||
$1 $2
|
||||
echo -e '\n'
|
||||
else
|
||||
echo $ does not exist
|
||||
echo $1 not found
|
||||
fi
|
||||
}
|
||||
|
||||
|
|
|
|||
556
src/math.f90
556
src/math.f90
|
|
@ -162,10 +162,8 @@ module math
|
|||
math_rotate_forward3333, &
|
||||
math_limit
|
||||
private :: &
|
||||
halton, &
|
||||
halton_memory, &
|
||||
halton_ndim_set, &
|
||||
halton_seed_set
|
||||
math_check, &
|
||||
halton
|
||||
|
||||
contains
|
||||
|
||||
|
|
@ -217,10 +215,6 @@ subroutine math_init
|
|||
write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest
|
||||
|
||||
call random_seed(put = randInit)
|
||||
|
||||
call halton_seed_set(int(randInit(1), pInt))
|
||||
call halton_ndim_set(3_pInt)
|
||||
|
||||
call math_check()
|
||||
|
||||
end subroutine math_init
|
||||
|
|
@ -284,7 +278,7 @@ subroutine math_check
|
|||
endif
|
||||
|
||||
! +++ check rotation sense of q and R +++
|
||||
call halton(3_pInt,v) ! random vector
|
||||
v = halton([2_pInt,8_pInt,5_pInt]) ! random vector
|
||||
R = math_qToR(q)
|
||||
if (any(abs(math_mul33x3(R,v) - math_qRot(q,v)) > tol_math_check)) then
|
||||
write (error_msg, '(a)' ) 'R(q)*v has different sense than q*v'
|
||||
|
|
@ -1238,7 +1232,7 @@ function math_qRand()
|
|||
real(pReal), dimension(4) :: math_qRand
|
||||
real(pReal), dimension(3) :: rnd
|
||||
|
||||
call halton(3_pInt,rnd)
|
||||
rnd = halton([8_pInt,4_pInt,9_pInt])
|
||||
math_qRand = [cos(2.0_pReal*PI*rnd(1))*sqrt(rnd(3)), &
|
||||
sin(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), &
|
||||
cos(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), &
|
||||
|
|
@ -1761,7 +1755,7 @@ function math_sampleRandomOri()
|
|||
implicit none
|
||||
real(pReal), dimension(3) :: math_sampleRandomOri, rnd
|
||||
|
||||
call halton(3_pInt,rnd)
|
||||
rnd = halton([1_pInt,7_pInt,3_pInt])
|
||||
math_sampleRandomOri = [rnd(1)*2.0_pReal*PI, &
|
||||
acos(2.0_pReal*rnd(2)-1.0_pReal), &
|
||||
rnd(3)*2.0_pReal*PI]
|
||||
|
|
@ -1770,118 +1764,96 @@ end function math_sampleRandomOri
|
|||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief draw a random sample from Gauss component with noise (in radians) half-width
|
||||
!> @brief draw a sample from an Gaussian distribution around given orientation and Full Width
|
||||
! at Half Maximum (FWHM)
|
||||
!> @details: A uniform misorientation (limited to 2*FWHM) is sampled followed by convolution with
|
||||
! a Gausian distribution
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function math_sampleGaussOri(center,noise)
|
||||
use prec, only: &
|
||||
tol_math_check
|
||||
function math_sampleGaussOri(center,FWHM)
|
||||
|
||||
implicit none
|
||||
real(pReal), intent(in) :: noise
|
||||
real(pReal), intent(in) :: FWHM
|
||||
real(pReal), dimension(3), intent(in) :: center
|
||||
real(pReal) :: cosScatter,scatter
|
||||
real(pReal), dimension(3) :: math_sampleGaussOri, disturb
|
||||
real(pReal), dimension(3), parameter :: ORIGIN = 0.0_pReal
|
||||
real(pReal), dimension(5) :: rnd
|
||||
real(pReal) :: angle
|
||||
real(pReal), dimension(3) :: math_sampleGaussOri, axis
|
||||
real(pReal), dimension(4) :: rnd
|
||||
real(pReal), dimension(3,3) :: R
|
||||
|
||||
noScatter: if (abs(noise) < tol_math_check) then
|
||||
if (FWHM < 0.1_pReal*INRAD) then
|
||||
math_sampleGaussOri = center
|
||||
else noScatter
|
||||
! Helming uses different distribution with Bessel functions
|
||||
! therefore the gauss scatter width has to be scaled differently
|
||||
scatter = 0.95_pReal * noise
|
||||
cosScatter = cos(scatter)
|
||||
else
|
||||
GaussConvolution: do
|
||||
rnd = halton([8_pInt,3_pInt,6_pInt,11_pInt])
|
||||
axis(1) = rnd(1)*2.0_pReal-1.0_pReal ! uniform on [-1,1]
|
||||
axis(2:3) = [sqrt(1.0-axis(1)**2.0_pReal)*cos(rnd(2)*2.0*PI),&
|
||||
sqrt(1.0-axis(1)**2.0_pReal)*sin(rnd(2)*2.0*PI)] ! random axis
|
||||
angle = (rnd(3)-0.5_pReal)*4.0_pReal*FWHM ! rotation by [0, +-2 FWHM]
|
||||
R = math_axisAngleToR(axis,angle)
|
||||
angle = math_EulerMisorientation([0.0_pReal,0.0_pReal,0.0_pReal],math_RtoEuler(R))
|
||||
if (rnd(4) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) exit ! rejection sampling (Gaussian)
|
||||
enddo GaussConvolution
|
||||
math_sampleGaussOri = math_RtoEuler(math_mul33x33(R,math_EulerToR(center)))
|
||||
endif
|
||||
|
||||
do
|
||||
call halton(5_pInt,rnd)
|
||||
rnd(1:3) = 2.0_pReal*rnd(1:3)-1.0_pReal ! expand 1:3 to range [-1,+1]
|
||||
disturb = [ scatter * rnd(1), & ! phi1
|
||||
sign(1.0_pReal,rnd(2))*acos(cosScatter+(1.0_pReal-cosScatter)*rnd(4)), & ! Phi
|
||||
scatter * rnd(3)] ! phi2
|
||||
if (rnd(5) <= exp(-1.0_pReal*(math_EulerMisorientation(ORIGIN,disturb)/scatter)**2_pReal)) exit
|
||||
enddo
|
||||
|
||||
math_sampleGaussOri = math_RtoEuler(math_mul33x33(math_EulerToR(disturb),math_EulerToR(center)))
|
||||
endif noScatter
|
||||
|
||||
end function math_sampleGaussOri
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief draw a random sample from Fiber component with noise (in radians)
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function math_sampleFiberOri(alpha,beta,noise)
|
||||
use prec, only: &
|
||||
tol_math_check
|
||||
!> @brief draw a sample from an Gaussian distribution around given fiber texture and Full Width
|
||||
! at Half Maximum (FWHM)
|
||||
!-------------------------------------------------------------------------------------------------
|
||||
function math_sampleFiberOri(alpha,beta,FWHM)
|
||||
|
||||
implicit none
|
||||
real(pReal), dimension(3) :: math_sampleFiberOri, fiberInC,fiberInS,axis
|
||||
real(pReal), dimension(2), intent(in) :: alpha,beta
|
||||
real(pReal), dimension(6) :: rnd
|
||||
real(pReal), dimension(3,3) :: oRot,fRot,pRot
|
||||
real(pReal) :: noise, scatter, cos2Scatter, angle
|
||||
integer(pInt), dimension(2,3), parameter :: ROTMAP = reshape([2_pInt,3_pInt,&
|
||||
3_pInt,1_pInt,&
|
||||
1_pInt,2_pInt],[2,3])
|
||||
integer(pInt) :: i
|
||||
real(pReal), intent(in) :: FWHM
|
||||
real(pReal), dimension(3) :: math_sampleFiberOri, &
|
||||
fInC,& !< fiber axis in crystal coordinate system
|
||||
fInS,& !< fiber axis in sample coordinate system
|
||||
u
|
||||
real(pReal), dimension(3) :: rnd
|
||||
real(pReal), dimension(:),allocatable :: a !< 2D vector to tilt
|
||||
integer(pInt), dimension(:),allocatable :: idx !< components of 2D vector
|
||||
real(pReal), dimension(3,3) :: R !< Rotation matrix (composed of three components)
|
||||
real(pReal):: angle,c
|
||||
integer(pInt):: j,& !< index of smallest component
|
||||
i
|
||||
|
||||
! Helming uses different distribution with Bessel functions
|
||||
! therefore the gauss scatter width has to be scaled differently
|
||||
scatter = 0.95_pReal * noise
|
||||
cos2Scatter = cos(2.0_pReal*scatter)
|
||||
fInC = [sin(alpha(1))*cos(alpha(2)), sin(alpha(1))*sin(alpha(2)), cos(alpha(1))]
|
||||
fInS = [sin(beta(1))*cos(beta(2)), sin(beta(1))*sin(beta(2)), cos(beta(1))]
|
||||
|
||||
! fiber axis in crystal coordinate system
|
||||
fiberInC = [ sin(alpha(1))*cos(alpha(2)) , &
|
||||
sin(alpha(1))*sin(alpha(2)), &
|
||||
cos(alpha(1))]
|
||||
! fiber axis in sample coordinate system
|
||||
fiberInS = [ sin(beta(1))*cos(beta(2)), &
|
||||
sin(beta(1))*sin(beta(2)), &
|
||||
cos(beta(1))]
|
||||
R = math_EulerAxisAngleToR(math_crossproduct(fInC,fInS),-acos(dot_product(fInC,fInS))) !< rotation to align fiber axis in crystal and sample system
|
||||
|
||||
! ---# rotation matrix from sample to crystal system #---
|
||||
angle = -acos(dot_product(fiberInC,fiberInS))
|
||||
if(abs(angle) > tol_math_check) then
|
||||
! rotation axis between sample and crystal system (cross product)
|
||||
forall(i=1_pInt:3_pInt) axis(i) = fiberInC(ROTMAP(1,i))*fiberInS(ROTMAP(2,i))-fiberInC(ROTMAP(2,i))*fiberInS(ROTMAP(1,i))
|
||||
oRot = math_EulerAxisAngleToR(math_crossproduct(fiberInC,fiberInS),angle)
|
||||
rnd = halton([7_pInt,10_pInt,3_pInt])
|
||||
R = math_mul33x33(R,math_EulerAxisAngleToR(fInS,rnd(1)*2.0_pReal*PI)) !< additional rotation (0..360deg) perpendicular to fiber axis
|
||||
|
||||
if (FWHM > 0.1_pReal*INRAD) then
|
||||
reducedTo2D: do i=1_pInt,3_pInt
|
||||
if (i /= minloc(abs(fInS),1)) then
|
||||
a=[a,fInS(i)]
|
||||
idx=[idx,i]
|
||||
else
|
||||
oRot = math_I3
|
||||
j = i
|
||||
endif
|
||||
enddo reducedTo2D
|
||||
GaussConvolution: do
|
||||
angle = (rnd(2)-0.5_pReal)*4.0_pReal*FWHM ! rotation by [0, +-2 FWHM]
|
||||
! solve cos(angle) = dot_product(fInS,u) under the assumption that their smallest component is the same
|
||||
c = cos(angle)-fInS(j)**2
|
||||
u(idx(2)) = -(2.0_pReal*c*a(2) + sqrt(4*((c*a(2))**2-sum(a**2)*(c**2-a(1)**2*(1-fInS(j)**2)))))/&
|
||||
(2*sum(a**2))
|
||||
u(idx(1)) = sqrt(1-u(idx(2))**2-fInS(j)**2)
|
||||
u(j) = fInS(j)
|
||||
|
||||
! ---# rotation matrix about fiber axis (random angle) #---
|
||||
do
|
||||
call halton(6_pInt,rnd)
|
||||
fRot = math_EulerAxisAngleToR(fiberInS,rnd(1)*2.0_pReal*pi)
|
||||
|
||||
! ---# rotation about random axis perpend to fiber #---
|
||||
! random axis pependicular to fiber axis
|
||||
axis(1:2) = rnd(2:3)
|
||||
if (abs(fiberInS(3)) > tol_math_check) then
|
||||
axis(3)=-(axis(1)*fiberInS(1)+axis(2)*fiberInS(2))/fiberInS(3)
|
||||
else if(abs(fiberInS(2)) > tol_math_check) then
|
||||
axis(3)=axis(2)
|
||||
axis(2)=-(axis(1)*fiberInS(1)+axis(3)*fiberInS(3))/fiberInS(2)
|
||||
else if(abs(fiberInS(1)) > tol_math_check) then
|
||||
axis(3)=axis(1)
|
||||
axis(1)=-(axis(2)*fiberInS(2)+axis(3)*fiberInS(3))/fiberInS(1)
|
||||
end if
|
||||
|
||||
! scattered rotation angle
|
||||
if (noise > 0.0_pReal) then
|
||||
angle = acos(cos2Scatter+(1.0_pReal-cos2Scatter)*rnd(4))
|
||||
if (rnd(5) <= exp(-1.0_pReal*(angle/scatter)**2.0_pReal)) exit
|
||||
else
|
||||
angle = 0.0_pReal
|
||||
rejectionSampling: if (rnd(3) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) then
|
||||
R = math_mul33x33(R,math_EulerAxisAngleToR(math_crossproduct(u,fInS),angle)) ! tilt around direction of smallest component
|
||||
exit
|
||||
endif rejectionSampling
|
||||
rnd = halton([7_pInt,10_pInt,3_pInt])
|
||||
enddo GaussConvolution
|
||||
endif
|
||||
enddo
|
||||
if (rnd(6) <= 0.5) angle = -angle
|
||||
|
||||
pRot = math_EulerAxisAngleToR(axis,angle)
|
||||
|
||||
! ---# apply the three rotations #---
|
||||
math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot)))
|
||||
math_sampleFiberOri = math_RtoEuler(R)
|
||||
|
||||
end function math_sampleFiberOri
|
||||
|
||||
|
|
@ -1903,18 +1875,17 @@ real(pReal) function math_sampleGaussVar(meanvalue, stddev, width)
|
|||
|
||||
if (abs(stddev) < tol_math_check) then
|
||||
math_sampleGaussVar = meanvalue
|
||||
return
|
||||
endif
|
||||
|
||||
else
|
||||
myWidth = merge(width,3.0_pReal,present(width)) ! use +-3*sigma as default value for scatter if not given
|
||||
|
||||
do
|
||||
call halton(2_pInt, rnd)
|
||||
rnd = halton([6_pInt,2_pInt])
|
||||
scatter = myWidth * (2.0_pReal * rnd(1) - 1.0_pReal)
|
||||
if (rnd(2) <= exp(-0.5_pReal * scatter ** 2.0_pReal)) exit ! test if scattered value is drawn
|
||||
enddo
|
||||
|
||||
math_sampleGaussVar = scatter * stddev
|
||||
endif
|
||||
|
||||
end function math_sampleGaussVar
|
||||
|
||||
|
|
@ -2285,172 +2256,36 @@ pure function math_invariantsSym33(m)
|
|||
end function math_invariantsSym33
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief computes the next element in the Halton sequence.
|
||||
!> @author John Burkardt
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine halton(ndim, r)
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: ndim !< dimension of the element
|
||||
real(pReal), intent(out), dimension(ndim) :: r !< next element of the current Halton sequence
|
||||
integer(pInt), dimension(ndim) :: base
|
||||
integer(pInt) :: seed
|
||||
integer(pInt), dimension(1) :: value_halton
|
||||
|
||||
call halton_memory ('GET', 'SEED', 1_pInt, value_halton)
|
||||
seed = value_halton(1)
|
||||
|
||||
call halton_memory ('GET', 'BASE', ndim, base)
|
||||
|
||||
call i_to_halton (seed, base, ndim, r)
|
||||
|
||||
value_halton(1) = 1_pInt
|
||||
call halton_memory ('INC', 'SEED', 1_pInt, value_halton)
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
contains
|
||||
|
||||
!-------------------------------------------------------------------------------------------------
|
||||
!> @brief computes an element of a Halton sequence.
|
||||
!> @details Only the absolute value of SEED is considered. SEED = 0 is allowed, and returns R = 0.
|
||||
!> @details Halton Bases should be distinct prime numbers. This routine only checks that each base
|
||||
!> @details is greater than 1.
|
||||
!> @author John Burkardt
|
||||
!> @author Martin Diehl
|
||||
!> @details Incrementally increasing elements of the Halton sequence for given bases (> 0)
|
||||
!> @details Reference:
|
||||
!> @details J.H. Halton: On the efficiency of certain quasi-random sequences of points in evaluating
|
||||
!> @details multi-dimensional integrals, Numerische Mathematik, Volume 2, pages 84-90, 1960.
|
||||
!> @author John Burkardt
|
||||
!-------------------------------------------------------------------------------------------------
|
||||
subroutine i_to_halton (seed, base, ndim, r)
|
||||
use IO, only: &
|
||||
IO_error
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ndim, & !< dimension of the sequence
|
||||
seed !< index of the desired element
|
||||
integer(pInt), intent(in), dimension(ndim) :: base !< Halton bases
|
||||
real(pReal), intent(out), dimension(ndim) :: r !< the SEED-th element of the Halton sequence for the given bases
|
||||
|
||||
real(pReal), dimension(ndim) :: base_inv
|
||||
integer(pInt), dimension(ndim) :: &
|
||||
digit, &
|
||||
seed2
|
||||
|
||||
seed2 = abs(seed)
|
||||
r = 0.0_pReal
|
||||
|
||||
if (any (base(1:ndim) <= 1_pInt)) call IO_error(error_ID=405_pInt)
|
||||
|
||||
base_inv(1:ndim) = 1.0_pReal / real (base(1:ndim), pReal)
|
||||
|
||||
do while ( any ( seed2(1:ndim) /= 0_pInt) )
|
||||
digit(1:ndim) = mod ( seed2(1:ndim), base(1:ndim))
|
||||
r(1:ndim) = r(1:ndim) + real ( digit(1:ndim), pReal) * base_inv(1:ndim)
|
||||
base_inv(1:ndim) = base_inv(1:ndim) / real ( base(1:ndim), pReal)
|
||||
seed2(1:ndim) = seed2(1:ndim) / base(1:ndim)
|
||||
enddo
|
||||
|
||||
end subroutine i_to_halton
|
||||
|
||||
|
||||
end subroutine halton
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief sets or returns quantities associated with the Halton sequence.
|
||||
!> @details If action_halton is 'SET' and action_halton is 'BASE', then NDIM is input, and
|
||||
!> @details is the number of entries in value_halton to be put into BASE.
|
||||
!> @details If action_halton is 'SET', then on input, value_halton contains values to be assigned
|
||||
!> @details to the internal variable.
|
||||
!> @details If action_halton is 'GET', then on output, value_halton contains the values of
|
||||
!> @details the specified internal variable.
|
||||
!> @details If action_halton is 'INC', then on input, value_halton contains the increment to
|
||||
!> @details be added to the specified internal variable.
|
||||
!> @author John Burkardt
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine halton_memory (action_halton, name_halton, ndim, value_halton)
|
||||
use IO, only: &
|
||||
IO_lc
|
||||
|
||||
implicit none
|
||||
character(len = *), intent(in) :: &
|
||||
action_halton, & !< desired action: GET the value of a particular quantity, SET the value of a particular quantity, INC the value of a particular quantity (only for SEED)
|
||||
name_halton !< name of the quantity: BASE: Halton base(s), NDIM: spatial dimension, SEED: current Halton seed
|
||||
integer(pInt), dimension(*), intent(inout) :: value_halton
|
||||
integer(pInt), allocatable, save, dimension(:) :: base
|
||||
logical, save :: first_call = .true.
|
||||
integer(pInt), intent(in) :: ndim !< dimension of the quantity
|
||||
integer(pInt), save :: ndim_save = 0_pInt, seed = 1_pInt
|
||||
integer(pInt) :: i
|
||||
|
||||
if (first_call) then
|
||||
ndim_save = 1_pInt
|
||||
allocate(base(ndim_save))
|
||||
base(1) = 2_pInt
|
||||
first_call = .false.
|
||||
endif
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! Set
|
||||
actionHalton: if(IO_lc(action_halton(1:1)) == 's') then
|
||||
|
||||
nameSet: if(IO_lc(name_halton(1:1)) == 'b') then
|
||||
if(ndim_save /= ndim) ndim_save = ndim
|
||||
base = value_halton(1:ndim)
|
||||
elseif(IO_lc(name_halton(1:1)) == 'n') then nameSet
|
||||
if(ndim_save /= value_halton(1)) then
|
||||
ndim_save = value_halton(1)
|
||||
base = [(prime(i),i=1_pInt,ndim_save)]
|
||||
else
|
||||
ndim_save = value_halton(1)
|
||||
endif
|
||||
elseif(IO_lc(name_halton(1:1)) == 's') then nameSet
|
||||
seed = value_halton(1)
|
||||
endif nameSet
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! Get
|
||||
elseif(IO_lc(action_halton(1:1)) == 'g') then actionHalton
|
||||
nameGet: if(IO_lc(name_halton(1:1)) == 'b') then
|
||||
if(ndim /= ndim_save) then
|
||||
ndim_save = ndim
|
||||
base = [(prime(i),i=1_pInt,ndim_save)]
|
||||
endif
|
||||
value_halton(1:ndim_save) = base(1:ndim_save)
|
||||
elseif(IO_lc(name_halton(1:1)) == 'n') then nameGet
|
||||
value_halton(1) = ndim_save
|
||||
elseif(IO_lc(name_halton(1:1)) == 's') then nameGet
|
||||
value_halton(1) = seed
|
||||
endif nameGet
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! Increment
|
||||
elseif(IO_lc(action_halton(1:1)) == 'i') then actionHalton
|
||||
if(IO_lc(name_halton(1:1)) == 's') seed = seed + value_halton(1)
|
||||
endif actionHalton
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
contains
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns any of the first 1500 prime numbers.
|
||||
!> @details n = 0 is legal, returning PRIME = 1.
|
||||
!> @details Reference:
|
||||
!> @details Reference for prime numbers:
|
||||
!> @details Milton Abramowitz and Irene Stegun: Handbook of Mathematical Functions,
|
||||
!> @details US Department of Commerce, 1964, pages 870-873.
|
||||
!> @details Daniel Zwillinger: CRC Standard Mathematical Tables and Formulae,
|
||||
!> @details 30th Edition, CRC Press, 1996, pages 95-98.
|
||||
!> @author John Burkardt
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
integer(pInt) function prime(n)
|
||||
use IO, only: &
|
||||
IO_error
|
||||
!-------------------------------------------------------------------------------------------------
|
||||
function halton(bases)
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: n !< index of the desired prime number
|
||||
integer(pInt), dimension(0:1500), parameter :: &
|
||||
npvec = int([&
|
||||
integer(pInt), intent(in), dimension(:):: &
|
||||
bases !< bases (prime number ID)
|
||||
real(pReal), dimension(size(bases)) :: &
|
||||
halton
|
||||
integer(pInt), save :: &
|
||||
current = 1_pInt
|
||||
real(pReal), dimension(size(bases)) :: &
|
||||
base_inv
|
||||
integer(pInt), dimension(size(bases)) :: &
|
||||
base, &
|
||||
t
|
||||
integer(pInt), dimension(0:1600), parameter :: &
|
||||
prime = int([&
|
||||
1, &
|
||||
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, &
|
||||
31, 37, 41, 43, 47, 53, 59, 61, 67, 71, &
|
||||
|
|
@ -2597,7 +2432,7 @@ subroutine halton_memory (action_halton, name_halton, ndim, value_halton)
|
|||
! 1301:1400
|
||||
10663, 10667, 10687, 10691, 10709, 10711, 10723, 10729, 10733, 10739, &
|
||||
10753, 10771, 10781, 10789, 10799, 10831, 10837, 10847, 10853, 10859, &
|
||||
10861, 10867, 10883, 10889, 10891, 10903, 10909, 19037, 10939, 10949, &
|
||||
10861, 10867, 10883, 10889, 10891, 10903, 10909, 10937, 10939, 10949, &
|
||||
10957, 10973, 10979, 10987, 10993, 11003, 11027, 11047, 11057, 11059, &
|
||||
11069, 11071, 11083, 11087, 11093, 11113, 11117, 11119, 11131, 11149, &
|
||||
11159, 11161, 11171, 11173, 11177, 11197, 11213, 11239, 11243, 11251, &
|
||||
|
|
@ -2615,58 +2450,34 @@ subroutine halton_memory (action_halton, name_halton, ndim, value_halton)
|
|||
12227, 12239, 12241, 12251, 12253, 12263, 12269, 12277, 12281, 12289, &
|
||||
12301, 12323, 12329, 12343, 12347, 12373, 12377, 12379, 12391, 12401, &
|
||||
12409, 12413, 12421, 12433, 12437, 12451, 12457, 12473, 12479, 12487, &
|
||||
12491, 12497, 12503, 12511, 12517, 12527, 12539, 12541, 12547, 12553],pInt)
|
||||
12491, 12497, 12503, 12511, 12517, 12527, 12539, 12541, 12547, 12553, &
|
||||
! 1501:1600
|
||||
12569, 12577, 12583, 12589, 12601, 12611, 12613, 12619, 12637, 12641, &
|
||||
12647, 12653, 12659, 12671, 12689, 12697, 12703, 12713, 12721, 12739, &
|
||||
12743, 12757, 12763, 12781, 12791, 12799, 12809, 12821, 12823, 12829, &
|
||||
12841, 12853, 12889, 12893, 12899, 12907, 12911, 12917, 12919, 12923, &
|
||||
12941, 12953, 12959, 12967, 12973, 12979, 12983, 13001, 13003, 13007, &
|
||||
13009, 13033, 13037, 13043, 13049, 13063, 13093, 13099, 13103, 13109, &
|
||||
13121, 13127, 13147, 13151, 13159, 13163, 13171, 13177, 13183, 13187, &
|
||||
13217, 13219, 13229, 13241, 13249, 13259, 13267, 13291, 13297, 13309, &
|
||||
13313, 13327, 13331, 13337, 13339, 13367, 13381, 13397, 13399, 13411, &
|
||||
13417, 13421, 13441, 13451, 13457, 13463, 13469, 13477, 13487, 13499],pInt)
|
||||
|
||||
if (n < size(npvec)) then
|
||||
prime = npvec(n)
|
||||
else
|
||||
call IO_error(error_ID=406_pInt)
|
||||
end if
|
||||
current = current + 1_pInt
|
||||
|
||||
end function prime
|
||||
base = prime(bases)
|
||||
base_inv = 1.0_pReal/real(base,pReal)
|
||||
|
||||
end subroutine halton_memory
|
||||
halton = 0.0_pReal
|
||||
t = current
|
||||
|
||||
do while (any( t /= 0_pInt) )
|
||||
halton = halton + real(mod(t,base), pReal) * base_inv
|
||||
base_inv = base_inv / real(base, pReal)
|
||||
t = t / base
|
||||
enddo
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief sets the dimension for a Halton sequence
|
||||
!> @author John Burkardt
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine halton_ndim_set(ndim)
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: ndim !< dimension of the Halton vectors
|
||||
integer(pInt) :: value_halton(1)
|
||||
|
||||
value_halton(1) = ndim
|
||||
call halton_memory ('SET', 'NDIM', 1_pInt, value_halton)
|
||||
|
||||
end subroutine halton_ndim_set
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief sets the seed for the Halton sequence.
|
||||
!> @details Calling HALTON repeatedly returns the elements of the Halton sequence in order,
|
||||
!> @details starting with element number 1.
|
||||
!> @details An internal counter, called SEED, keeps track of the next element to return. Each time
|
||||
!> @details is computed, and then SEED is incremented by 1.
|
||||
!> @details To restart the Halton sequence, it is only necessary to reset SEED to 1. It might also
|
||||
!> @details be desirable to reset SEED to some other value. This routine allows the user to specify
|
||||
!> @details any value of SEED.
|
||||
!> @details The default value of SEED is 1, which restarts the Halton sequence.
|
||||
!> @author John Burkardt
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine halton_seed_set(seed)
|
||||
implicit none
|
||||
|
||||
integer(pInt), parameter :: NDIM = 1_pInt
|
||||
integer(pInt), intent(in) :: seed !< seed for the Halton sequence.
|
||||
integer(pInt) :: value_halton(ndim)
|
||||
|
||||
value_halton(1) = seed
|
||||
call halton_memory ('SET', 'SEED', NDIM, value_halton)
|
||||
|
||||
end subroutine halton_seed_set
|
||||
end function halton
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -2821,4 +2632,135 @@ real(pReal) pure function math_limit(a, left, right)
|
|||
|
||||
end function math_limit
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Modified Bessel I function of order 0
|
||||
!> @author John Burkardt
|
||||
!> @details original version available on https://people.sc.fsu.edu/~jburkardt/f_src/toms715/toms715.html
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
real(pReal) function bessel_i0 (x)
|
||||
use, intrinsic :: IEEE_ARITHMETIC
|
||||
|
||||
implicit none
|
||||
real(pReal), intent(in) :: x
|
||||
integer(pInt) :: i
|
||||
real(pReal) :: sump_p, sump_q, xAbs, xx
|
||||
real(pReal), parameter, dimension(15) :: p_small = real( &
|
||||
[-5.2487866627945699800e-18, -1.5982226675653184646e-14, -2.6843448573468483278e-11, &
|
||||
-3.0517226450451067446e-08, -2.5172644670688975051e-05, -1.5453977791786851041e-02, &
|
||||
-7.0935347449210549190e+00, -2.4125195876041896775e+03, -5.9545626019847898221e+05, &
|
||||
-1.0313066708737980747e+08, -1.1912746104985237192e+10, -8.4925101247114157499e+11, &
|
||||
-3.2940087627407749166e+13, -5.5050369673018427753e+14, -2.2335582639474375249e+15], pReal)
|
||||
real(pReal), parameter, dimension(5) :: q_small = real( &
|
||||
[-3.7277560179962773046e+03, 6.5158506418655165707e+06, -6.5626560740833869295e+09, &
|
||||
3.7604188704092954661e+12, -9.7087946179594019126e+14], pReal)
|
||||
real(pReal), parameter, dimension(8) :: p_large = real( &
|
||||
[-3.9843750000000000000e-01, 2.9205384596336793945e+00, -2.4708469169133954315e+00, &
|
||||
4.7914889422856814203e-01, -3.7384991926068969150e-03, -2.6801520353328635310e-03, &
|
||||
9.9168777670983678974e-05, -2.1877128189032726730e-06], pReal)
|
||||
real(pReal), parameter, dimension(7) :: q_large = real( &
|
||||
[-3.1446690275135491500e+01, 8.5539563258012929600e+01, -6.0228002066743340583e+01, &
|
||||
1.3982595353892851542e+01, -1.1151759188741312645e+00, 3.2547697594819615062e-02, &
|
||||
-5.5194330231005480228e-04], pReal)
|
||||
|
||||
|
||||
xAbs = abs(x)
|
||||
|
||||
argRange: if (xAbs < 5.55e-17_pReal) then
|
||||
bessel_i0 = 1.0_pReal
|
||||
else if (xAbs < 15.0_pReal) then argRange
|
||||
xx = xAbs**2.0_pReal
|
||||
sump_p = p_small(1)
|
||||
do i = 2, 15
|
||||
sump_p = sump_p * xx + p_small(i)
|
||||
end do
|
||||
xx = xx - 225.0_pReal
|
||||
sump_q = ((((xx+q_small(1))*xx+q_small(2))*xx+q_small(3))*xx+q_small(4))*xx+q_small(5)
|
||||
bessel_i0 = sump_p / sump_q
|
||||
else if (xAbs <= 713.986_pReal) then argRange
|
||||
xx = 1.0_pReal / xAbs - 2.0_pReal/30.0_pReal
|
||||
sump_p = ((((((p_large(1)*xx+p_large(2))*xx+p_large(3))*xx+p_large(4))*xx+ &
|
||||
p_large(5))*xx+p_large(6))*xx+p_large(7))*xx+p_large(8)
|
||||
sump_q = ((((((xx+q_large(1))*xx+q_large(2))*xx+q_large(3))*xx+ &
|
||||
q_large(4))*xx+q_large(5))*xx+q_large(6))*xx+q_large(7)
|
||||
bessel_i0 = sump_p / sump_q
|
||||
|
||||
avoidOverflow: if (xAbs > 698.986_pReal) then
|
||||
bessel_i0 = ((bessel_i0*exp(xAbs-40.0_pReal)-p_large(1)*exp(xAbs-40.0_pReal))/sqrt(xAbs))*exp(40.0)
|
||||
else avoidOverflow
|
||||
bessel_i0 = ((bessel_i0*exp(xAbs)-p_large(1)*exp(xAbs))/sqrt(xAbs))
|
||||
endif avoidOverflow
|
||||
|
||||
else argRange
|
||||
bessel_i0 = IEEE_value(bessel_i0,IEEE_positive_inf)
|
||||
end if argRange
|
||||
|
||||
end function bessel_i0
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Modified Bessel I function of order 1
|
||||
!> @author John Burkardt
|
||||
!> @details original version available on https://people.sc.fsu.edu/~jburkardt/f_src/toms715/toms715.html
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
real(pReal) function bessel_i1 (x)
|
||||
use, intrinsic :: IEEE_ARITHMETIC
|
||||
|
||||
implicit none
|
||||
real(pReal), intent(in) :: x
|
||||
integer(pInt) :: i
|
||||
real(pReal) :: sump_p, sump_q, xAbs, xx
|
||||
real(pReal), dimension(15), parameter :: p_small = real( &
|
||||
[-1.9705291802535139930e-19, -6.5245515583151902910e-16, -1.1928788903603238754e-12, &
|
||||
-1.4831904935994647675e-09, -1.3466829827635152875e-06, -9.1746443287817501309e-04, &
|
||||
-4.7207090827310162436e-01, -1.8225946631657315931e+02, -5.1894091982308017540e+04, &
|
||||
-1.0588550724769347106e+07, -1.4828267606612366099e+09, -1.3357437682275493024e+11, &
|
||||
-6.9876779648010090070e+12, -1.7732037840791591320e+14, -1.4577180278143463643e+15], pReal)
|
||||
real(pReal), dimension(5), parameter :: q_small = real( &
|
||||
[-4.0076864679904189921e+03, 7.4810580356655069138e+06, -8.0059518998619764991e+09, &
|
||||
4.8544714258273622913e+12, -1.3218168307321442305e+15], pReal)
|
||||
real(pReal), dimension(8), parameter :: p_large = real( &
|
||||
[-6.0437159056137600000e-02, 4.5748122901933459000e-01, -4.2843766903304806403e-01, &
|
||||
9.7356000150886612134e-02, -3.2457723974465568321e-03, -3.6395264712121795296e-04, &
|
||||
1.6258661867440836395e-05, -3.6347578404608223492e-07], pReal)
|
||||
real(pReal), dimension(6), parameter :: q_large = real( &
|
||||
[-3.8806586721556593450e+00, 3.2593714889036996297e+00, -8.5017476463217924408e-01, &
|
||||
7.4212010813186530069e-02, -2.2835624489492512649e-03, 3.7510433111922824643e-05], pReal)
|
||||
real(pReal), parameter :: pbar = 3.98437500e-01
|
||||
|
||||
|
||||
xAbs = abs(x)
|
||||
|
||||
argRange: if (xAbs < 5.55e-17_pReal) then
|
||||
bessel_i1 = 0.5_pReal * xAbs
|
||||
else if (xAbs < 15.0_pReal) then argRange
|
||||
xx = xAbs**2.0_pReal
|
||||
sump_p = p_small(1)
|
||||
do i = 2, 15
|
||||
sump_p = sump_p * xx + p_small(i)
|
||||
end do
|
||||
xx = xx - 225.0_pReal
|
||||
sump_q = ((((xx+q_small(1))*xx+q_small(2))*xx+q_small(3))*xx+q_small(4)) * xx + q_small(5)
|
||||
bessel_i1 = (sump_p / sump_q) * xAbs
|
||||
else if (xAbs <= 713.986_pReal) then argRange
|
||||
xx = 1.0_pReal / xAbs - 2.0_pReal/30.0_pReal
|
||||
sump_p = ((((((p_large(1)*xx+p_large(2))*xx+p_large(3))*xx+p_large(4))*xx+&
|
||||
p_large(5))*xx+p_large(6))*xx+p_large(7))*xx+p_large(8)
|
||||
sump_q = (((((xx+q_large(1))*xx+q_large(2))*xx+q_large(3))*xx+ q_large(4))*xx+q_large(5))*xx+q_large(6)
|
||||
bessel_i1 = sump_p / sump_q
|
||||
|
||||
avoidOverflow: if (xAbs > 698.986_pReal) then
|
||||
bessel_i1 = ((bessel_i1 * exp(xAbs-40.0_pReal) + pbar * exp(xAbs-40.0_pReal)) / sqrt(xAbs)) * exp(40.0_pReal)
|
||||
else avoidOverflow
|
||||
bessel_i1 = ((bessel_i1 * exp(xAbs) + pbar * exp(xAbs)) / sqrt(xAbs))
|
||||
endif avoidOverflow
|
||||
|
||||
else argRange
|
||||
bessel_i1 = IEEE_value(bessel_i1,IEEE_positive_inf)
|
||||
end if argRange
|
||||
|
||||
if (x < 0.0_pReal) bessel_i1 = -bessel_i1
|
||||
|
||||
end function bessel_i1
|
||||
|
||||
end module math
|
||||
|
|
|
|||
Loading…
Reference in New Issue