      FUNCTION gasdev(idum)
      INTEGER idum
      double precision gasdev

C     USES ran1

      INTEGER iset
      double precision fac,gset,rsq,v1,v2,ran1
      SAVE iset,gset
      DATA iset/0/
      if (idum.lt.0)  iset=0
      if (iset.eq.0) then
1        v1=2.*ran1(idum)-1
         v2=2.*ran1(idum)-1
c         write(*,*) 'v1=',v1
c         write(*,*) 'v2=',v2
         rsq=v1**2+v2**2
         if(rsq.ge.1..or.rsq.eq.0.) goto 1
         fac=sqrt(-2.*log(rsq)/rsq)
         gset=v1*fac
         gasdev=v2*fac
         iset=1
      else
         gasdev=gset
         iset=0
      endif
      return
      END

      FUNCTION ran1(idum)
      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      double precision ran1,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     &       NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)

      INTEGER j,k,iv(NTAB),iy/0/
      SAVE iv,iy
      DATA iv /NTAB*0/, iy/0/
      if (idum.le.0.or.iy.eq.0) then
          idum=max(-idum,1)
          do j=NTAB+8,1,-1
             k=idum/IQ
             idum=IA*(idum-k*IQ)-IR*k
             if (idum.lt.0) idum=idum+IM
             if (j.le.NTAB) iv(j)=idum
          enddo
          iy=iv(1)
       endif
       k=idum/IQ
       idum=IA*(idum-k*IQ)-IR*k
       if(idum.lt.0) idum=idum+IM
       j=1+iy/NDIV
       iy=iv(j)
       iv(j)=idum
       ran1=min(AM*iy,RNMX)
       return
       END
