!***********************************************************************
!*****************************  KEK, High Energy Accelerator Research  *
!*****************************  Organization                           *
!*********** d e t ***********                                         *
!*****************************    EGS5.0 USER CODE -  09 Aug 2024/1030 *
!***********************************************************************
!     A frame of EGS5 CG calculation ver. 25Jul2024                    *
!***********************************************************************
!                                                                      *
!  PROGRAMMERS:  H. Hirayama                                           *
!                Applied Research 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-5451                          *
!                Fax:         +81-29-864-4051                          *
!                                                                      *
!                Y. Namito                                             *
!                Radiation Science Center                              *
!                Applied Research Laboratory                           *
!                KEK, High Energy Accelerator Research Organization    *
!                1-1, Oho, Tsukuba, Ibaraki, 305-0801                  *
!                Japan                                                 *
!                                                                      *
!                E-mail:      yoshihito.namito@kek.jp                  *
!                Telephone:   +81-29-864-5489                          *
!                Fax:         +81-29-864-1993                          *
!                                                                      *
!***********************************************************************
!***********************************************************************
!***********************************************************************
!                                                                      *
!             -----------------------                                  *
!             cg Geometry (dat0)                                       *
!             -----------------------                                  *
!                                                                      *
!                                                                      *
!                              R                                       *
!                              ~                                       *
!                              |                                       *
!                              |                                       *
!                 Outer vacuum region                                  *
!                     +----+--------+------+-----+ R=15.0 cm           *
!                     |    Outer Gap (vacuume)   |                     *
!                     |        +-------+         +R=10.0 cm            *
!                     |        |       |         |                     *
!                     |        | NaI   |         |                     *
!                     |        |       |         |                     *
!                     |        |       |         |                     *
! 1.00 MEVV           |        |       |         |                     *
!      ======================>-+-------+---+-----+--------> Z          *
!     photons       -10.0     0.0    5.0        10.0                   *
!                                                                      *
!                                                                      *
!                                                                      *
!                                                                      *
!***********************************************************************
!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/egs5_bounds.f'
      include 'include/egs5_brempr.f'
      include 'include/egs5_edge.f'
      include 'include/egs5_media.f'
      include 'include/egs5_misc.f'
      include 'include/egs5_thresh.f'
      include 'include/egs5_uphiot.f'
      include 'include/egs5_useful.f'
      include 'include/egs5_usersc.f'
      include 'include/egs5_userxt.f'
      include 'include/randomm.f'

!     ----------------------
!     Auxiliary-code COMMONs
!     ----------------------
      include 'auxcommons/aux_h.f'   ! Auxiliary-code "header" file

      include 'auxcommons/edata.f'
      include 'auxcommons/etaly1.f'
      include 'auxcommons/instuf.f'
      include 'auxcommons/lines.f'
      include 'auxcommons/nfac.f'
      include 'auxcommons/watch.f'

!     ------------------
!     cg related COMMONs
!     ------------------
      include 'auxcommons/geom_common.f' ! geom-common file
      integer irinn

      common/totals/                          ! Variables to score
     * maxpict,escore,de
      integer maxpict
      real*8 escore(0:3),de

      real*8 emin,emax,debin,e1,e2,dehist(10000)
      integer nbin,ie

!**** real*8                                                 ! Arguments
      real*8 totke
      real*8 rnnow,etot
      real*8 esumt
      
      real*8                                           ! Local variables
     * availke,wtin,wtsum

      real
     * tarray(2),tt,tt0,tt1,cputime,etime

      integer
     * i,idin,ifti,ifto,ii,iiz,imed,ireg,isam,
     * izn,nlist,j,k,n,ner

      character*24 medarr(MXMED)

      integer NELEM
      parameter(NELEM = 20)
      character*24 denarr(MXMED)
      character*2  elem(MXMED,NELEM)
      character*4  subst(MXMED)   ! substance (ELEM, COMP, or MIXT)
      integer IRAYL,IAPRIM,IBOUND,INCOH,ICPROF,IMPACT,IGSDST
      real*8  rho(MXMED),gasp(MXMED)
      integer ne(MXMED)

