c FILE : random.f
c ----------------------------------------------------------------------
c ------ Uniform Distribution Random Number Generator ------------------
c ----------------------------------------------------------------------
c  From Numerical Recipies.

	FUNCTION ran3(idum)
	INTEGER*2 i,ii,k,inext,inextp,iff
	INTEGER*4 idum,mbig,mseed,mz,mj,ma(55),mk
	REAL*4 fac, ran3
	parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
	save inext,inextp,ma,iff
	data iff /0/
C                                                    
	if(idum.lt.0.or.iff.eq.0)then
	      iff=1
	      mj=mseed-iabs(idum)
	      mj=mod(mj,mbig)
	      ma(55)=mj
	      mk=1
	      do i=1,54
		   ii=mod(21*i,55)
		   ma(ii)=mk
		   mk=mj-mk
		   if(mk.lt.mz) mk=mk+mbig
		   mj=ma(ii)
	      end do
	      do k=1,4
		   do i=1,55
			ma(i)=ma(i)-ma(1+mod(i+30,55))
			if(ma(i).lt.mz)ma(i)=ma(i)+mbig
		   end do
	      end do
	      inext=0
	      inextp=31
	      idum=1
	endif
	inext=inext+1
	if(inext.eq.56)inext=1
	inextp=inextp+1
	if(inextp.eq.56)inextp=1
	mj=ma(inext)-ma(inextp)
	if(mj.lt.mz)mj=mj+mbig
	ma(inext)=mj
	ran3=real(mj)*fac
	return
	end


	FUNCTION ran4(idum)
	INTEGER*2 i,ii,k,inext,inextp,iff
	INTEGER*4 idum,mbig,mseed,mz,mj,ma(55),mk
	REAL*4 fac, ran4
	parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
	save inext,inextp,ma,iff
	data iff /0/
C                                                    
	if(idum.lt.0.or.iff.eq.0)then
	      iff=1
	      mj=mseed-iabs(idum)
	      mj=mod(mj,mbig)
	      ma(55)=mj
	      mk=1
	      do i=1,54
		   ii=mod(21*i,55)
		   ma(ii)=mk
		   mk=mj-mk
		   if(mk.lt.mz) mk=mk+mbig
		   mj=ma(ii)
	      end do
	      do k=1,4
		   do i=1,55
			ma(i)=ma(i)-ma(1+mod(i+30,55))
			if(ma(i).lt.mz)ma(i)=ma(i)+mbig
		   end do
	      end do
	      inext=0
	      inextp=31
	      idum=1
	endif
	inext=inext+1
	if(inext.eq.56)inext=1
	inextp=inextp+1
	if(inextp.eq.56)inextp=1
	mj=ma(inext)-ma(inextp)
	if(mj.lt.mz)mj=mj+mbig
	ma(inext)=mj
	ran4=real(mj)*fac
	return
	end


c ---------------------------------------------------------------------
c -------- Normally Distributed Random number generator ---------------
c ---------------------------------------------------------------------
c From numerical recipies
c Mean=0, Stdev=1

	FUNCTION normdis(idum)
	INTEGER iset, idum
	REAL    v1, v2, r, fac, normdis, ran4, gset	
	SAVE gset, iset
	DATA iset/0/
c	idum=0
	IF (iset.eq.0) THEN
1		v1=2.*RAN4(idum)-1.
		v2=2.*RAN4(idum)-1.
		r=v1**2+v2**2
		IF (r.ge.1.) GOTO 1
		fac=SQRT(-2.*LOG(r)/r)
		gset=v1*fac
		normdis=v2*fac
		iset=1
	ELSE
		normdis=gset
		iset=0
	END IF
	RETURN
	END

c ---------------------------------------------------------------

	FUNCTION rbit()
	INTEGER	rbit, isd
	REAL r, ran3
	r=ran3(isd)
	IF (r.le.0.5) THEN
		rbit=1
	ELSE
		rbit=-1
	ENDIF
	RETURN
	END