!***********************************************************************
!************************ High Energy Accelerator Research Organization*
!** u c r a n d m _ t e s t **                                         *
!*****************************     EGS5.0 USER CODE - 27 Jul 2015/1315 *
!* Check moller sampling.                                              *
!***********************************************************************
!                                                                      *
!  PROGRAMMERS:  H. Hirayama and Y. Namito                             *
!                Radiation Science Center                              *
!                Applied Science Laboratory                            *
!                KEK, High Energy Accelerator Research Organization    *
!                1-1, Oho, Tsukuba, Ibaraki, 305-0801                  *
!                Japan                                                 *
!                                                                      *
!                E-mail:      hideo.hirayama@kek.jp                    *
!                Telephone:   +81-29-864-5489                          *
!                Fax:         +81-29-864-1993                          *
!                Based on ucrtz_homog by Nelson and James.             *
!                                                                      *
!***********************************************************************
!***********************************************************************
! This user code to test random number generator accuracy.             *
! Use Ranlux random number generator.                                  *
!***********************************************************************
!23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
!-----------------------------------------------------------------------
!------------------------------- main code -----------------------------
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!     Step 1. Initialization
!-----------------------------------------------------------------------

      implicit none

!     ------------
!     EGS5 COMMONs
!     ------------
!     include 'include/egs5_h.f'                ! Main EGS "header" file

      include 'include/randomm.f'

!     ----------------------
!     Auxiliary-code COMMONs
!     ----------------------
      include 'auxcommons/lines.f'


      real*8 srandoms(100),srandom2s(100)

      real*8 rnnow,rnnow1,rupp,rlow,srandomerr,around       ! Local variables

      real tarray(2)

      integer i,ie,j,ncases,ipreci

!      write(6,*) ' 0:ordinary randumnumber, 1:doubre precision'
!      read(5,*) ipreci
       ipreci=0

!      write(6,*) ' 0.0:test near 0, 0.1-1:test near specified value'
!      read(5,*) around
      around=1.0

!     ----------
!     Open files
!     ----------
      open(UNIT= 6,FILE='egs5job.out',STATUS='unknown')

      if(ipreci.eq.0) then
        write(6,'(A)') ' Single precision'
      else
        write(6,'(A)') ' Double precision'
      end if

!-----------------------------------------------------------------------
! Step 2: pegs5-call
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Step 3: Pre-hatch-call-initialization
!-----------------------------------------------------------------------

!     --------------------------------------------------------
!     Random number seeds.  Must be defined before call hatch.
!     ins (1- 2^31)
!     --------------------------------------------------------
      inseed=1
      luxlev=1

      write(6,'(/A,I12,5X,A)') ' inseed=',inseed,
     *    ' (seed for generating unique sequences of Ranlux)'

!     =============
      call rluxinit   ! Initialize the Ranlux random-number generator
!     =============

!-----------------------------------------------------------------------
! Step 4:  Determination-of-incident-particle-parameters
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Step 5:   hatch-call
!-----------------------------------------------------------------------

! ----------------------------------------------------------
! Print various data associated with each media (not region)
! ----------------------------------------------------------

!-----------------------------------------------------------------------
! Step 6:  Initialization-for-howfar
!-----------------------------------------------------------------------


!-----------------------------------------------------------------------
! Step 7:  Initialization-for-ausgab
!-----------------------------------------------------------------------

      do ie=1,100
        srandoms(ie)=0.D0
        srandom2s(ie)=0.D0
      end do

      ncount=0
      ncases=1000000
      ncases=1000000

!-----------------------------------------------------------------------
! Step 8:  Shower-call
!-----------------------------------------------------------------------

      do j=1,1000
      do i=1,ncases

        call randomset(rnnow)
        if(ipreci.ne.0) then
          call randomset(rnnow1)
          rnnow=rnnow+twom24*rnnow1
        end if

        if(around.ne.0.0) then
!          if(rnnow.ge.around-1e-6.and.rnnow.le.around) then
         if(rnnow.ge.around-1e-6) then
            ncount=ncount+1
            ie=(around-rnnow)/1e-8+1
            if(ie.gt.100) stop
            if(ie.le.10) then
              write(6,601) rnnow,1-rnnow,1-rnnow-twom24,ie
601           format(' rnnow,ie=',G22.12,2G12.5,I5)
            end if
            srandoms(ie)=srandoms(ie)+ 1.d0
            srandom2s(ie)=srandom2s(ie)+1.d0*1.d0
          end if
        else
          if(rnnow.le.1e-6) then
            ncount=ncount+1
            ie=(rnnow)/1e-8+1
            if(ie.gt.100) stop
            if(ie.le.10) then
              write(6,*) 'rnnow, ie=',rnnow,ie
            end if
            srandoms(ie)=srandoms(ie)+ 1.d0
            srandom2s(ie)=srandom2s(ie)+1.d0*1.d0
          end if
        end if
      end do
      end do

!-----------------------------------------------------------------------
! Step 9:  Output-of-results
!-----------------------------------------------------------------------

      write(6,'(/A,I5)') ' Sampled result: ncount =',ncount
      do ie=1,100
        if(around.ne.0.0) then
          rupp=around-(ie-1)*1.0D-8
          rlow=around-ie*1.0E-8
        else
          rupp=(ie-1)*1.0D-8
          rlow=ie*1.0E-8
        end if
        srandoms(ie)=srandoms(ie)/ncount
        srandom2s(ie)=srandom2s(ie)/ncount
        srandomerr=dsqrt((srandom2s(ie)-srandoms(ie)*srandoms(ie))
     *            /ncount)
        write(6,250) rupp,srandoms(ie),srandomerr
250       format(' Upper random (',G15.8,')=',G15.5,'+-',G15.5)
      end do

      stop
      end

!-------------------------last line of main code------------------------

!-------------------------------ausgab.f--------------------------------
! Version:   040830-1800
! Reference: SLAC-730, KEK-2004-5 (Appendix 2)
!-----------------------------------------------------------------------
!23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
! ----------------------------------------------------------------------
! Required subroutine for use with the EGS5 Code System
! ----------------------------------------------------------------------
! A simple AUSGAB to:
!
!   1) Score energy deposition
!   2) Print out stack information
!   3) Print out particle transport information (if switch is turned on)
!

! ----------------------------------------------------------------------

      subroutine ausgab(iarg)


      return
      end

!--------------------------last line of ausgab.f------------------------

!-------------------------------howfar.f--------------------------------
! Version:   040902-1630
! Reference: SLAC-730, KEK-2004-5 (Appendix 2)
!-----------------------------------------------------------------------
!23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
! ----------------------------------------------------------------------
! Required (geometry) subroutine for use with the EGS5 Code System
! ----------------------------------------------------------------------
! This is a 1-dimensional plane geometry. 
! ----------------------------------------------------------------------

      subroutine howfar

      return
      end

!--------------------------last line of howfar.f------------------------
