SUBROUTINE gasdev_s(harvest)
	USE nrtype
	USE nr, ONLY : ran1
	IMPLICIT NONE
	REAL(SP), INTENT(OUT) :: harvest
	REAL(SP) :: rsq,v1,v2
	REAL(SP), SAVE :: g
	LOGICAL, SAVE :: gaus_stored=.false.
	if (gaus_stored) then
		harvest=g
		gaus_stored=.false.
	else
		do
			call ran1(v1)
			call ran1(v2)
			v1=2.0_sp*v1-1.0_sp
			v2=2.0_sp*v2-1.0_sp
			rsq=v1**2+v2**2
			if (rsq > 0.0 .and. rsq < 1.0) exit
		end do
		rsq=sqrt(-2.0_sp*log(rsq)/rsq)
		harvest=v1*rsq
		g=v2*rsq
		gaus_stored=.true.
	end if
	END SUBROUTINE gasdev_s

	SUBROUTINE gasdev_v(harvest)
	USE nrtype; USE nrutil, ONLY : array_copy
	USE nr, ONLY : ran1
	IMPLICIT NONE
	REAL(SP), DIMENSION(:), INTENT(OUT) :: harvest
	REAL(SP), DIMENSION(size(harvest)) :: rsq,v1,v2
	REAL(SP), ALLOCATABLE, DIMENSION(:), SAVE :: g
	INTEGER(I4B) :: n,ng,nn,m
	INTEGER(I4B), SAVE :: last_allocated=0
	LOGICAL, SAVE :: gaus_stored=.false.
	LOGICAL, DIMENSION(size(harvest)) :: mask
	n=size(harvest)
	if (n /= last_allocated) then
		if (last_allocated /= 0) deallocate(g)
		allocate(g(n))
		last_allocated=n
		gaus_stored=.false.
	end if
	if (gaus_stored) then
		harvest=g
		gaus_stored=.false.
	else
		ng=1
		do
			if (ng > n) exit
			call ran1(v1(ng:n))
			call ran1(v2(ng:n))
			v1(ng:n)=2.0_sp*v1(ng:n)-1.0_sp
			v2(ng:n)=2.0_sp*v2(ng:n)-1.0_sp
			rsq(ng:n)=v1(ng:n)**2+v2(ng:n)**2
			mask(ng:n)=(rsq(ng:n)>0.0 .and. rsq(ng:n)<1.0)
			call array_copy(pack(v1(ng:n),mask(ng:n)),v1(ng:),nn,m)
			v2(ng:ng+nn-1)=pack(v2(ng:n),mask(ng:n))
			rsq(ng:ng+nn-1)=pack(rsq(ng:n),mask(ng:n))
			ng=ng+nn
		end do
		rsq=sqrt(-2.0_sp*log(rsq)/rsq)
		harvest=v1*rsq
		g=v2*rsq
		gaus_stored=.true.
	end if
	END SUBROUTINE gasdev_v