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

polishing

This commit is contained in:
Martin Diehl 2020-01-10 02:11:19 +01:00
parent 8f43f05437
commit 87c7a5d5a3
1 changed files with 17 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
src/math.f90
View File

@ -398,19 +398,19 @@ pure function math_exp33(A,n)
real(pReal), dimension(3,3), intent(in) :: A
real(pReal), dimension(3,3) :: B, math_exp33
real(pReal) :: invFac
integer :: order
B = math_I3 ! init
invFac = 1.0_pReal ! 0!
math_exp33 = B ! A^0 = eye2
integer :: n_
if (present(n)) then
order = n
n_ = n
else
order = 5
n_ = 5
endif
do i = 1, order
invFac = 1.0_pReal ! 0!
B = math_I3
math_exp33 = math_I3 ! A^0 = I
do i = 1, n_
invFac = invFac/real(i,pReal) ! invfac = 1/(i!)
B = matmul(B,A)
math_exp33 = math_exp33 + invFac*B ! exp = SUM (A^i)/(i!)
@ -882,16 +882,20 @@ real(pReal) function math_sampleGaussVar(meanvalue, stddev, width)
real(pReal), intent(in), optional :: width ! width of considered values as multiples of standard deviation
real(pReal), dimension(2) :: rnd ! random numbers
real(pReal) :: scatter, & ! normalized scatter around meanvalue
myWidth
width_
if (abs(stddev) < tol_math_check) then
math_sampleGaussVar = meanvalue
else
myWidth = merge(width,3.0_pReal,present(width)) ! use +-3*sigma as default value for scatter if not given
if (present(width)) then
width_ = width
else
width_ = 3.0_pReal ! use +-3*sigma as default scatter
endif
do
call random_number(rnd)
scatter = myWidth * (2.0_pReal * rnd(1) - 1.0_pReal)
scatter = width_ * (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