!     ----------
!     Open files
!     ----------
!----------------------------------------------------------------
!     Units 7-31 are used in pegs and closed.  It is better not
!     to use as output file. If they are used, they must be opened
!     after getcg etc. Unit for pict must be 39.
!----------------------------------------------------------------

      open(4,FILE='egs5job.inp',STATUS='old')

      ifti = 4
! Write CG input data
      write(ifti,*)"RCC", 1, 0.,0.,0., 0.,0.,5, 10
      write(ifti,*)"RCC", 2, 0.,0.,-10, 0.,0., 20, 15
      write(ifti,*)"END"
      write(ifti,*)"Z1 +1"
      write(ifti,*)"Z2 +2 -1"
      write(ifti,*)"Z3 -2"
      write(ifti,*)"END"
      write(ifti,*)"1 2 0"
      rewind(ifti)

      nmed=2

      subst(1)='COMP'
      medarr(1)='NAI                     '
      denarr(1)=medarr(1)
      rho(1) = 3.667        ! DENSITY (g/cm^3)
      ne(1) = 2
      pz(1,1) = 1.0         ! ATOM FRACTIONs
      pz(1,2) = 1.0
      elem(1,1)='NA'
      elem(1,2)='I '

      subst(2)='MIXT'
      medarr(2) ='AIR AT NTP              '
      denarr(2) ='AIR-GAS                 '
      rho(2) = 1.2050d-3
      gasp(2) = 1.0
      ne(2) = 3
      rhoz(2,1) = 0.75575 
      rhoz(2,2) = 0.23143
      rhoz(2,3) = 0.01282
      elem(2,1)='N '
      elem(2,2)='O '
      elem(2,3)='AR'

      open(39,FILE='egs5job.pic',STATUS='unknown')
      open(95,FILE='dehist.out',STATUS='unknown')

!     ====================
      call counters_out(0)
!     ====================

!-----------------------------------------------------------------------
! Step 2: pegs5-call
!-----------------------------------------------------------------------
!     Prepare calculation
      do i=1,nmed
        AE(i)=0.521
        AP(i)=0.01
        UE(i)=2.511
        UP(i)=2.0
      end do

!     ---------------------------------
!     Define media before calling PEGS5
!     ---------------------------------

      if(nmed.gt.MXMED) then
        write(99,'(A,I4,A,I4,A/A)')  
     *     ' nmed (',nmed,') larger than MXMED (',MXMED,')',
     *     ' MXMED in iclude/egs5_h.f must be increased.'
        stop
      end if

!     ==============
      call block_set                 ! Initialize some general variables
!     ==============

      medarr(1)='PLASTIC SCINT.          '
      medarr(2)='AIR AT NTP              '

      do j=1,nmed
        do i=1,24
          media(i,j)=medarr(j)(i:i)
        end do
      end do  

      chard(1) = 5.0d0       !  automatic step-size control
      chard(2) = 5.0d0      

      write(99,fmt="('chard =',5e12.5)") (chard(j),j=1,nmed)

!     ~~~~~~~ FOR PEGS ~~~~~~~
      gasp(1)=0
      gasp(2)=1
! ~~~~~~~ Pegs COPTIONs INPUT ~~~~~~~
      IRAYL  = 1  
      IAPRIM = 1  
      IBOUND = 0  
      INCOH  = 0  
      ICPROF = 0 
      IMPACT = 0  

      open(25,FILE='pgs5job.pegs5inp')

      do ii = 1, nmed
         write(25,1)subst(ii)
         write(25,1)"&INP"
         if(subst(ii).eq.'COMP')then
            write(25,1)"RHO="
            write(25,2) rho(ii)
            write(25,1)"NE="
            write(25,3) ne(ii)
            write(25,1)"PZ="
            write(25,2)(pz(ii,i),i=1,ne(ii))
         elseif(subst(ii).eq.'MIXT')then
            write(25,1)"RHO="
            write(25,2) rho(ii)
            write(25,1)"NE="
            write(25,3) ne(ii)
            write(25,1)"RHOZ="
            write(25,2)(rhoz(ii,i),i=1,ne(ii))
         endif
         if(gasp(ii).ne.0)then
            write(25,1)"GASP="
            write(25,2) gasp(ii)
         endif
         write(25,1)"IRAYL="
         write(25,3)IRAYL
         write(25,1)"IAPRIM="
         write(25,3)IAPRIM
         write(25,1)"IBOUND="
         write(25,3)IBOUND
         write(25,1)"INCOH="
         write(25,3)INCOH 
         write(25,1)"ICPROF="
         write(25,3)ICPROF
         write(25,1)"IMPACT="
         write(25,3)IMPACT
         write(25,1)"&END"

         write(25,4)medarr(ii),denarr(ii)
         write(25,5)(elem(ii,i),i=1,ne(ii))

         write(25,1)"ENER"
         write(25,1)"&INP"
         write(25,1)"AE="
         write(25,2)AE(ii)
         write(25,1)"AP="
         write(25,2)AP(ii)
         write(25,1)"UE="
         write(25,2)UE(ii)
         write(25,1)"UP="
         write(25,2)UP(ii)
         write(25,1)"&END"
         write(25,1)"PWLF"
         write(25,*)" &INP  &END"
         write(25,1)"DECK"
         write(25,*)" &INP  &END"
      enddo
