var/da/da_tools/da_gaus_noise.inc

   1 real function da_gaus_noise(KDUM)
   2
   3    !---------------------------------------------------------------------
   4    ! Purpose: return a normally(gaussian) distributed random variable ctions
   5    !     with 0 mean and 1 standard deviation
   6    !
   7    !   method :
   8    !   ------
   9    !
  10    !   input :
  11    !   -----
  12    !      argument
  13    !         kdum             : seed for random number generator
  14    !                           (set to a large negative number on entry)
  15    !      common              : none
  16    !
  17    !   output :
  18    !   ------
  19    !      argument
  20    !         gauss_noise      : gaussian random betwen
  21    !      common              : none
  22    !
  23    !   workspace :              none
  24    !   ---------
  25    !      local :
  26    !
  27    !   external :
  28    !   --------
  29    !      unifva
  30    !
  31    !   reference :
  32    !   ---------
  33    !      numerical recipes in fortran. the art of scientific computing.
  34    !      second edition.  cambridge university press.
  35    !      press et. al., 1986.
  36    !
  37    !   modifications:
  38    !   --------------
  39    !       original  : 95-01(f. vandenberghe)
  40    !       addition  : 96-06(a. weaver)
  41    !---------------------------------------------------------------------
  42
  43    implicit none
  44
  45    integer, intent(inout) :: kdum
  46
  47    integer             :: niset
  48    real                :: gset
  49
  50    save niset,gset
  51    data niset/0/
  52
  53    real zfac,zrsq,zv1,zv2
  54
  55    if (trace_use) call da_trace_entry("da_gaus_noise")
  56
  57    if (niset.eq.0) then
  58
  59 1000  continue
  60       zv1   = 2.0*da_unifva(kdum) - 1.0
  61       zv2   = 2.0*da_unifva(kdum) - 1.0
  62       zrsq  = zv1**2 + zv2**2
  63
  64       if ((zrsq>=1.0).or.(zrsq==0.0)) goto 1000
  65
  66       zfac  = sqrt(-2.0*log(zrsq)/zrsq)
  67       gset  = zv1*zfac
  68       da_gaus_noise = zv2*zfac
  69       niset = 1
  70    else
  71       da_gaus_noise = gset
  72       niset = 0
  73    end if
  74
  75    if (trace_use) call da_trace_exit("da_gaus_noise")
  76
  77 end function da_gaus_noise
  78
  79