223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

using proper Gauss sampling also for the fiber components 

function most probably still contains a bu
This commit is contained in:
Martin Diehl 2018-02-22 00:42:06 +01:00
parent 5c908e44ec
commit c6c66bb653
1 changed files with 7 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
src/math.f90
View File

@ -1828,8 +1828,8 @@ function math_sampleFiberOri(alpha,beta,FWHM)
1_pInt,2_pInt],[2,3]) 1_pInt,2_pInt],[2,3])
integer(pInt) :: i integer(pInt) :: i
! Helming uses different distribution with Bessel functions ! MD: this is a leftover from using the wrongly scaled Gaussian distribution.
! therefore the gauss scatter width has to be scaled differently ! I'm relatively sure that it is not correct to scale by any constant factor here
scatter = 0.95_pReal * FWHM scatter = 0.95_pReal * FWHM
cos2Scatter = cos(2.0_pReal*scatter) cos2Scatter = cos(2.0_pReal*scatter)
@ -1853,9 +1853,9 @@ function math_sampleFiberOri(alpha,beta,FWHM)
end if end if
! ---# rotation matrix about fiber axis (random angle) #--- ! ---# rotation matrix about fiber axis (random angle) #---
do GaussConvolution: do
rnd = halton(6_pInt) 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 #--- ! ---# rotation about random axis perpend to fiber #---
! random axis pependicular to fiber axis ! random axis pependicular to fiber axis
@ -1872,13 +1872,13 @@ function math_sampleFiberOri(alpha,beta,FWHM)
! scattered rotation angle ! scattered rotation angle
if (FWHM > 0.0_pReal) then if (FWHM > 0.0_pReal) then
angle = acos(cos2Scatter+(1.0_pReal-cos2Scatter)*rnd(4)) angle = acos(cos2Scatter+(1.0_pReal-cos2Scatter)*rnd(4)) ! MD: This is probably not ok
if (rnd(5) <= exp(-1.0_pReal*(angle/scatter)**2.0_pReal)) exit if (rnd(5) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2.0_pReal)) exit
else else
angle = 0.0_pReal angle = 0.0_pReal
exit exit
end if end if
enddo enddo GaussConvolution
if (rnd(6) <= 0.5) angle = -angle if (rnd(6) <= 0.5) angle = -angle
pRot = math_EulerAxisAngleToR(axis,angle) pRot = math_EulerAxisAngleToR(axis,angle)