1       format(a)
2       format(50(1pe10.4,1X))
3       format(i3)
4       format(a24,6x,a24)
5       format(50(a2,1x))

      close(25)
!     -----------------------------------
!     Run KEK PEGS5 before calling HATCH
!     -----------------------------------
      write(99,100)
100   FORMAT('PEGS5-call comes next'/)

!     ==========
      call pegs5
!     ==========

!-----------------------------------------------------------------------
! Step 3: Pre-hatch-call-initialization
!-----------------------------------------------------------------------
      write(99,*) 'Read cg-related data'

!-----------------------------------------------
!     Initialize CG related parameters 
!-----------------------------------------------
      npreci=3     ! PICT data mode for CGView in free format

      ifti = 4     ! Input unit number for cg-data
      ifto = 39    ! Output unit number for PICT

      write(99,fmt="(' CG data')")
      call geomgt(ifti,6)  ! Read in CG data
      write(99,fmt="(' End of CG data',/)")

      if(npreci.eq.3) write(ifto,fmt="('CSTA-FREE-TIME')")
      if(npreci.eq.2) write(ifto,fmt="('CSTA-TIME')")

      rewind ifti
      call geomgt(ifti,ifto)! Dummy call to write geom info for ifto
      write(ifto,110)
110   FORMAT('CEND')

!--------------------------------
!     Get nreg from cg input data
!--------------------------------
      nreg=izonin

!   Read material for each region from egs5job.data
      read(4,*) (med(i),i=1,nreg)

!   Set option except vacuum region
      do i=1,nreg-1
        if(med(i).ne.0) then
          iphter(i) = 0    ! Switches for PE-angle sampling
          iedgfl(i) = 1    ! K & L-edge fluorescence
          iauger(i) = 0    ! K & L-Auger
          iraylr(i) = 1    ! Rayleigh scattering
          lpolar(i) = 0    ! Linearly-polarized photon scattering
          incohr(i) = 0    ! S/Z rejection
          iprofr(i) = 0    ! Doppler broadening
          impacr(i) = 0    ! Electron impact ionization
        end if
      end do

!     --------------------------------------------------------
!     Random number seeds.  Must be defined before call hatch
!     or defaults will be used.  inseed (1- 2^31)
!     --------------------------------------------------------
      luxlev = 1
      inseed=1
      write(99,120) inseed
120   FORMAT(/,' inseed=',I12,5X,
     *         ' (seed for generating unique sequences of Ranlux)')

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

!-----------------------------------------------------------------------
! Step 4:  Determination-of-incident-particle-parameters
!-----------------------------------------------------------------------
! Define initial variables for incident particle normally incident
! on the slab
      iqin=0             ! Incident particle charge - photons
      ekein=1.0          ! Incident particle kinetic energy
      xin=0.0            ! Source position
      yin=0.0
      zin=-1.0
      uin=0.0            ! Moving along z axis
      vin=0.0
      win=1.0
      irin=0             ! Starting region (0: Automatic search in CG)
      wtin=1.0           ! Weight = 1 since no variance reduction used

!-----------------------------------------
!     Get source region from cg input data
!-----------------------------------------
!
      if(irin.le.0.or.irin.gt.nreg) then
        call srzone(xin,yin,zin,iqin+2,0,irin)
        if(irin.le.0.or.irin.ge.nreg) then
          write(99,fmt="(' Stopped in MAIN. irin = ',i5)")irin
          stop
        end if
        call rstnxt(iqin+2,0,irin)
      end if

