1C     program gasdev1
2C     double precision gks_gasdev,gasdev,seed
3C      seed = -1234
4C      write(*,*)  'before gasdev'
5C      gasdev=gks_gasdev(seed)
6C      write(*,*)  'Afer gasdev'
7C     end
8C    rundom number generators from gaussian distribution,
9C    from Numerical Recipies
10C    it is initialized with -IDUM
11
12      double precision FUNCTION GKS_GASDEV(IDUM)
13      save
14c     implicit undefined (a-z)
15      integer idum,iset
16      double precision v1,v2,r,fac,gset,gks_ran1
17      DATA ISET /0/
18c      write(*,*) 'inside GKS_GASDEV'
19      IF(ISET.EQ.0) THEN
20   10   V1=2.0E+00*GKS_RAN1(IDUM)-1.0E+00
21        V2=2.0E+00*GKS_RAN1(IDUM)-1.0E+00
22        R=V1**2+V2**2
23        IF(R.GE.1.0E+00) GO TO 10
24        FAC=SQRT(-2.0E+00*LOG(R)/R)
25        GSET      =V1*FAC
26        GKS_GASDEV=V2*FAC
27        ISET=1
28      ELSE
29        GKS_GASDEV=GSET
30        ISET=0
31      ENDIF
32c      write(*,*) idum,iset,gks_gasdev
33C9999 format(' ---gasdev--- = ',2i15,f10.6)
34      RETURN
35      END
36ccccccccc/ccccccccc/ccccccccc/ccccccccc/ccccccccc/ccccccccc/ccccccccc/ccccccccc/
37      double precision function gks_ran1(idum)
38c
39c     N.R. FUNCTION gks_ran2(idum)
40c
41      INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
42      double precision AM,EPS,RNMX
43      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
44     1   IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
45     2   NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
46      INTEGER idum2,j,k,iv(NTAB),iy
47      SAVE iv,iy,idum2
48      DATA idum2/123456789/, iv/NTAB*0/, iy/0/
49      if (idum.le.0) then
50        idum=max(-idum,1)
51        idum2=idum
52        do 11 j=NTAB+8,1,-1
53          k=idum/IQ1
54          idum=IA1*(idum-k*IQ1)-k*IR1
55          if (idum.lt.0) idum=idum+IM1
56          if (j.le.NTAB) iv(j)=idum
5711      continue
58        iy=iv(1)
59      endif
60      k=idum/IQ1
61      idum=IA1*(idum-k*IQ1)-k*IR1
62      if (idum.lt.0) idum=idum+IM1
63      k=idum2/IQ2
64      idum2=IA2*(idum2-k*IQ2)-k*IR2
65      if (idum2.lt.0) idum2=idum2+IM2
66      j=1+iy/NDIV
67      iy=iv(j)-idum2
68      iv(j)=idum
69      if(iy.lt.1)iy=iy+IMM1
70      gks_ran1=min(AM*iy,RNMX)
71      return
72      END
73ccccccc
74
75c $Id$
76