rnorm Function

public function rnorm() result(fn_val)

Generate a random normal deviate using the polar method.

reference: marsaglia,g. & bray,t.a. ‘a convenient method for generating normal variables’, siam rev., vol.6, 260-264, 1964.

Arguments

None

Return Value real(kind=R8)


Called by

proc~~rnorm~~CalledByGraph proc~rnorm rnorm proc~rnorm_vec rnorm_vec proc~rnorm_vec->proc~rnorm program~test_stat test_stat program~test_stat->proc~rnorm_vec

Source Code

   function rnorm() result(fn_val)
   !================================================================================================
   !! Generate a random normal deviate using the polar method.
   !!
   !! *reference*: marsaglia,g. & bray,t.a. 'a convenient method for generating
   !!              normal variables', siam rev., vol.6, 260-264, 1964.
   !!
   !------------------------------------------------------------------------------------------------
   implicit none
   real(kind=R8) :: fn_val

      real(kind=R8)            :: u, sum
      real(kind=R8),    save   :: v, sln
      logical(kind=I4), save   :: second = .false.
      real(kind=R8), parameter :: one = 1.0_R8, vsmall = tiny( one )

      if (second) then
      ! if second, use the second random number generated on last call

         second = .false.
         fn_val = v*sln

      else
      ! first call; generate a pair of random normals

         second = .true.
         do
            call random_number( u )
            call random_number( v )
            u = scale( u, 1 ) - one
            v = scale( v, 1 ) - one
            sum = u*u + v*v + vsmall         ! vsmall added to prevent log(zero) / zero
            if( sum < one ) exit
         enddo

         sln = sqrt(-scale(log(sum),1)/sum)
         fn_val = u*sln

      endif

   return
   endfunction rnorm