c pickcard1.for			J. C. Lahr  July 2006
c This program will repeat this sequence one million times:
c	Shuffel the deck
c	Deal cards until all four suits are represented and save that number of cards
c Write the final set of numbers to the file cards.txt
      dimension inonrand(52), irand(52), ihist(52)
      data inonrand /13*1, 13*2, 13*3, 13*4/
      logical used(52), suits(4)



      write(6, '(a)') ' Begin pickcard1'
      open(UNIT=9,FILE='hist.txt',ACTION='write',STATUS='new',
     *  IOSTAT=ierror)
      if(ierror .ne. 0) then
	open(UNIT=9,FILE='cards.txt',ACTION='write',STATUS='old')
	rewind(9)
      endif
      iseed = 100
      dum = ran3(-iseed)
      print *, 'Finished seeding ran3'
      iseed = 0
      n = 0

      do i = 1, 52
        ihist(i) = 0
      enddo

15    do i = 1, 52
        used(i) = .false.
      enddo
      do i = 1, 4
	suits(i) = .false.
      enddo


c shuffle the deck irand(52)
c produce a random real number between 0 and 51.999
      do i = 1, 52
20      dum = 51.999*ran3(iseed)
c add one, so range becomes 1 to 52.999
	dum = dum + 1.0
c truncate to produce and integer between 1 and 52.
	icard = dum
        if(used(icard)) goto 20
        irand(i) = inonrand(icard)
	used(icard) = .true.
      enddo
c deal cards until all suits are represented
      do i = 1, 52
        suits(irand(i)) = .true.
        if(suits(1) .and. suits(2) .and. suits(3) .and. suits(4))
     *    goto 30
      enddo

30    n = n + 1
      ihist(i) = ihist(i) + 1
c     write(9, '(i2)') i
      if(n .lt. 1000000) goto 15

      do i = 1, 52
	write(9, '(i6)') ihist(i)
      enddo
      stop
      end


c ran3.for    []
      function ran3(idum)
      implicit real(m)
      save inextp, inext, ma, mj, mz, fac
c portable random number generator from press and others, "numerical
c recipes, the art of scientific computing"
      parameter (mbig=4000000.,mseed=1618033.,mz=0.,fac=2.5e-7)
c         parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
      dimension ma(55)
      data iff /0/
      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 11 i=1,54
          ii=mod(21*i,55)
          ma(ii)=mk
          mk=mj-mk
          if(mk.lt.mz)mk=mk+mbig
          mj=ma(ii)
11      continue
        do 13 k=1,4
          do 12 i=1,55
            ma(i)=ma(i)-ma(1+mod(i+30,55))
            if(ma(i).lt.mz)ma(i)=ma(i)+mbig
12        continue
13      continue
        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=mj*fac
      return
      end
c end ran3