!-----------------------------------------------------------------------
! Step 5:   hatch-call
!-----------------------------------------------------------------------
      emaxe = 0.D0 ! dummy value to extract min(UE,UP+RM).

      write(99,130)
130   format(/' Call hatch to get cross-section data')

!     ------------------------------
!     Open files (before HATCH call)
!     ------------------------------
      open(UNIT=KMPI,FILE='pgs5job.pegs5dat',STATUS='old')
      open(UNIT=KMPO,FILE='egs5job.dummy',STATUS='unknown')

      write(99,140)
140   FORMAT(/,' HATCH-call comes next',/)

!     ==========
      call hatch
!     ==========

!     ------------------------------
!     Close files (after HATCH call)
!     ------------------------------
      close(UNIT=KMPI)
      close(UNIT=KMPO)

! ----------------------------------------------------------
! Print various data associated with each media (not region)
! ----------------------------------------------------------
      write(99,150)
150   FORMAT(/,' Quantities associated with each MEDIA:')
      do j=1,nmed
        write(99,160) (media(i,j),i=1,24)
160     FORMAT(/,1X,24A1)
        write(99,170) rhom(j),rlcm(j)
170     FORMAT(5X,' rho=',G15.7,' g/cu.cm     rlc=',G15.7,' cm')
        write(99,180) ae(j),ue(j)
180     FORMAT(5X,' ae=',G15.7,' MeV    ue=',G15.7,' MeV')
        write(99,190) ap(j),up(j)
190     FORMAT(5X,' ap=',G15.7,' MeV    up=',G15.7,' MeV',/)
      end do

! -------------------------------------------------------
! Print media and cutoff energies assigned to each region
! -------------------------------------------------------
      do i=1,nreg
        if (med(i) .eq. 0) then
          write(99,200) i
200       FORMAT(' medium(',I3,')=vacuum')
        else
          write(99,210) i,(media(ii,med(i)),ii=1,24),ecut(i),pcut(i)
210       FORMAT(' medium(',I3,')=',24A1,
     *           'ecut=',G10.5,' MeV, pcut=',G10.5,' MeV')
!        -----------------------------------------------
!        Print out energy information of K- and L-X-rays
!        -----------------------------------------------
          if (iedgfl(i) .ne. 0) then             ! Output X-ray energy
            ner = nne(med(i))
            do iiz=1,ner
              izn = zelem(med(i),iiz)  ! Atomic number of this element
              write(99,220) izn
220           FORMAT('   X-ray information for Z=',I3)
              write(99,230) (ekx(ii,izn),ii=1,10)
230           FORMAT('   K-X-ray energy in keV',/,
     *               4G15.5,/,4G15.5,/,2G15.5)
              write(99,240) (elx1(ii,izn),ii=1,8)
240           FORMAT('   L-1 X-ray in keV',/,4G15.5,/,4G15.5)
              write(99,250) (elx2(ii,izn),ii=1,5)
250           FORMAT('   L-2 X-ray in keV',/,5G15.5)
              write(99,260) (elx3(ii,izn),ii=1,7)
260           FORMAT('   L-3 X-ray in keV',/,4G15.5,/,3G15.5)
            end do
          end if
        end if
      end do

      write(39,fmt="('MSTA')")
      write(39,fmt="(i4)") nreg
      write(39,fmt="(15i4)") (med(i),i=1,nreg)
      write(39,fmt="('MEND')")

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

      ncount = 0
      ilines = 0
      nwrite = 10
      nlines = 10
      idin = -1
      totke = 0.
      wtsum = 0.
      iwatch=0

!     =========================
      call ecnsv1(0,nreg,totke)
      call ntally(0,nreg)
!     =========================

      write(99,270)
270   FORMAT(/,' Energy/Coordinates/Direction cosines/etc.',/,
     *        6X,'e',14X,'x',14X,'y',14X,'z',
     *        14X,'u',14X,'v',14X,'w',11X,'iq',3X,'ir',1X,'iarg',/)

!     Set histories
      ncases=10000
!     Set maximum number for pict
      maxpict=50

      tt=etime(tarray)
      tt0=tarray(1)

