c     f77 -o mckpts mckpts.f; mckpts > mckpts.2000
      PROGRAM MCKPTS
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NUM_KPTS = 2000)
      DIMENSION b(3,3)
      
      b(1,1) = 1.d0
      b(2,1) = 1./sqrt(3.d0)
      b(3,1) = 0.d0
      
      b(1,2) = 0.d0
      b(2,2) = 2/sqrt(3.d0)
      b(3,2) = 0.d0
      
      b(1,3) = 0.d0
      b(2,3) = 0.d0
      b(3,3) = sqrt(3.d0/8.d0)
      
      iseed = 10
      write(6,*) NUM_KPTS 
      do k=1,NUM_KPTS 
      c1 = ran1(iseed)
      c2 = ran1(iseed)
      c3 = ran1(iseed)
      fkx = 2*(c1*b(1,1) + c2*b(1,2) + c3*b(1,3))
      fky = 2*(c1*b(2,1) + c2*b(2,2) + c3*b(2,3))
      fkz = 2*(c1*b(3,1) + c2*b(3,2) + c3*b(3,3))
      write(6,10) fkx, fky, fkz, 1.d0/NUM_KPTS
10    format(4f18.14)
      enddo
      
      stop
      end
      
      FUNCTION RAN1(IDUM)
C     ---------------------------------------------------------
      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
      SAVE IV,IY
      DATA IV /NTAB*0/, IY /0/
      IF (IDUM.LE.0.OR.IY.EQ.0) THEN
        IDUM=MAX(-IDUM,1)
        DO 11 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
11      CONTINUE
        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
C     ---------------------------------------------------------