rang Subroutine

        subroutine rang(x, y, avg, sigma)
        !! sample a 2D Guassian distribution

            !> mean of the gaussian to sample from
            real(kind=wp), intent(IN)  :: avg
            !> \(\sigma\) of the guassian to sample from.
            real(kind=wp), intent(IN)  :: sigma
            !> first value to return
            real(kind=wp), intent(OUT) :: x
            !> 2nd value to return
            real(kind=wp), intent(OUT) :: y            
            
            real(kind=wp) :: s, tmp

            s = 1._wp

            do while(s >= 1._wp)
                x = ranu(-1._wp, 1._wp)
                y = ranu(-1._wp, 1._wp)
                s = y**2 + x**2
            end do

            tmp = x*sqrt(-2._wp*log(s)/s)
            x = avg + sigma*tmp

            tmp = y*sqrt(-2._wp*log(s)/s)
            y = avg + sigma*tmp

        end subroutine rang