!-----------------------------------------------------------------------
! Step 8:  Shower-call
!-----------------------------------------------------------------------
!     Write batch number 
      write(39,fmt="('0    1')")

!                     ========================
      if(iwatch.gt.0) call swatch(-99,iwatch)
!                     ========================

                                             ! -------------------------
      do i=1,ncases                          ! Start of shower call-loop
                                             ! -------------------------
        de=0.d0

!       ----------------------
!       Select incident energy
!       ----------------------

        wtin = 1.0

        wtsum = wtsum + wtin               ! Keep running sum of weights
        etot = ekein + iabs(iqin)*RM        ! Incident total energy (MeV)
        escore(0) = escore(0) + etot
        if(iqin.eq.1) then              ! Available K.E. (MeV) in system 
          availke = ekein + 2.0*RM       ! for positron
        else                            ! Available K.E. (MeV) in system
          availke = ekein                ! for photon and electron
        end if

        totke = totke + availke                 ! Keep running sum of KE

!       ----------------------
!       Select incident angle
!       ----------------------

!       ---------------------------------------------------
!       Print first NWRITE or NLINES, whichever comes first
!       ---------------------------------------------------
        if (ncount .le. nwrite .and. ilines .le. nlines) then
          ilines = ilines + 1
          write(99,280) etot,xin,yin,zin,uin,vin,win,iqin,irin,idin
280       FORMAT(7G15.7,3I5)
        end if

!       -----------------------------------------------------------
!       Compare maximum energy of material data and incident energy
!       -----------------------------------------------------------
        if(etot+(1-iabs(iqin))*RM.gt.emaxe) then 
          write(99,fmt="(' Stopped in MAIN.',
     1    ' (Incident kinetic energy + RM) > min(UE,UP+RM).')")
          stop
        end if

!       ----------------------------------------------------
!       Verify the normalization of source direction cosines
!       ----------------------------------------------------
        if(abs(uin*uin+vin*vin+win*win-1.0).gt.1.e-6) then
          write(99,fmt="(' Following source direction cosines are not',
     1    ' normalized.',3e12.5)")uin,vin,win
          stop
        end if

!       ==========================================================
        call shower (iqin,etot,xin,yin,zin,uin,vin,win,irin,wtin)
!       ==========================================================

!      If some energy is deposited inside detector add pulse-height
!      and efficiency.

        ncount = ncount + 1         ! Count total number of actual cases

!                       ========================
        if(iwatch.gt.0) call swatch(-1,iwatch)
!                       ========================
! HISTOGRAM dE
        nbin = 256
        emin = 0
        emax = ekein+abs(iqin)*2*RM
        debin = (emax - emin)/real(nbin - 1)

        if (de .gt. 0d0 ) then
          ie = de/debin +1
          if (ie .gt. nbin) ie = nbin
          dehist(ie) = dehist(ie) + wtin
        endif

                                             ! -----------------------
      end do                                 ! End of CALL SHOWER loop
                                             ! -----------------------
!/////////////////////////..
!     POST PROCESSING
!/////////////////////////..
! WRITE histogram
      e1 = emin
      e2 = e1 + debin
      do i=1,nbin+1
        write(95,'(3G15.4)') e1, e2, dehist(i)/real(ncases)
        e1 = e2
        e2 = e1 + debin
      enddo

      write(6,*)"escore results"
      write(6,*)escore(0)
      write(6,*)escore(1)
      write(6,*)escore(2)
      write(6,*)escore(3)

!                     ========================
      if(iwatch.gt.0) call swatch(-88,iwatch)
!                     ========================

      call plotxyz(99,0,0,0.D0,0.D0,0.D0,0.D0,0,0.D0,0.D0) 

      write(39,fmt="('9')")          ! Set end of batch for CG View

      tt=etime(tarray)
      tt1=tarray(1)
      cputime=tt1-tt0
      write(99,300) cputime
300   format(' Elapsed Time (sec)=',G15.5)

!-----------------------------------------------------------------------
! Step 9:  Output-of-results
!-----------------------------------------------------------------------
      write(99,310) ncount,ncases,totke
310   FORMAT(/,' Ncount=',I10,' (actual cases run)',/,
     *       ' Ncases=',I10,' (number of cases requested)',/,
     *       ' TotKE =',G15.5,' (total KE (MeV) in run)')

      if (totke .le. 0.D0) then
        write(99,320) totke,availke,ncount
320     FORMAT(//,' Stopped in MAIN with TotKE=',G15.5,/,
     *         ' AvailKE=',G15.5, /,' Ncount=',I10)
        stop
      end if

      nlist=1
      

!     =============================
      call ecnsv1(nlist,nreg,totke)
      call ntally(nlist,nreg)
!     =============================

      stop

      end

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

!-------------------------------ausgab.f--------------------------------
! Version:   080708-1600
! Reference: SLAC-265 (p.19-20, Appendix 2)
!-----------------------------------------------------------------------
!23456789|123456789|123456789|123456789|123456789|123456789|123456789|12

! ----------------------------------------------------------------------
! Required subroutine for use with the EGS5 Code System
! ----------------------------------------------------------------------
! A AUSGAB to:
!
!   1) Score energy deposition
!   2) Score particle information enter to detector from outside
!   3) Print out particle transport information 
!   4) call plotxyz if imode=0

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

      subroutine ausgab(iarg)

      implicit none

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

      include 'include/egs5_epcont.f'    ! COMMONs required by EGS5 code
      include 'include/egs5_misc.f'
      include 'include/egs5_stack.f'
      include 'include/egs5_useful.f'

      include 'auxcommons/aux_h.f'   ! Auxiliary-code "header" file

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

      common/totals/                          ! Variables to score
     * maxpict,escore,de
      integer maxpict
      real*8 escore(0:3),de


      integer                                                ! Arguments
     * iarg

      real*8                                           ! Local variables
     * edepwt
     
      integer 
     * ie,iql,irl

!     ------------------------
!     Set some local variables
!     ------------------------
      irl = ir(np)
      iql = iq(np)
      edepwt = edep*wt(np)

!     -----------------------------------------------------------
!     Keep track of energy deposition (for conservation purposes)
!     -----------------------------------------------------------
      if (iarg .lt. 5) then
        esum(iql+2,irl,iarg+1) = esum(iql+2,irl,iarg+1) + edepwt
        nsum(iql+2,irl,iarg+1) = nsum(iql+2,irl,iarg+1) + 1
      end if

!     -----------------------------------------------------------------
!     Print out particle transport information (if switch is turned on)
!     -----------------------------------------------------------------
!                        ========================
      if (iwatch .gt. 0) call swatch(iarg,iwatch)
!                        ========================

      if ( iarg.le.4 )then
        if( ir(np).eq.1 )then
          escore(1) = escore(1) + edep
          de = de + edep
        elseif( ir(np).eq.2 )then
          escore(2) = escore(2) + edep
        elseif( ir(np).eq.3 )then
          escore(3) = escore(3) + edep
        endif
      endif

      if(iarg .ge. 5) return

!     ----------------------------------------------------------------
!     Print out stack information (for limited number cases and lines)
!     ----------------------------------------------------------------
      if (ncount .le. nwrite .and. ilines .le. nlines) then
        ilines = ilines + 1
        write(99,100) e(np),x(np),y(np),z(np),u(np),v(np),w(np),
     *               iql,irl,iarg
 100    FORMAT(7G15.7,3I5)
      end if

!     ------------------------------------
!     Output particle information for plot
!     ------------------------------------
      if (ncount.le.maxpict) then
        call plotxyz(iarg,np,iq(np),x(np),y(np),z(np),e(np),ir(np),
     *       wt(np),time(np))
      end if

      return

      end

!--------------------------last line of ausgab.f------------------------
!-------------------------------howfar.f--------------------------------
! Version:   070627-1600
! Reference: T. Torii and T. Sugita, "Development of PRESTA-CG 
! Incorporating Combinatorial Geometry in EGS4/PRESTA", JNC TN1410 2002-201,
! Japan Nuclear Cycle Development Institute (2002).
! Improved version is provided by T. Sugita. 7/27/2004
!-----------------------------------------------------------------------
!23456789|123456789|123456789|123456789|123456789|123456789|123456789|12

! ----------------------------------------------------------------------
! Required (geometry) subroutine for use with the EGS5 Code System
! ----------------------------------------------------------------------
! This is a CG-HOWFAR. 
! ----------------------------------------------------------------------

      subroutine howfar
      implicit none
c
      include 'include/egs5_h.f'       ! Main EGS "header" file
      include 'include/egs5_epcont.f'  ! COMMONs required by EGS5 code
      include 'include/egs5_stack.f'
      include 'auxcommons/geom_common.f' ! geom-common file
c
c
      integer i,j,jjj,ir_np,nozone,jty,kno
      integer irnear,irnext,irlold,irlfg,itvlfg,ihitcg
      double precision xidd,yidd,zidd,x_np,y_np,z_np,u_np,v_np,w_np
      double precision tval,tval0,tval00,tval10,tvalmn,delhow
      double precision atvaltmp
      integer iq_np
c
      ir_np = ir(np)
      iq_np = iq(np) + 2
c
      if(ir_np.le.0) then
        write(99,*) 'Stopped in howfar with ir(np) <=0'
        stop
      end if
c
      if(ir_np.gt.izonin) then
        write(99,*) 'Stopped in howfar with ir(np) > izonin'
        stop
      end if
c
      if(ir_np.EQ.izonin) then
        idisc=1
        return
      end if
c
      tval=1.d+30
      itvalm=0
c
c     body check
      u_np=u(np)
      v_np=v(np)
      w_np=w(np)
      x_np=x(np)
      y_np=y(np)
      z_np=z(np)
c
      do i=1,nbbody(ir_np)
        nozone=ABS(nbzone(i,ir_np))
        jty=itblty(nozone)
        kno=itblno(nozone)
c     rpp check
        if(jty.eq.ityknd(1)) then
          if(kno.le.0.or.kno.gt.irppin) go to 190
          call rppcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
c     sph check
        elseif(jty.eq.ityknd(2)) then
          if(kno.le.0.or.kno.gt.isphin) go to 190
          call sphcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
c     rcc check
        elseif(jty.eq.ityknd(3)) then
          if(kno.le.0.or.kno.gt.irccin) go to 190
          call rcccg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
c     trc check
        elseif(jty.eq.ityknd(4)) then
          if(kno.le.0.or.kno.gt.itrcin) go to 190
          call trccg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
c     tor check
        elseif(jty.eq.ityknd(5)) then
          if(kno.le.0.or.kno.gt.itorin) go to 190
          call torcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
c     rec check
        elseif(jty.eq.ityknd(6)) then
          if(kno.le.0.or.kno.gt.irecin) go to 190
          call reccg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
c     ell check
        elseif(jty.eq.ityknd(7)) then
          if(kno.le.0.or.kno.gt.iellin) go to 190
          call ellcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
c     wed check
        elseif(jty.eq.ityknd(8)) then
          if(kno.le.0.or.kno.gt.iwedin) go to 190
          call wedcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
c     box check
        elseif(jty.eq.ityknd(9)) then
          if(kno.le.0.or.kno.gt.iboxin) go to 190
          call boxcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
c     arb check
        elseif(jty.eq.ityknd(10)) then
          if(kno.le.0.or.kno.gt.iarbin) go to 190
          call arbcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
c     hex check
        elseif(jty.eq.ityknd(11)) then
          if(kno.le.0.or.kno.gt.ihexin) go to 190
          call hexcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
c     haf check
        elseif(jty.eq.ityknd(12)) then
          if(kno.le.0.or.kno.gt.ihafin) go to 190
          call hafcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
c     tec check
        elseif(jty.eq.ityknd(13)) then
          if(kno.le.0.or.kno.gt.itecin) go to 190
          call teccg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
c     gel check
        elseif(jty.eq.ityknd(14)) then
          if(kno.le.0.or.kno.gt.igelin) go to 190
          call gelcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
c
c**** add new geometry in here
c
       end if
  190  continue
      end do
c
      irnear=ir_np
      if(itvalm.eq.0) then
        tval0=cgeps1
        xidd=x_np+tval0*u_np
        yidd=y_np+tval0*v_np
        zidd=z_np+tval0*w_np
  310   continue
          if(x_np.ne.xidd.or.y_np.ne.yidd.or.z_np.ne.zidd) goto 320
          tval0=tval0*10.d0
          xidd=x_np+tval0*u_np
          yidd=y_np+tval0*v_np
          zidd=z_np+tval0*w_np
          go to 310
  320   continue
c       write(*,*) 'srzone:1'
        call srzone(xidd,yidd,zidd,iq_np,ir_np,irnext)
