betarnd Function

function betarnd(a,b) result(ans)
  real(8) :: a,b,ans,u,v
  if ((a <= 0d0) .or. (b <= 0d0)) then
     write(*,*) "BETARND: Beta parameters must be positive"
  end if
  !    ## There are more efficient methods for generating beta samples.
  !    ## However such methods are a little more efficient and much more complicated.
  !    ## For an explanation of why the following method works, see
  !    ## http://www.johndcook.com/distribution_chart.html#gamma_beta
  u = gammarnd(a,1d0)
  v = gammarnd(b,1d0)
  ans = u / (u + v)
end function betarnd