Commit 4062a958 authored by ronesy's avatar ronesy
Browse files

Change to normal distribution for random perturbations

parent ddf4d9f1
......@@ -31,7 +31,7 @@ module mod_perturb
implicit none
private
public :: perturb_flux, perturb_obs, randnum
public :: perturb_flux, perturb_obs, randnum, random_normal
contains
......@@ -116,6 +116,11 @@ module mod_perturb
write(100,fmt='(F14.4)') perturb(i)
end do
close(100)
open(100,file=trim(files%path_output)//'rvec.txt',status='replace',action='write')
do i = 1, nvar
write(100,fmt='(F10.8)') rvec(i)
end do
close(100)
end subroutine perturb_flux
......@@ -196,7 +201,8 @@ module mod_perturb
end do
call random_seed(put=seedvec(1:ssize))
do i = 1, num
call random_number(ri)
! call random_number(ri)
ri = random_normal()
rvec(i) = ri
end do
deallocate(seedvec)
......@@ -210,6 +216,61 @@ module mod_perturb
end subroutine randnum
! --------------------------------------------------
! random_normal
! --------------------------------------------------
!> random_normal
!!
!! Adapted from the following Fortran 77 code
!! ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
!! THIS WORK PUBLISHED IN TRANSACTIONS ON
!! MATHEMATICAL SOFTWARE, VOL. 18, NO. 4,
!! DECEMBER, 1992, PP. 434-435.
!!
!! Purpose:
!! The function random_normal() returns a normally
!! distributed pseudo-random number with zero mean
!! and unit variance.
!!
!! The algorithm uses the ratio of uniforms method
!! of A.J. Kinderman and J.F. Monahan augmented
!! with quadratic bounding curves.
!!
! --------------------------------------------------
real function random_normal()
implicit none
real :: s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472, &
r1 = 0.27597, r2 = 0.27846, u, v, x, y, q, half = 0.5
! Generate P = (u,v) uniform in rectangle enclosing acceptance region
do
call random_number(u)
call random_number(v)
v = 1.7156 * (v - half)
! Evaluate the quadratic form
x = u - s
y = abs(v) - t
q = x**2 + y*(a*y - b*x)
! Accept P if inside inner ellipse
if (q < r1) exit
! Reject P if outside outer ellipse
if (q > r2) cycle
! Reject P if outside acceptance region
if (v**2 < -4.0*log(u)*u**2) exit
end do
! Return ratio of P's coordinates as the normal deviate
random_normal = v/u
return
end function random_normal
! --------------------------------------------------
end module mod_perturb
......
......@@ -45,7 +45,7 @@ module mod_perturb
implicit none
private
public :: perturb_state, perturb_obs, randnum
public :: perturb_state, perturb_obs, randnum, random_normal
contains
......@@ -158,7 +158,8 @@ module mod_perturb
end do
call random_seed(put=seedvec(1:ssize))
do i = 1, num
call random_number(ri)
! call random_number(ri)
ri = random_normal()
rvec(i) = ri
end do
deallocate(seedvec)
......@@ -172,6 +173,61 @@ module mod_perturb
end subroutine randnum
! --------------------------------------------------
! random_normal
! --------------------------------------------------
!> random_normal
!!
!! Adapted from the following Fortran 77 code
!! ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
!! THIS WORK PUBLISHED IN TRANSACTIONS ON
!! MATHEMATICAL SOFTWARE, VOL. 18, NO. 4,
!! DECEMBER, 1992, PP. 434-435.
!!
!! Purpose:
!! The function random_normal() returns a normally
!! distributed pseudo-random number with zero mean
!! and unit variance.
!!
!! The algorithm uses the ratio of uniforms method
!! of A.J. Kinderman and J.F. Monahan augmented
!! with quadratic bounding curves.
!!
! --------------------------------------------------
real function random_normal()
implicit none
real :: s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472, &
r1 = 0.27597, r2 = 0.27846, u, v, x, y, q, half = 0.5
! Generate P = (u,v) uniform in rectangle enclosing acceptance region
do
call random_number(u)
call random_number(v)
v = 1.7156 * (v - half)
! Evaluate the quadratic form
x = u - s
y = abs(v) - t
q = x**2 + y*(a*y - b*x)
! Accept P if inside inner ellipse
if (q < r1) exit
! Reject P if outside outer ellipse
if (q > r2) cycle
! Reject P if outside acceptance region
if (v**2 < -4.0*log(u)*u**2) exit
end do
! Return ratio of P's coordinates as the normal deviate
random_normal = v/u
return
end function random_normal
! --------------------------------------------------
end module mod_perturb
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment