use the Box-Muller transform instead of random sampling
still needs testing.
This commit is contained in:
parent
430b947c7e
commit
f40d731fe1
45
src/math.f90
45
src/math.f90
|
@ -961,39 +961,36 @@ pure function math_3333toVoigt66(m3333)
|
||||||
end function math_3333toVoigt66
|
end function math_3333toVoigt66
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
!> @brief draw a random sample from Gauss variable
|
!> @brief Draw a sample from a normal distribution.
|
||||||
|
!> @details https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
|
||||||
|
!> https://masuday.github.io/fortran_tutorial/random.html
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
real(pReal) function math_sampleGaussVar(mu, sigma, width)
|
impure elemental subroutine math_normal(x,mu,sigma)
|
||||||
|
|
||||||
real(pReal), intent(in) :: mu, & !< mean
|
real(pReal), intent(out) :: x
|
||||||
sigma !< standard deviation
|
real(pReal), intent(in), optional :: mu, sigma
|
||||||
real(pReal), intent(in), optional :: width !< cut off as multiples of standard deviation
|
|
||||||
|
|
||||||
real(pReal), dimension(2) :: rnd ! random numbers
|
real(pReal) :: sigma_, mu_
|
||||||
real(pReal) :: scatter, & ! normalized scatter around mean
|
real(pReal), dimension(2) :: rnd
|
||||||
width_
|
|
||||||
|
|
||||||
if (abs(sigma) < tol_math_check) then
|
|
||||||
math_sampleGaussVar = mu
|
if (present(mu)) then
|
||||||
|
mu_ = mu
|
||||||
else
|
else
|
||||||
if (present(width)) then
|
mu_ = 0.0_pReal
|
||||||
width_ = width
|
end if
|
||||||
else
|
|
||||||
width_ = 3.0_pReal ! use +-3*sigma as default scatter
|
if (present(sigma)) then
|
||||||
endif
|
sigma_ = sigma
|
||||||
|
else
|
||||||
|
sigma_ = 1.0_pReal
|
||||||
|
end if
|
||||||
|
|
||||||
do
|
|
||||||
call random_number(rnd)
|
call random_number(rnd)
|
||||||
scatter = width_ * (2.0_pReal * rnd(1) - 1.0_pReal)
|
x = mu_ + sigma_ * sqrt(-2.0_pReal*log(1.0_pReal-rnd(1)))*cos(2.0_pReal*PI*(1.0_pReal - rnd(2)))
|
||||||
if (rnd(2) <= exp(-0.5_pReal * scatter**2)) exit ! test if scattered value is drawn
|
|
||||||
enddo
|
|
||||||
|
|
||||||
math_sampleGaussVar = scatter * sigma
|
end subroutine math_normal
|
||||||
endif
|
|
||||||
|
|
||||||
end function math_sampleGaussVar
|
|
||||||
|
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
|
|
@ -1592,21 +1592,15 @@ subroutine stateInit(ini,phase,Nentries)
|
||||||
stt%rhoSglMobile(s,e) = densityBinning
|
stt%rhoSglMobile(s,e) = densityBinning
|
||||||
end do
|
end do
|
||||||
else ! homogeneous distribution with noise
|
else ! homogeneous distribution with noise
|
||||||
do e = 1, Nentries
|
|
||||||
do f = 1,size(ini%N_sl,1)
|
do f = 1,size(ini%N_sl,1)
|
||||||
from = 1 + sum(ini%N_sl(1:f-1))
|
from = 1 + sum(ini%N_sl(1:f-1))
|
||||||
upto = sum(ini%N_sl(1:f))
|
upto = sum(ini%N_sl(1:f))
|
||||||
do s = from,upto
|
call math_normal(stt%rho_sgl_mob_edg_pos(from:upto,:),ini%rho_u_ed_pos_0(f),ini%sigma_rho_u)
|
||||||
noise = [math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u), &
|
call math_normal(stt%rho_sgl_mob_edg_neg(from:upto,:),ini%rho_u_ed_neg_0(f),ini%sigma_rho_u)
|
||||||
math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u)]
|
call math_normal(stt%rho_sgl_mob_scr_pos(from:upto,:),ini%rho_u_sc_pos_0(f),ini%sigma_rho_u)
|
||||||
stt%rho_sgl_mob_edg_pos(s,e) = ini%rho_u_ed_pos_0(f) + noise(1)
|
call math_normal(stt%rho_sgl_mob_scr_neg(from:upto,:),ini%rho_u_sc_neg_0(f),ini%sigma_rho_u)
|
||||||
stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1)
|
stt%rho_dip_edg(from:upto,:) = ini%rho_d_ed_0(f)
|
||||||
stt%rho_sgl_mob_scr_pos(s,e) = ini%rho_u_sc_pos_0(f) + noise(2)
|
stt%rho_dip_scr(from:upto,:) = ini%rho_d_sc_0(f)
|
||||||
stt%rho_sgl_mob_scr_neg(s,e) = ini%rho_u_sc_neg_0(f) + noise(2)
|
|
||||||
end do
|
|
||||||
stt%rho_dip_edg(from:upto,e) = ini%rho_d_ed_0(f)
|
|
||||||
stt%rho_dip_scr(from:upto,e) = ini%rho_d_sc_0(f)
|
|
||||||
end do
|
|
||||||
end do
|
end do
|
||||||
end if
|
end if
|
||||||
|
|
||||||
|
|
Loading…
Reference in New Issue