From 230ffe8f35c10d8d4e265a3b84a47dac6bb978c8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 14 Feb 2018 13:03:50 +0100 Subject: [PATCH 01/19] Bessel functions for Gaussian scatter --- src/math.f90 | 131 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 131 insertions(+) diff --git a/src/math.f90 b/src/math.f90 index ee228834a..e88a5ed03 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -2821,5 +2821,136 @@ real(pReal) pure function math_limit(a, left, right) math_limit = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_limit, 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 From 9b1e72e7bc4176d0898db90a0887ca138356449b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 15 Feb 2018 11:26:19 +0100 Subject: [PATCH 02/19] more sensible cutoff for scatter --- src/math.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/math.f90 b/src/math.f90 index e88a5ed03..3ede28beb 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1783,7 +1783,7 @@ function math_sampleGaussOri(center,noise) real(pReal), dimension(3), parameter :: ORIGIN = 0.0_pReal real(pReal), dimension(5) :: rnd - noScatter: if (abs(noise) < tol_math_check) then + noScatter: if (noise < 0.5_pReal*INRAD) then math_sampleGaussOri = center else noScatter ! Helming uses different distribution with Bessel functions From 04d873c9c2d2de2812d6acc03578b514caa9e7ab Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 15 Feb 2018 11:26:41 +0100 Subject: [PATCH 03/19] typo caused useless report --- DAMASK_prerequisites.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DAMASK_prerequisites.sh b/DAMASK_prerequisites.sh index 12dd9bd07..9c850aa4d 100755 --- a/DAMASK_prerequisites.sh +++ b/DAMASK_prerequisites.sh @@ -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 } From 6b1c1af609c06db38e3fab5797c1da073609a4a7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Feb 2018 08:55:33 +0100 Subject: [PATCH 04/19] no early return required here --- src/math.f90 | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 3ede28beb..2eda31097 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1902,19 +1902,18 @@ real(pReal) function math_sampleGaussVar(meanvalue, stddev, width) if (abs(stddev) < tol_math_check) then math_sampleGaussVar = meanvalue - return + 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) + 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 - 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) - 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 - end function math_sampleGaussVar From a52f54a9a00a8b5996fb3dac61d62e372ec960b6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Feb 2018 14:37:02 +0100 Subject: [PATCH 05/19] fixed typo in prime number list and extended to 1600 values source https://people.sc.fsu.edu/~jburkardt/f_src/halton/halton.f90 --- src/math.f90 | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 2eda31097..b6e3658e9 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -2447,7 +2447,7 @@ subroutine halton_memory (action_halton, name_halton, ndim, value_halton) implicit none integer(pInt), intent(in) :: n !< index of the desired prime number - integer(pInt), dimension(0:1500), parameter :: & + integer(pInt), dimension(0:1600), parameter :: & npvec = int([& 1, & 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, & @@ -2595,7 +2595,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, & @@ -2613,8 +2613,19 @@ 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 From 98df2d14278f592008d761748570d87753b6cc80 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Feb 2018 16:49:29 +0100 Subject: [PATCH 06/19] better description and names --- src/math.f90 | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index b6e3658e9..2fab73955 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1769,26 +1769,27 @@ 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) !-------------------------------------------------------------------------------------------------- -function math_sampleGaussOri(center,noise) +function math_sampleGaussOri(center,FWHM) use prec, only: & tol_math_check 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 - noScatter: if (noise < 0.5_pReal*INRAD) then + noScatter: if (FWHM < 0.5_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 + scatter = 0.95_pReal * FWHM cosScatter = cos(scatter) do @@ -1807,18 +1808,20 @@ end function math_sampleGaussOri !-------------------------------------------------------------------------------------------------- -!> @brief draw a random sample from Fiber component with noise (in radians) -!-------------------------------------------------------------------------------------------------- -function math_sampleFiberOri(alpha,beta,noise) +!> @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) use prec, only: & tol_math_check implicit none real(pReal), dimension(3) :: math_sampleFiberOri, fiberInC,fiberInS,axis real(pReal), dimension(2), intent(in) :: alpha,beta + real(pReal), intent(in) :: FWHM real(pReal), dimension(6) :: rnd real(pReal), dimension(3,3) :: oRot,fRot,pRot - real(pReal) :: noise, scatter, cos2Scatter, angle + real(pReal) :: 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]) @@ -1826,7 +1829,7 @@ function math_sampleFiberOri(alpha,beta,noise) ! Helming uses different distribution with Bessel functions ! therefore the gauss scatter width has to be scaled differently - scatter = 0.95_pReal * noise + scatter = 0.95_pReal * FWHM cos2Scatter = cos(2.0_pReal*scatter) ! fiber axis in crystal coordinate system @@ -1867,7 +1870,7 @@ function math_sampleFiberOri(alpha,beta,noise) end if ! scattered rotation angle - if (noise > 0.0_pReal) then + if (FWHM > 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 From ae27660e86827939c1bb6074bda302a0e29bcdb9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Feb 2018 18:47:39 +0100 Subject: [PATCH 07/19] simplified halto procedure (still needs testing) --- src/math.f90 | 620 +++++++++++++++++++-------------------------------- 1 file changed, 232 insertions(+), 388 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 2fab73955..09e0b36ba 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -161,10 +161,8 @@ module math math_rotate_forward3333, & math_limit private :: & - halton, & - halton_memory, & - halton_ndim_set, & - halton_seed_set + math_check, & + halton contains @@ -217,8 +215,7 @@ subroutine math_init call random_seed(put = randInit) - call halton_seed_set(int(randInit(1), pInt)) - call halton_ndim_set(3_pInt) + randTest = halton(4,seed = randTest(1)) call math_check() @@ -283,7 +280,7 @@ subroutine math_check endif ! +++ check rotation sense of q and R +++ - call halton(3_pInt,v) ! random vector + v = halton(3_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' @@ -1237,7 +1234,7 @@ function math_qRand() real(pReal), dimension(4) :: math_qRand real(pReal), dimension(3) :: rnd - call halton(3_pInt,rnd) + rnd = halton(3_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)), & @@ -1760,7 +1757,7 @@ function math_sampleRandomOri() implicit none real(pReal), dimension(3) :: math_sampleRandomOri, rnd - call halton(3_pInt,rnd) + rnd = halton(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] @@ -1793,7 +1790,7 @@ function math_sampleGaussOri(center,FWHM) cosScatter = cos(scatter) do - call halton(5_pInt,rnd) + rnd = halton(5_pInt) 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 @@ -1853,7 +1850,7 @@ function math_sampleFiberOri(alpha,beta,FWHM) ! ---# rotation matrix about fiber axis (random angle) #--- do - call halton(6_pInt,rnd) + rnd = halton(6_pInt) fRot = math_EulerAxisAngleToR(fiberInS,rnd(1)*2.0_pReal*pi) ! ---# rotation about random axis perpend to fiber #--- @@ -1909,7 +1906,7 @@ real(pReal) function math_sampleGaussVar(meanvalue, stddev, width) 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(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 @@ -2286,157 +2283,238 @@ pure function math_invariantsSym33(m) end function math_invariantsSym33 -!-------------------------------------------------------------------------------------------------- -!> @brief computes the next element in the Halton sequence. +!------------------------------------------------------------------------------------------------- +!> @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. +!> @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 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. - !> @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 +!------------------------------------------------------------------------------------------------- +function halton(ndim, seed) - 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 + integer(pInt), intent(in) :: & + ndim !< dimension of the sequence + real(pReal), intent(in), optional :: & + seed !< seed value [0.0,1.0] + real(pReal), dimension(ndim) :: & + halton + integer(pInt), save :: & + start, & + current + real(pReal), dimension(ndim) :: & + base_inv + integer(pInt), dimension(ndim) :: & + base, & + t + integer(pInt) :: & + i + 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, & + 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, & + 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, & + 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, & + 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, & + 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, & + 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, & + 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, & + 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, & + ! 101:200 + 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, & + 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, & + 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, & + 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, & + 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, & + 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, & + 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, & + 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, & + 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, & + 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, & + ! 201:300 + 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, & + 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, & + 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, & + 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, & + 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, & + 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, & + 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, & + 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, & + 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, & + 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, & + ! 301:400 + 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, & + 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, & + 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, & + 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, & + 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, & + 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, & + 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, & + 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, & + 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, & + 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, & + ! 401:500 + 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, & + 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, & + 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, & + 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, & + 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, & + 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, & + 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, & + 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, & + 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, & + 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, & + ! 501:600 + 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, & + 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, & + 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, & + 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, & + 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, & + 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, & + 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, & + 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, & + 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, & + 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, & + ! 601:700 + 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, & + 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, & + 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, & + 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, & + 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, & + 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, & + 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, & + 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, & + 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, & + 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, & + ! 701:800 + 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, & + 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, & + 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, & + 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639, & + 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, & + 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, & + 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, & + 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, & + 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, & + 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, & + ! 801:900 + 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, & + 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, & + 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, & + 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, & + 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, & + 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, & + 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, & + 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, & + 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, & + 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, & + ! 901:1000 + 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, & + 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, & + 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, & + 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, & + 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, & + 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, & + 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, & + 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, & + 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, & + 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919, & + ! 1001:1100 + 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017, & + 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, & + 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219, & + 8221, 8231, 8233, 8237, 8243, 8263, 8269, 8273, 8287, 8291, & + 8293, 8297, 8311, 8317, 8329, 8353, 8363, 8369, 8377, 8387, & + 8389, 8419, 8423, 8429, 8431, 8443, 8447, 8461, 8467, 8501, & + 8513, 8521, 8527, 8537, 8539, 8543, 8563, 8573, 8581, 8597, & + 8599, 8609, 8623, 8627, 8629, 8641, 8647, 8663, 8669, 8677, & + 8681, 8689, 8693, 8699, 8707, 8713, 8719, 8731, 8737, 8741, & + 8747, 8753, 8761, 8779, 8783, 8803, 8807, 8819, 8821, 8831, & + ! 1101:1200 + 8837, 8839, 8849, 8861, 8863, 8867, 8887, 8893, 8923, 8929, & + 8933, 8941, 8951, 8963, 8969, 8971, 8999, 9001, 9007, 9011, & + 9013, 9029, 9041, 9043, 9049, 9059, 9067, 9091, 9103, 9109, & + 9127, 9133, 9137, 9151, 9157, 9161, 9173, 9181, 9187, 9199, & + 9203, 9209, 9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283, & + 9293, 9311, 9319, 9323, 9337, 9341, 9343, 9349, 9371, 9377, & + 9391, 9397, 9403, 9413, 9419, 9421, 9431, 9433, 9437, 9439, & + 9461, 9463, 9467, 9473, 9479, 9491, 9497, 9511, 9521, 9533, & + 9539, 9547, 9551, 9587, 9601, 9613, 9619, 9623, 9629, 9631, & + 9643, 9649, 9661, 9677, 9679, 9689, 9697, 9719, 9721, 9733, & + ! 1201:1300 + 9739, 9743, 9749, 9767, 9769, 9781, 9787, 9791, 9803, 9811, & + 9817, 9829, 9833, 9839, 9851, 9857, 9859, 9871, 9883, 9887, & + 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973, 10007, & + 10009, 10037, 10039, 10061, 10067, 10069, 10079, 10091, 10093, 10099, & + 10103, 10111, 10133, 10139, 10141, 10151, 10159, 10163, 10169, 10177, & + 10181, 10193, 10211, 10223, 10243, 10247, 10253, 10259, 10267, 10271, & + 10273, 10289, 10301, 10303, 10313, 10321, 10331, 10333, 10337, 10343, & + 10357, 10369, 10391, 10399, 10427, 10429, 10433, 10453, 10457, 10459, & + 10463, 10477, 10487, 10499, 10501, 10513, 10529, 10531, 10559, 10567, & + 10589, 10597, 10601, 10607, 10613, 10627, 10631, 10639, 10651, 10657, & + ! 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, 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, & + 11257, 11261, 11273, 11279, 11287, 11299, 11311, 11317, 11321, 11329, & + 11351, 11353, 11369, 11383, 11393, 11399, 11411, 11423, 11437, 11443, & + 11447, 11467, 11471, 11483, 11489, 11491, 11497, 11503, 11519, 11527, & + 11549, 11551, 11579, 11587, 11593, 11597, 11617, 11621, 11633, 11657, & + ! 1401:1500 + 11677, 11681, 11689, 11699, 11701, 11717, 11719, 11731, 11743, 11777, & + 11779, 11783, 11789, 11801, 11807, 11813, 11821, 11827, 11831, 11833, & + 11839, 11863, 11867, 11887, 11897, 11903, 11909, 11923, 11927, 11933, & + 11939, 11941, 11953, 11959, 11969, 11971, 11981, 11987, 12007, 12011, & + 12037, 12041, 12043, 12049, 12071, 12073, 12097, 12101, 12107, 12109, & + 12113, 12119, 12143, 12149, 12157, 12161, 12163, 12197, 12203, 12211, & + 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, & + ! 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 (first_call) then - ndim_save = 1_pInt - allocate(base(ndim_save)) - base(1) = 2_pInt - first_call = .false. + + if (present(seed)) then + start = int(seed * 1599.0_pReal,pInt) + 1_pInt ! select one prime number to start with + current = 1_pInt + else + current = current + 1_pInt 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 + base = [(prime(mod(start,1600_pInt)+i), i=1_pInt, ndim)] + base_inv = 1.0_pReal/real(base,pReal) -!-------------------------------------------------------------------------------------------------- - contains + 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 returns any of the first 1500 prime numbers. + !> @brief returns any of the first 1600 prime numbers. !> @details n = 0 is legal, returning PRIME = 1. + !> @details 0 > n > 1600 returns -1 !> @details Reference: !> @details Milton Abramowitz and Irene Stegun: Handbook of Mathematical Functions, !> @details US Department of Commerce, 1964, pages 870-873. @@ -2444,241 +2522,7 @@ subroutine halton_memory (action_halton, name_halton, ndim, value_halton) !> @details 30th Edition, CRC Press, 1996, pages 95-98. !> @author John Burkardt !-------------------------------------------------------------------------------------------------- - integer(pInt) function prime(n) - use IO, only: & - IO_error - - implicit none - integer(pInt), intent(in) :: n !< index of the desired prime number - integer(pInt), dimension(0:1600), parameter :: & - npvec = int([& - 1, & - 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, & - 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, & - 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, & - 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, & - 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, & - 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, & - 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, & - 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, & - 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, & - 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, & - ! 101:200 - 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, & - 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, & - 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, & - 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, & - 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, & - 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, & - 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, & - 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, & - 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, & - 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, & - ! 201:300 - 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, & - 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, & - 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, & - 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, & - 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, & - 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, & - 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, & - 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, & - 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, & - 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, & - ! 301:400 - 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, & - 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, & - 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, & - 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, & - 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, & - 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, & - 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, & - 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, & - 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, & - 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, & - ! 401:500 - 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, & - 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, & - 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, & - 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, & - 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, & - 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, & - 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, & - 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, & - 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, & - 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, & - ! 501:600 - 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, & - 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, & - 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, & - 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, & - 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, & - 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, & - 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, & - 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, & - 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, & - 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, & - ! 601:700 - 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, & - 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, & - 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, & - 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, & - 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, & - 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, & - 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, & - 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, & - 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, & - 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, & - ! 701:800 - 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, & - 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, & - 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, & - 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639, & - 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, & - 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, & - 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, & - 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, & - 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, & - 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, & - ! 801:900 - 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, & - 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, & - 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, & - 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, & - 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, & - 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, & - 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, & - 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, & - 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, & - 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, & - ! 901:1000 - 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, & - 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, & - 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, & - 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, & - 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, & - 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, & - 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, & - 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, & - 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, & - 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919, & - ! 1001:1100 - 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017, & - 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, & - 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219, & - 8221, 8231, 8233, 8237, 8243, 8263, 8269, 8273, 8287, 8291, & - 8293, 8297, 8311, 8317, 8329, 8353, 8363, 8369, 8377, 8387, & - 8389, 8419, 8423, 8429, 8431, 8443, 8447, 8461, 8467, 8501, & - 8513, 8521, 8527, 8537, 8539, 8543, 8563, 8573, 8581, 8597, & - 8599, 8609, 8623, 8627, 8629, 8641, 8647, 8663, 8669, 8677, & - 8681, 8689, 8693, 8699, 8707, 8713, 8719, 8731, 8737, 8741, & - 8747, 8753, 8761, 8779, 8783, 8803, 8807, 8819, 8821, 8831, & - ! 1101:1200 - 8837, 8839, 8849, 8861, 8863, 8867, 8887, 8893, 8923, 8929, & - 8933, 8941, 8951, 8963, 8969, 8971, 8999, 9001, 9007, 9011, & - 9013, 9029, 9041, 9043, 9049, 9059, 9067, 9091, 9103, 9109, & - 9127, 9133, 9137, 9151, 9157, 9161, 9173, 9181, 9187, 9199, & - 9203, 9209, 9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283, & - 9293, 9311, 9319, 9323, 9337, 9341, 9343, 9349, 9371, 9377, & - 9391, 9397, 9403, 9413, 9419, 9421, 9431, 9433, 9437, 9439, & - 9461, 9463, 9467, 9473, 9479, 9491, 9497, 9511, 9521, 9533, & - 9539, 9547, 9551, 9587, 9601, 9613, 9619, 9623, 9629, 9631, & - 9643, 9649, 9661, 9677, 9679, 9689, 9697, 9719, 9721, 9733, & - ! 1201:1300 - 9739, 9743, 9749, 9767, 9769, 9781, 9787, 9791, 9803, 9811, & - 9817, 9829, 9833, 9839, 9851, 9857, 9859, 9871, 9883, 9887, & - 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973, 10007, & - 10009, 10037, 10039, 10061, 10067, 10069, 10079, 10091, 10093, 10099, & - 10103, 10111, 10133, 10139, 10141, 10151, 10159, 10163, 10169, 10177, & - 10181, 10193, 10211, 10223, 10243, 10247, 10253, 10259, 10267, 10271, & - 10273, 10289, 10301, 10303, 10313, 10321, 10331, 10333, 10337, 10343, & - 10357, 10369, 10391, 10399, 10427, 10429, 10433, 10453, 10457, 10459, & - 10463, 10477, 10487, 10499, 10501, 10513, 10529, 10531, 10559, 10567, & - 10589, 10597, 10601, 10607, 10613, 10627, 10631, 10639, 10651, 10657, & - ! 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, 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, & - 11257, 11261, 11273, 11279, 11287, 11299, 11311, 11317, 11321, 11329, & - 11351, 11353, 11369, 11383, 11393, 11399, 11411, 11423, 11437, 11443, & - 11447, 11467, 11471, 11483, 11489, 11491, 11497, 11503, 11519, 11527, & - 11549, 11551, 11579, 11587, 11593, 11597, 11617, 11621, 11633, 11657, & - ! 1401:1500 - 11677, 11681, 11689, 11699, 11701, 11717, 11719, 11731, 11743, 11777, & - 11779, 11783, 11789, 11801, 11807, 11813, 11821, 11827, 11831, 11833, & - 11839, 11863, 11867, 11887, 11897, 11903, 11909, 11923, 11927, 11933, & - 11939, 11941, 11953, 11959, 11969, 11971, 11981, 11987, 12007, 12011, & - 12037, 12041, 12043, 12049, 12071, 12073, 12097, 12101, 12107, 12109, & - 12113, 12119, 12143, 12149, 12157, 12161, 12163, 12197, 12203, 12211, & - 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, & - ! 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 - - end function prime - -end subroutine halton_memory - - -!-------------------------------------------------------------------------------------------------- -!> @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 !-------------------------------------------------------------------------------------------------- From 9173d12d149679dd43d767e01729f34adbaeb5b3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Feb 2018 20:32:52 +0100 Subject: [PATCH 08/19] correct algorithm for sampling of uniform orientations and fix for Halton series Halton series gives strange results for large prime numbers, now always starting with 2 for first dimension, 3 for second etc. Consecutive Halton numbers for rejection sampling seem to cause problems (i.e. introduce patterns). Algorithm for uniformly distributed orientations with FWHM specified is taken from https://math.stackexchange.com/questions/131336. WIP: Gauss filtering is currently not implemented! --- src/math.f90 | 61 ++++++++++++++++++++++++++-------------------------- 1 file changed, 30 insertions(+), 31 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 09e0b36ba..fd4a8c665 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -214,9 +214,6 @@ subroutine math_init write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest call random_seed(put = randInit) - - randTest = halton(4,seed = randTest(1)) - call math_check() end subroutine math_init @@ -1776,34 +1773,45 @@ function math_sampleGaussOri(center,FWHM) implicit none 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) :: scatter, omega + real(pReal), dimension(3) :: math_sampleGaussOri, axis + real(pReal), dimension(2) :: rnd + real(pReal), dimension(3,3) :: R noScatter: if (FWHM < 0.5_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 + ! therefore the gauss scatter width has to be scaled differentlyuu scatter = 0.95_pReal * FWHM - cosScatter = cos(scatter) - do - rnd = halton(5_pInt) - 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 + GaussConvolution: do + selectiveSampling: if (FWHM*INRAD < 90.0_pReal) then + rnd = halton(2_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)] + do + call random_number(rnd) ! no halton to avoid correlation + omega = (rnd(1)*(2.0_pReal*FWHM)**3.0_pReal)**(1.0_pReal/3.0_pReal) + if (rnd(2) < sin(omega)**2.0_pReal/omega**2.0_pReal) exit + enddo + R = math_axisAngleToR(axis,omega) + else selectiveSampling + R = math_EulerToR(math_sampleRandomOri()) + endif selectiveSampling + exit + !if (rnd(3) <= exp(-(math_EulerMisorientation([0.0_pReal,0.0_pReal,0.0_pReal],& + ! math_RtoEuler(R))/scatter)**2_pReal)) exit + enddo GaussConvolution - math_sampleGaussOri = math_RtoEuler(math_mul33x33(math_EulerToR(disturb),math_EulerToR(center))) + math_sampleGaussOri = math_RtoEuler(math_mul33x33(R,math_EulerToR(center))) endif noScatter end function math_sampleGaussOri + !-------------------------------------------------------------------------------------------------- !> @brief draw a sample from an Gaussian distribution around given fiber texture and Full Width ! at Half Maximum (FWHM) @@ -2293,18 +2301,15 @@ end function math_invariantsSym33 !> @details multi-dimensional integrals, Numerische Mathematik, Volume 2, pages 84-90, 1960. !> @author John Burkardt !------------------------------------------------------------------------------------------------- -function halton(ndim, seed) +function halton(ndim) implicit none integer(pInt), intent(in) :: & ndim !< dimension of the sequence - real(pReal), intent(in), optional :: & - seed !< seed value [0.0,1.0] real(pReal), dimension(ndim) :: & halton integer(pInt), save :: & - start, & - current + current = 1_pInt real(pReal), dimension(ndim) :: & base_inv integer(pInt), dimension(ndim) :: & @@ -2491,15 +2496,9 @@ function halton(ndim, seed) 13313, 13327, 13331, 13337, 13339, 13367, 13381, 13397, 13399, 13411, & 13417, 13421, 13441, 13451, 13457, 13463, 13469, 13477, 13487, 13499],pInt) - - if (present(seed)) then - start = int(seed * 1599.0_pReal,pInt) + 1_pInt ! select one prime number to start with - current = 1_pInt - else - current = current + 1_pInt - endif + current = current + 1_pInt - base = [(prime(mod(start,1600_pInt)+i), i=1_pInt, ndim)] + base = [(prime(i), i=1_pInt, ndim)] base_inv = 1.0_pReal/real(base,pReal) halton = 0.0_pReal From 5c908e44ecd0d22fc0af314686fc2f525a457abd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Feb 2018 00:10:38 +0100 Subject: [PATCH 09/19] correct scaling for FWHM FWHM was wrong by a factor of approx 2. Analytic expression can be found in https://en.wikipedia.org/wiki/Gaussian_function --- src/math.f90 | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index fd4a8c665..4e40d0e43 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1773,7 +1773,7 @@ function math_sampleGaussOri(center,FWHM) implicit none real(pReal), intent(in) :: FWHM real(pReal), dimension(3), intent(in) :: center - real(pReal) :: scatter, omega + real(pReal) :: omega real(pReal), dimension(3) :: math_sampleGaussOri, axis real(pReal), dimension(2) :: rnd real(pReal), dimension(3,3) :: R @@ -1781,10 +1781,6 @@ function math_sampleGaussOri(center,FWHM) noScatter: if (FWHM < 0.5_pReal*INRAD) then math_sampleGaussOri = center else noScatter - ! Helming uses different distribution with Bessel functions - ! therefore the gauss scatter width has to be scaled differentlyuu - scatter = 0.95_pReal * FWHM - GaussConvolution: do selectiveSampling: if (FWHM*INRAD < 90.0_pReal) then rnd = halton(2_pInt) @@ -1800,9 +1796,9 @@ function math_sampleGaussOri(center,FWHM) else selectiveSampling R = math_EulerToR(math_sampleRandomOri()) endif selectiveSampling - exit - !if (rnd(3) <= exp(-(math_EulerMisorientation([0.0_pReal,0.0_pReal,0.0_pReal],& - ! math_RtoEuler(R))/scatter)**2_pReal)) exit + call random_number(rnd) ! no halton to avoid correlation + omega = math_EulerMisorientation([0.0_pReal,0.0_pReal,0.0_pReal],math_RtoEuler(R)) + if (rnd(1) <= exp(-4.0_pReal*log(2.0_pReal)*(omega/FWHM)**2_pReal)) exit enddo GaussConvolution math_sampleGaussOri = math_RtoEuler(math_mul33x33(R,math_EulerToR(center))) From c6c66bb653b5af95919d24332b9b8e2ffa0f8aab Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Feb 2018 00:42:06 +0100 Subject: [PATCH 10/19] using proper Gauss sampling also for the fiber components function most probably still contains a bu --- src/math.f90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 4e40d0e43..4d2a124d2 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1828,8 +1828,8 @@ function math_sampleFiberOri(alpha,beta,FWHM) 1_pInt,2_pInt],[2,3]) integer(pInt) :: i -! Helming uses different distribution with Bessel functions -! therefore the gauss scatter width has to be scaled differently +! MD: this is a leftover from using the wrongly scaled Gaussian distribution. +! I'm relatively sure that it is not correct to scale by any constant factor here scatter = 0.95_pReal * FWHM cos2Scatter = cos(2.0_pReal*scatter) @@ -1853,9 +1853,9 @@ function math_sampleFiberOri(alpha,beta,FWHM) end if ! ---# rotation matrix about fiber axis (random angle) #--- - do + GaussConvolution: do rnd = halton(6_pInt) - fRot = math_EulerAxisAngleToR(fiberInS,rnd(1)*2.0_pReal*pi) + fRot = math_EulerAxisAngleToR(fiberInS,rnd(1)*2.0_pReal*PI) ! ---# rotation about random axis perpend to fiber #--- ! random axis pependicular to fiber axis @@ -1872,13 +1872,13 @@ function math_sampleFiberOri(alpha,beta,FWHM) ! scattered rotation angle if (FWHM > 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 + angle = acos(cos2Scatter+(1.0_pReal-cos2Scatter)*rnd(4)) ! MD: This is probably not ok + if (rnd(5) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2.0_pReal)) exit else angle = 0.0_pReal exit end if - enddo + enddo GaussConvolution if (rnd(6) <= 0.5) angle = -angle pRot = math_EulerAxisAngleToR(axis,angle) From e51de7ffd8566b4bb910bc525837a0bd8e04c40a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Feb 2018 14:16:36 +0100 Subject: [PATCH 11/19] explicitly select halton bases in call --- src/math.f90 | 46 +++++++++++++++++++++++----------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 4d2a124d2..a4fead8e0 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -277,7 +277,7 @@ subroutine math_check endif ! +++ check rotation sense of q and R +++ - v = halton(3_pInt) ! 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' @@ -1231,7 +1231,7 @@ function math_qRand() real(pReal), dimension(4) :: math_qRand real(pReal), dimension(3) :: rnd - rnd = halton(3_pInt) + 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)), & @@ -1754,7 +1754,7 @@ function math_sampleRandomOri() implicit none real(pReal), dimension(3) :: math_sampleRandomOri, rnd - rnd = halton(3_pInt) + 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] @@ -1783,12 +1783,12 @@ function math_sampleGaussOri(center,FWHM) else noScatter GaussConvolution: do selectiveSampling: if (FWHM*INRAD < 90.0_pReal) then - rnd = halton(2_pInt) + rnd = halton([3_pInt,6_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)] do - call random_number(rnd) ! no halton to avoid correlation + rnd = halton([14_pInt,10_pInt]) omega = (rnd(1)*(2.0_pReal*FWHM)**3.0_pReal)**(1.0_pReal/3.0_pReal) if (rnd(2) < sin(omega)**2.0_pReal/omega**2.0_pReal) exit enddo @@ -1796,7 +1796,7 @@ function math_sampleGaussOri(center,FWHM) else selectiveSampling R = math_EulerToR(math_sampleRandomOri()) endif selectiveSampling - call random_number(rnd) ! no halton to avoid correlation + rnd = halton([8_pInt,11_pInt]) omega = math_EulerMisorientation([0.0_pReal,0.0_pReal,0.0_pReal],math_RtoEuler(R)) if (rnd(1) <= exp(-4.0_pReal*log(2.0_pReal)*(omega/FWHM)**2_pReal)) exit enddo GaussConvolution @@ -1822,7 +1822,7 @@ function math_sampleFiberOri(alpha,beta,FWHM) real(pReal), intent(in) :: FWHM real(pReal), dimension(6) :: rnd real(pReal), dimension(3,3) :: oRot,fRot,pRot - real(pReal) :: scatter, cos2Scatter, angle + real(pReal) :: cos2Scatter, angle integer(pInt), dimension(2,3), parameter :: ROTMAP = reshape([2_pInt,3_pInt,& 3_pInt,1_pInt,& 1_pInt,2_pInt],[2,3]) @@ -1830,8 +1830,7 @@ function math_sampleFiberOri(alpha,beta,FWHM) ! MD: this is a leftover from using the wrongly scaled Gaussian distribution. ! I'm relatively sure that it is not correct to scale by any constant factor here - scatter = 0.95_pReal * FWHM - cos2Scatter = cos(2.0_pReal*scatter) + cos2Scatter = cos(2.0_pReal*0.95_pReal * FWHM) ! fiber axis in crystal coordinate system fiberInC = [ sin(alpha(1))*cos(alpha(2)) , & @@ -1852,13 +1851,13 @@ function math_sampleFiberOri(alpha,beta,FWHM) oRot = math_I3 end if -! ---# rotation matrix about fiber axis (random angle) #--- GaussConvolution: do - rnd = halton(6_pInt) + rnd = halton(int([14,5,10,3,9,17],pInt)) + + ! random rotation around fiber axis fRot = math_EulerAxisAngleToR(fiberInS,rnd(1)*2.0_pReal*PI) -! ---# rotation about random axis perpend to fiber #--- -! random axis pependicular to fiber axis + ! 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) @@ -1870,9 +1869,10 @@ function math_sampleFiberOri(alpha,beta,FWHM) axis(1)=-(axis(2)*fiberInS(2)+axis(3)*fiberInS(3))/fiberInS(1) end if -! scattered rotation angle + ! rotate by XXX amount around axis perpendicular to fiber axis + ! Proabably https://math.stackexchange.com/questions/56784 is the right approach if (FWHM > 0.0_pReal) then - angle = acos(cos2Scatter+(1.0_pReal-cos2Scatter)*rnd(4)) ! MD: This is probably not ok + angle = acos(cos2Scatter+(1.0_pReal-cos2Scatter)*rnd(4)) if (rnd(5) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2.0_pReal)) exit else angle = 0.0_pReal @@ -1910,7 +1910,7 @@ real(pReal) function math_sampleGaussVar(meanvalue, stddev, width) myWidth = merge(width,3.0_pReal,present(width)) ! use +-3*sigma as default value for scatter if not given do - rnd = halton(2_pInt) + 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 @@ -2297,18 +2297,18 @@ end function math_invariantsSym33 !> @details multi-dimensional integrals, Numerische Mathematik, Volume 2, pages 84-90, 1960. !> @author John Burkardt !------------------------------------------------------------------------------------------------- -function halton(ndim) +function halton(bases) implicit none - integer(pInt), intent(in) :: & - ndim !< dimension of the sequence - real(pReal), dimension(ndim) :: & + integer(pInt), intent(in), dimension(:):: & + bases !< dimension of the sequence + real(pReal), dimension(size(bases)) :: & halton integer(pInt), save :: & current = 1_pInt - real(pReal), dimension(ndim) :: & + real(pReal), dimension(size(bases)) :: & base_inv - integer(pInt), dimension(ndim) :: & + integer(pInt), dimension(size(bases)) :: & base, & t integer(pInt) :: & @@ -2494,7 +2494,7 @@ function halton(ndim) current = current + 1_pInt - base = [(prime(i), i=1_pInt, ndim)] + base = prime(bases) base_inv = 1.0_pReal/real(base,pReal) halton = 0.0_pReal From 4f4fa5daf8f61d62cc62a9d58b0519e5c4c55ea0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Feb 2018 14:54:49 +0100 Subject: [PATCH 12/19] simplified --- src/math.f90 | 24 +++++++----------------- 1 file changed, 7 insertions(+), 17 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index a4fead8e0..32e96c4be 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1817,16 +1817,12 @@ function math_sampleFiberOri(alpha,beta,FWHM) tol_math_check implicit none - real(pReal), dimension(3) :: math_sampleFiberOri, fiberInC,fiberInS,axis real(pReal), dimension(2), intent(in) :: alpha,beta real(pReal), intent(in) :: FWHM + real(pReal), dimension(3) :: math_sampleFiberOri, fiberInC,fiberInS,axis real(pReal), dimension(6) :: rnd real(pReal), dimension(3,3) :: oRot,fRot,pRot real(pReal) :: 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 ! MD: this is a leftover from using the wrongly scaled Gaussian distribution. ! I'm relatively sure that it is not correct to scale by any constant factor here @@ -1841,15 +1837,8 @@ function math_sampleFiberOri(alpha,beta,FWHM) sin(beta(1))*sin(beta(2)), & cos(beta(1))] -! ---# 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) - else - oRot = math_I3 - end if + oRot = merge(math_EulerAxisAngleToR(math_crossproduct(fiberInC,fiberInS),angle), + math_I3, acos(dot_product(fiberInC,fiberInS)) > tol_math_check) GaussConvolution: do rnd = halton(int([14,5,10,3,9,17],pInt)) @@ -1860,14 +1849,15 @@ function math_sampleFiberOri(alpha,beta,FWHM) ! 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) + axis(3)=-(axis(1)*fiberInS(1)+axis(2)*fiberInS(2)) 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) + axis(2)=-(axis(1)*fiberInS(1)+axis(3)*fiberInS(3)) 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) + axis(1)=-(axis(2)*fiberInS(2)+axis(3)*fiberInS(3)) end if + axis = axis/norm2(axis) ! rotate by XXX amount around axis perpendicular to fiber axis ! Proabably https://math.stackexchange.com/questions/56784 is the right approach From 0ccce7faccab175deb7b79629b26f5445707503e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 23 Feb 2018 11:33:21 +0100 Subject: [PATCH 13/19] fixed missing angle initialization, simplified and commented --- src/math.f90 | 142 +++++++++++++++++++++------------------------------ 1 file changed, 59 insertions(+), 83 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 32e96c4be..80c678164 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1519,7 +1519,7 @@ pure function math_axisAngleToR(axis,omega) norm = norm2(axis) wellDefined: if (norm > 1.0e-8_pReal) then - n = axis/norm ! normalize axis to be sure + n = axis/norm ! normalize axis to be sure s = sin(omega) c = cos(omega) @@ -1765,20 +1765,26 @@ end function math_sampleRandomOri !-------------------------------------------------------------------------------------------------- !> @brief draw a sample from an Gaussian distribution around given orientation and Full Width ! at Half Maximum (FWHM) +!> @details: for very small FWHM values the given orientation is returned. +!> @details: for intermediate FWHM values, an orientation is picked from uniformly distributed +!> @details: oreintations around the nominal orientation with maximum misorientation of 2*FWHM +!> @deatils: according to https://math.stackexchange.com/questions/13133 followed by +!> @details: the application of a Gaussian filter. +!> @details: for large FWHM values, a random orientation from a uniform distribution is picked +!> @details: followed by tge aookucatuib if a Gaussian filter. Additionally, the misorientation is +!> @details: limited to 2*FWHM, !-------------------------------------------------------------------------------------------------- function math_sampleGaussOri(center,FWHM) - use prec, only: & - tol_math_check implicit none real(pReal), intent(in) :: FWHM real(pReal), dimension(3), intent(in) :: center - real(pReal) :: omega + real(pReal) :: angle real(pReal), dimension(3) :: math_sampleGaussOri, axis real(pReal), dimension(2) :: rnd real(pReal), dimension(3,3) :: R - noScatter: if (FWHM < 0.5_pReal*INRAD) then + noScatter: if (FWHM < 0.1_pReal*INRAD) then math_sampleGaussOri = center else noScatter GaussConvolution: do @@ -1786,19 +1792,20 @@ function math_sampleGaussOri(center,FWHM) rnd = halton([3_pInt,6_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)] + sqrt(1.0-axis(1)**2.0_pReal)*sin(rnd(2)*2.0*PI)] ! random axis do rnd = halton([14_pInt,10_pInt]) - omega = (rnd(1)*(2.0_pReal*FWHM)**3.0_pReal)**(1.0_pReal/3.0_pReal) - if (rnd(2) < sin(omega)**2.0_pReal/omega**2.0_pReal) exit + angle = (rnd(1)*(2.0_pReal*FWHM)**3.0_pReal)**(1.0_pReal/3.0_pReal) ! maximum misorientation of 2*FWHM + if (rnd(2) < sin(angle)**2.0_pReal/angle**2.0_pReal) exit ! rejection sampling enddo - R = math_axisAngleToR(axis,omega) + R = math_axisAngleToR(axis,angle) else selectiveSampling R = math_EulerToR(math_sampleRandomOri()) endif selectiveSampling rnd = halton([8_pInt,11_pInt]) - omega = math_EulerMisorientation([0.0_pReal,0.0_pReal,0.0_pReal],math_RtoEuler(R)) - if (rnd(1) <= exp(-4.0_pReal*log(2.0_pReal)*(omega/FWHM)**2_pReal)) exit + angle = math_EulerMisorientation([0.0_pReal,0.0_pReal,0.0_pReal],math_RtoEuler(R)) + if (rnd(1) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal) .and. & ! rejection sampling (Gaussian) + angle < 2.0_pReal * FWHM) exit ! limit (in case of non-selective orientation selection enddo GaussConvolution math_sampleGaussOri = math_RtoEuler(math_mul33x33(R,math_EulerToR(center))) @@ -1807,74 +1814,52 @@ function math_sampleGaussOri(center,FWHM) end function math_sampleGaussOri - !-------------------------------------------------------------------------------------------------- !> @brief draw a sample from an Gaussian distribution around given fiber texture and Full Width ! at Half Maximum (FWHM) +!>@details: vector in cone around axis is uniformly distributed according to +! https://math.stackexchange.com/questions/56784 !------------------------------------------------------------------------------------------------- function math_sampleFiberOri(alpha,beta,FWHM) - use prec, only: & - tol_math_check implicit none real(pReal), dimension(2), intent(in) :: alpha,beta real(pReal), intent(in) :: FWHM - real(pReal), dimension(3) :: math_sampleFiberOri, fiberInC,fiberInS,axis - real(pReal), dimension(6) :: rnd - real(pReal), dimension(3,3) :: oRot,fRot,pRot - real(pReal) :: cos2Scatter, angle + real(pReal), dimension(3) :: math_sampleFiberOri, & + fInC,& !< fiber axis in crystal coordinate system + fInS,& !< fiber axis in sample coordinate system + axis + real(pReal), dimension(5) :: rnd + real(pReal), dimension(3,3) :: & + R_o, & !< rotation to aling fiber axis in crystal and sample system + R_f, & !< random rotation along fiber axis [0, 2*Pi[ + R_p !< deviation of axis alingment, bound by 2*FWHM + real(pReal) :: angle -! MD: this is a leftover from using the wrongly scaled Gaussian distribution. -! I'm relatively sure that it is not correct to scale by any constant factor here - cos2Scatter = cos(2.0_pReal*0.95_pReal * FWHM) + 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_o = math_EulerAxisAngleToR(math_crossproduct(fInC,fInS),-acos(dot_product(fInC,fInS))) - oRot = merge(math_EulerAxisAngleToR(math_crossproduct(fiberInC,fiberInS),angle), - math_I3, acos(dot_product(fiberInC,fiberInS)) > tol_math_check) - - GaussConvolution: do - rnd = halton(int([14,5,10,3,9,17],pInt)) - - ! random rotation around fiber axis - fRot = math_EulerAxisAngleToR(fiberInS,rnd(1)*2.0_pReal*PI) - - ! 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)) - else if(abs(fiberInS(2)) > tol_math_check) then - axis(3)=axis(2) - axis(2)=-(axis(1)*fiberInS(1)+axis(3)*fiberInS(3)) - else if(abs(fiberInS(1)) > tol_math_check) then - axis(3)=axis(1) - axis(1)=-(axis(2)*fiberInS(2)+axis(3)*fiberInS(3)) - end if - axis = axis/norm2(axis) - - ! rotate by XXX amount around axis perpendicular to fiber axis - ! Proabably https://math.stackexchange.com/questions/56784 is the right approach - if (FWHM > 0.0_pReal) then - angle = acos(cos2Scatter+(1.0_pReal-cos2Scatter)*rnd(4)) - if (rnd(5) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2.0_pReal)) exit - else - angle = 0.0_pReal - exit - end if - enddo GaussConvolution - if (rnd(6) <= 0.5) angle = -angle + if (FWHM > 0.0_pReal) then + GaussConvolution: do + rnd = halton(int([5,10,3,9,17],pInt)) + rnd(1:2) = [cos(FWHM*2.0_pReal),-1.0_pReal] + rnd(1:2)*[1.0_pReal - cos(FWHM*2.0_pReal),2.0_pReal] + axis = [sqrt(1.0_pReal - rnd(2)**2.0_pReal)*sin(rnd(1)),& + sqrt(1.0_pReal - rnd(2)**2.0_pReal)*cos(rnd(1)),rnd(2)] + angle = acos(dot_product([0.0_pReal,0.0_pReal,1.0_pReal],axis)) + if (rnd(3) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2.0_pReal)) exit + enddo GaussConvolution + if (rnd(4) <= 0.5) angle = -angle + R_p = math_EulerAxisAngleToR(math_crossproduct(axis,[0.0_pReal,0.0_pReal,1.0_pReal]),angle) + else + R_p = math_I3 + rnd = halton(int([5,10,3,9,17],pInt)) + endif - pRot = math_EulerAxisAngleToR(axis,angle) + R_f = math_EulerAxisAngleToR(fInS,rnd(5)*2.0_pReal*PI) -! ---# apply the three rotations #--- - math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))) + math_sampleFiberOri = math_RtoEuler(math_mul33x33(R_p,math_mul33x33(R_f,R_o))) end function math_sampleFiberOri @@ -2279,19 +2264,23 @@ end function math_invariantsSym33 !------------------------------------------------------------------------------------------------- !> @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 +!> @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. !------------------------------------------------------------------------------------------------- function halton(bases) implicit none integer(pInt), intent(in), dimension(:):: & - bases !< dimension of the sequence + bases !< bases (prime number ID) real(pReal), dimension(size(bases)) :: & halton integer(pInt), save :: & @@ -2301,8 +2290,6 @@ function halton(bases) integer(pInt), dimension(size(bases)) :: & base, & t - integer(pInt) :: & - i integer(pInt), dimension(0:1600), parameter :: & prime = int([& 1, & @@ -2496,17 +2483,6 @@ function halton(bases) t = t / base enddo - !-------------------------------------------------------------------------------------------------- - !> @brief returns any of the first 1600 prime numbers. - !> @details n = 0 is legal, returning PRIME = 1. - !> @details 0 > n > 1600 returns -1 - !> @details Reference: - !> @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 - !-------------------------------------------------------------------------------------------------- end function halton From 8fa8427c993a8d842c5e2dea6cef5418f825a454 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 8 Mar 2018 17:35:59 +0100 Subject: [PATCH 14/19] simple test, only ensures that scatter is not bigger than 2 FWHM (hard max) --- .gitlab-ci.yml | 7 +++++++ PRIVATE | 2 +- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index df8991900..f38bdb99b 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -324,6 +324,13 @@ HybridIA: - master - release +TextureComponents: + stage: spectral + script: TextureComponents/test.py + except: + - master + - release + ################################################################################################### Marc_compileIfort2014: stage: compileMarc2014 diff --git a/PRIVATE b/PRIVATE index 8de4f792a..9064217d4 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 8de4f792a46d98d98418dbf2d3e621b22c13b18a +Subproject commit 9064217d4b2753715b896b2b074eef96173cd751 From 44e5644e78795568d11d8db0b736f44d5e5a647e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 9 Mar 2018 16:46:16 +0100 Subject: [PATCH 15/19] fixed random Gaussian sampling sampling needs to be performed from unfiform misorientation, NOT uniformly distributed rotations for Fiber, compute uniform tilt of Fiber axis --- src/math.f90 | 111 ++++++++++++++++++++++++--------------------------- 1 file changed, 52 insertions(+), 59 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 80c678164..70879cf77 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1765,14 +1765,8 @@ end function math_sampleRandomOri !-------------------------------------------------------------------------------------------------- !> @brief draw a sample from an Gaussian distribution around given orientation and Full Width ! at Half Maximum (FWHM) -!> @details: for very small FWHM values the given orientation is returned. -!> @details: for intermediate FWHM values, an orientation is picked from uniformly distributed -!> @details: oreintations around the nominal orientation with maximum misorientation of 2*FWHM -!> @deatils: according to https://math.stackexchange.com/questions/13133 followed by -!> @details: the application of a Gaussian filter. -!> @details: for large FWHM values, a random orientation from a uniform distribution is picked -!> @details: followed by tge aookucatuib if a Gaussian filter. Additionally, the misorientation is -!> @details: limited to 2*FWHM, +!> @details: A uniform misorientation (limited to 2*FWHM) is sampled followed by convolution with +! a Gausian distribution !-------------------------------------------------------------------------------------------------- function math_sampleGaussOri(center,FWHM) @@ -1781,35 +1775,25 @@ function math_sampleGaussOri(center,FWHM) real(pReal), dimension(3), intent(in) :: center real(pReal) :: angle real(pReal), dimension(3) :: math_sampleGaussOri, axis - real(pReal), dimension(2) :: rnd + real(pReal), dimension(4) :: rnd real(pReal), dimension(3,3) :: R - noScatter: if (FWHM < 0.1_pReal*INRAD) then - math_sampleGaussOri = center - else noScatter + if (FWHM < 0.1_pReal*INRAD) then + R = math_I3 + else GaussConvolution: do - selectiveSampling: if (FWHM*INRAD < 90.0_pReal) then - rnd = halton([3_pInt,6_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 - do - rnd = halton([14_pInt,10_pInt]) - angle = (rnd(1)*(2.0_pReal*FWHM)**3.0_pReal)**(1.0_pReal/3.0_pReal) ! maximum misorientation of 2*FWHM - if (rnd(2) < sin(angle)**2.0_pReal/angle**2.0_pReal) exit ! rejection sampling - enddo - R = math_axisAngleToR(axis,angle) - else selectiveSampling - R = math_EulerToR(math_sampleRandomOri()) - endif selectiveSampling - rnd = halton([8_pInt,11_pInt]) - angle = math_EulerMisorientation([0.0_pReal,0.0_pReal,0.0_pReal],math_RtoEuler(R)) - if (rnd(1) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal) .and. & ! rejection sampling (Gaussian) - angle < 2.0_pReal * FWHM) exit ! limit (in case of non-selective orientation selection + 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 + endif - math_sampleGaussOri = math_RtoEuler(math_mul33x33(R,math_EulerToR(center))) - endif noScatter + math_sampleGaussOri = math_RtoEuler(math_mul33x33(R,math_EulerToR(center))) end function math_sampleGaussOri @@ -1817,8 +1801,6 @@ end function math_sampleGaussOri !-------------------------------------------------------------------------------------------------- !> @brief draw a sample from an Gaussian distribution around given fiber texture and Full Width ! at Half Maximum (FWHM) -!>@details: vector in cone around axis is uniformly distributed according to -! https://math.stackexchange.com/questions/56784 !------------------------------------------------------------------------------------------------- function math_sampleFiberOri(alpha,beta,FWHM) @@ -1828,38 +1810,49 @@ function math_sampleFiberOri(alpha,beta,FWHM) real(pReal), dimension(3) :: math_sampleFiberOri, & fInC,& !< fiber axis in crystal coordinate system fInS,& !< fiber axis in sample coordinate system - axis - real(pReal), dimension(5) :: rnd - real(pReal), dimension(3,3) :: & - R_o, & !< rotation to aling fiber axis in crystal and sample system - R_f, & !< random rotation along fiber axis [0, 2*Pi[ - R_p !< deviation of axis alingment, bound by 2*FWHM - real(pReal) :: angle + 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 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))] - R_o = math_EulerAxisAngleToR(math_crossproduct(fInC,fInS),-acos(dot_product(fInC,fInS))) + R = math_EulerAxisAngleToR(math_crossproduct(fInC,fInS),-acos(dot_product(fInC,fInS))) !< rotation to align fiber axis in crystal and sample system - if (FWHM > 0.0_pReal) then + 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=[b,i] + else + j = i + endif + enddo reducedTo2D GaussConvolution: do - rnd = halton(int([5,10,3,9,17],pInt)) - rnd(1:2) = [cos(FWHM*2.0_pReal),-1.0_pReal] + rnd(1:2)*[1.0_pReal - cos(FWHM*2.0_pReal),2.0_pReal] - axis = [sqrt(1.0_pReal - rnd(2)**2.0_pReal)*sin(rnd(1)),& - sqrt(1.0_pReal - rnd(2)**2.0_pReal)*cos(rnd(1)),rnd(2)] - angle = acos(dot_product([0.0_pReal,0.0_pReal,1.0_pReal],axis)) - if (rnd(3) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2.0_pReal)) exit - enddo GaussConvolution - if (rnd(4) <= 0.5) angle = -angle - R_p = math_EulerAxisAngleToR(math_crossproduct(axis,[0.0_pReal,0.0_pReal,1.0_pReal]),angle) - else - R_p = math_I3 - rnd = halton(int([5,10,3,9,17],pInt)) - endif - - R_f = math_EulerAxisAngleToR(fInS,rnd(5)*2.0_pReal*PI) + 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) - math_sampleFiberOri = math_RtoEuler(math_mul33x33(R_p,math_mul33x33(R_f,R_o))) + 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 + math_sampleFiberOri = math_RtoEuler(R) end function math_sampleFiberOri From 005c9f80c588e867edee4eaed1971620d52cb721 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 13 Mar 2018 03:37:28 +0100 Subject: [PATCH 16/19] include updated test --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 9064217d4..563d327d9 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 9064217d4b2753715b896b2b074eef96173cd751 +Subproject commit 563d327d92b3e22e6845080f9971e400888a8c9d From 2ceb3aeea08b4a4fa33f63faa98419eb149c956a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 13 Mar 2018 03:56:02 +0100 Subject: [PATCH 17/19] not all variables were consistently renamed --- src/math.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/math.f90 b/src/math.f90 index 70879cf77..742ab2f2c 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1831,7 +1831,7 @@ function math_sampleFiberOri(alpha,beta,FWHM) reducedTo2D: do i=1_pInt,3_pInt if (i /= minloc(abs(fInS),1)) then a=[a,fInS(i)] - idx=[b,i] + idx=[idx,i] else j = i endif From 8509b13f168c1fc3301da4abe3f70df1c7e37690 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 13 Mar 2018 09:21:41 +0100 Subject: [PATCH 18/19] avoid numerical noise --- src/math.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 0c6f427ab..8eac3e304 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1779,7 +1779,7 @@ function math_sampleGaussOri(center,FWHM) real(pReal), dimension(3,3) :: R if (FWHM < 0.1_pReal*INRAD) then - R = math_I3 + math_sampleGaussOri = center else GaussConvolution: do rnd = halton([8_pInt,3_pInt,6_pInt,11_pInt]) @@ -1791,9 +1791,9 @@ function math_sampleGaussOri(center,FWHM) 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 - math_sampleGaussOri = math_RtoEuler(math_mul33x33(R,math_EulerToR(center))) end function math_sampleGaussOri From 0473d06c6fc40ffb1a54fb02580c6dcc93969119 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 10 Apr 2018 10:27:03 +0200 Subject: [PATCH 19/19] using updated test --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index af8516892..7c69abfc5 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit af851689285b8c1a633495219abd9dbbd5a11c69 +Subproject commit 7c69abfc5bf54c083b9096511abde7d74b806b7f