c
        if(irnext.ne.ir_np) then
          tval=0.0d0
          irnear=irnext
        else
          tval00=0.0d0
          tval10=10.0d0*tval0
          irlold=ir_np
          irlfg=0
  330     continue
          if(irlfg.eq.1) go to 340
            tval00=tval00+tval10
            if(tval00.gt.1.0d+06) then
              write(99,9000) iq(np),ir(np),x(np),y(np),z(np),
     &                      u(np),v(np),w(np),tval00
 9000 format(' TVAL00 ERROR : iq,ir,x,y,z,u,v,w,tval=',
     &       2I3,1P7E12.5)
              stop
            end if
            xidd=x_np+tval00*u_np
            yidd=y_np+tval00*v_np
            zidd=z_np+tval00*w_np
            call srzold(xidd,yidd,zidd,irlold,irlfg)
            go to 330
  340     continue
c
          tval=tval00
          do j=1,10
            xidd=x_np+tval00*u_np
            yidd=y_np+tval00*v_np
            zidd=z_np+tval00*w_np
c           write(*,*) 'srzone:2'
            call srzone(xidd,yidd,zidd,iq_np,irlold,irnext)
            if(irnext.ne.irlold) then
              tval=tval00
              irnear=irnext
            end if
            tval00=tval00-tval0
          end do
          if(ir_np.eq.irnear) then
            write(0,*) 'ir(np),tval=',ir_np,tval
          end if
        end if
      else
        do j=1,itvalm-1
          do i=j+1,itvalm
            if(atval(i).lt.atval(j)) then
              atvaltmp=atval(i)
              atval(i)=atval(j)
              atval(j)=atvaltmp
            endif
          enddo
        enddo
        itvlfg=0
        tvalmn=tval
        do jjj=1,itvalm
          if(tvalmn.gt.atval(jjj)) then
            tvalmn=atval(jjj)
          end if
          delhow=cgeps2
          tval0=atval(jjj)+delhow
          xidd=x_np+tval0*u_np
          yidd=y_np+tval0*v_np
          zidd=z_np+tval0*w_np
  410     continue
          if(x_np.ne.xidd.or.y_np.ne.yidd.or.z_np.ne.zidd) go to 420
            delhow=delhow*10.d0
            tval0=atval(jjj)+delhow
            xidd=x_np+tval0*u_np
            yidd=y_np+tval0*v_np
            zidd=z_np+tval0*w_np
          go to 410
  420     continue
c         write(*,*) 'srzone:3'
          call srzone(xidd,yidd,zidd,iq_np,ir_np,irnext)
          if((irnext.ne.ir_np.or.atval(jjj).ge.1.).and.
     &        tval.gt.atval(jjj)) THEN
            tval=atval(jjj)
            irnear=irnext
            itvlfg=1
            goto 425
          end if
        end do
  425   continue
        if(itvlfg.eq.0) then
          tval0=cgmnst
          xidd=x_np+tval0*u_np
          yidd=y_np+tval0*v_np
          zidd=z_np+tval0*w_np
  430     continue
          if(x_np.ne.xidd.or.y_np.ne.yidd.or.z_np.ne.zidd) go to 440
            tval0=tval0*10.d0
            xidd=x_np+tval0*u_np
            yidd=y_np+tval0*v_np
            zidd=z_np+tval0*w_np
            go to 430
  440     continue
          if(tvalmn.gt.tval0) then
            tval=tvalmn
          else
            tval=tval0
          end if
        end if
      end if
      ihitcg=0
      if(tval.le.ustep) then
        ustep=tval
        ihitcg=1
      end if
      if(ihitcg.eq.1) THEN
        if(irnear.eq.0) THEN
          write(99,9200) iq(np),ir(np),x(np),y(np),z(np),
     &                  u(np),v(np),w(np),tval
 9200 format(' TVAL ERROR : iq,ir,x,y,z,u,v,w,tval=',2I3,1P7E12.5)
          idisc=1
          itverr=itverr+1
          if(itverr.ge.100) then
            stop
          end if
          return
        end if
        irnew=irnear
        if(irnew.ne.ir_np) then
          call rstnxt(iq_np,ir_np,irnew)
        endif
      end if
      return
      end
!--------------------last line of subroutine howfar---------------------
