!***********************************************************************
!*****************************  KEK, High Energy Accelerator Research  *
!*****************************  Organization                           *
!** u c p h a n t o m c g v **                                         *
!*****************************    EGS5.0 USER CODE -  09 Jul 2013/1700 *
!                                 Add ambient dose equivalent          *
!                                 Add explanation of dose calculation  *
!***********************************************************************
!***********************************************************************
!                                                                      *
!  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                          *
!                                                                      *
!***********************************************************************
!***********************************************************************
! The ucphantomcgv.f User Code requires a cg-input file only           *
! (e.g., ucphantomcgv.data).                                           *
! The following shows the geometry for ucphantomcgv.data.              *
! Cg-data can be checked by CGview.                                    *
! This user code corresponds to ucphantomcgp.mor for egs4.             *
! Use Ranlux random number generator.                                  *
!***********************************************************************
!                                                                      *
!             --------------------------                               *
!             cg Geometry (ucphantomcgv)                               *
!             ---------------------------                              *
!                                                                      *
!                                                                      *
!                X                                                     *
!                ^                                                     *
!                |                                                     *
!           +----+----+----+----   -+---+----+----+--- 20.0            *
!           |                                     |                    *
!           |     Outer vacuum region             |                    *
!           +    +----+----+----   -+---+----+    +--- 15.0            *
!           |    |    |                 |    |    |                    *
!           |    +    +  Water (H2O)    |    |    |                    *
!           |    |    |                 |    |    |                    *
!           +    +    +---+---+-   -+---+    +    +--   0.5            *
!           |    |    |   |   |     |   |    |    |                    *
!           |    |Air |H2O|H2O|     |H2O|Air |    |                    *
! 1.253MeV  |    |    |   |   |     |   |    |    |                    *
!      ==========>+---+---+---+-   -+---+----+----+--> Z               *
! photons -20   -5.0  0  1.0 2.0  19.0 20.0     40.0                   *
!                                          25.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_edge.f'
      include 'include/egs5_elecin.f'
      include 'include/egs5_media.f'
      include 'include/egs5_misc.f'
      include 'include/egs5_stack.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
     * depe(20),faexp,fexps,fambde,sambde,maxpict,ndet
      real*8 depe,faexp,fexps,fambde,sambde
      integer maxpict,ndet

!**** real*8                                                 ! Arguments
      real*8 etot,totke
      integer ins

!**** real*8                                           ! Local variables
      real*8
     * area,availke,depthl,depths,dis,disair,ei0,elow,eup,
     * phai0,phai,radma2,rnnow,sinth,sposi,tnum,vol,w0,wimin,wtin,
     * wtsum,xhbeam,xpf,yhbeam,ypf

      real*8 bsfa,bsferr,faexps,faexp2s,faexrr,fexpss,fexps2s,fexerr,
     *       faexpa,fexpsa,fambdes,fambde2s,sambdes,sambde2s,fambdeq,
     *       famberr,sambdeq,samberr

      real*8
     * depeh(20),depeh2(20),dose(20),dose2(20),doseun(20)

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

      integer
     * i,ii,icases,idin,ie,ifti,ifto,imed,ireg,isam,
     * ixtype,j,k,kdet,nnn

      character*24 medarr(MXMED)

!     ----------
!     Open files
!     ----------
!-----------------------------------------------------------------
!     Units 7-26 are used in pegs and closed.  It is better not
!     to use as output file. If they are used must be re-open after
!     call pegs5.   Unit for pict must be 39.
!-----------------------------------------------------------------

      open(6,file='egs5job.out',status='unknown')
      open(4,FILE='egs5job.inp',STATUS='old')
      open(39,FILE='egs5job.pic',STATUS='unknown')

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

!-----------------------------------------------------------------------
! Step 2: pegs5-call
!-----------------------------------------------------------------------
      nmed=2
      if(nmed.gt.MXMED) then
        write(6,'(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
!     ==============

!     ---------------------------------
!     define media before calling PEGS5
!     ---------------------------------

      medarr(1)='WATER                   '
      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) = 1.0d0       !  automatic step-size control
      chard(2) = 1.0d0
      write(6,fmt="('chard =',5e12.5)") (chard(j),j=1,nmed)

!     ------------------------------
!     Run PEGS5 before calling HATCH
!     ------------------------------
      write(6,*) 'PEGS5-call comes next'

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

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

!-----------------------------------------------
!     Define pict data mode.
!-----------------------------------------------
!     npreci 1: for PICT32
!            2: for CGview
!            3: for CGview in free format
      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(6,fmt="(' CG data')")
      call geomgt(ifti,6)  ! Read in CG data
      write(6,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=2,nreg-2
        if(med(i).ne.0) then
          iphter(i) = 1    ! 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(6,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 source position from phantom surface.
!-------------------------------------------------
!     Source position from phantom surface in cm.
      sposi=10.0

      iqin=0             ! Incident charge - photons
      ekein=1.253        ! Kinetic energy of source photon
      etot=ekein + abs(iqin)*RM
      xin=0.D0    
      yin=0.D0    
      zin=-sposi 
      uin=0.D0    
      vin=0.D0    
      win=1.D0    
      irin=0      ! Starting region (0: Automatic search in CG)

!---------------------------------------------
!     Half width and height at phantom surface
!---------------------------------------------
!     X-direction half width of beam at phantom surface in cm.
      xhbeam=1.0
!     Y-direction half height of beam at phantom surface in cm.
      yhbeam=1.0
      radma2=xhbeam*xhbeam+yhbeam*yhbeam
      wimin=sposi/dsqrt(sposi*sposi+radma2)

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

      write(6,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(6,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(6,150)
150   FORMAT(/,' Quantities associated with each MEDIA:')
      do j=1,nmed
        write(6,160) (media(i,j),i=1,24)
160     FORMAT(/,1X,24A1)
        write(6,170) rhom(j),rlcm(j)
170     FORMAT(5X,' rho=',G15.7,' g/cu.cm     rlc=',G15.7,' cm')
        write(6,180) ae(j),ue(j)
180     FORMAT(5X,' ae=',G15.7,' MeV    ue=',G15.7,' MeV')
        write(6,190) ap(j),up(j)
190     FORMAT(5X,' ap=',G15.7,' MeV    up=',G15.7,' MeV',/)
      end do

      write(6,200)
200   FORMAT(/' Information of medium and cut-off for each region')
      do i=1,nreg
        if (med(i).eq.0) then
          write(6,210) i
210       FORMAT(' Medium(region:',I5,')= Vacuum')
        else
          write(6,220) i,(media(ii,med(i)),ii=1,24),
     *                 ecut(i),pcut(i),rhor(i)
220       FORMAT(' Medium(region:',I5,
     *           ')=',24A1,/5X,'ECUT=',G10.5,' MeV, PCUT=',
     *           G10.5, ' MeV, density=',F10.3)
        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 = 25
      idin = -1
      totke = 0.
      wtsum = 0.

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

!---------------------
!     Clear variables
!---------------------
      do nnn=1,20
        depe(nnn)=0.D0
        depeh(nnn)=0.D0
        depeh2(nnn)=0.D0
      end do  

      faexp=0.D0
      faexps=0.D0
      faexp2s=0.D0
      fexps=0.D0
      fexpss=0.D0
      fexps2s=0.D0
      fambde=0.d0
      fambdes=0.d0
      fambde2s=0.d0
      sambde=0.d0
      sambdes=0.d0
      sambde2s=0.d0

!-----------------------------
!     Detector number to score
!-----------------------------
      ndet=20

      write(6,230)
230   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',/)

!--------------------------
!     History number 
!--------------------------
!     History number
      ncases=100000
!     Maximum history number to write trajectory data
      maxpict=100
      iwatch=0

      write(39,fmt="('0    1')")

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

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

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

                                             ! -------------------------
      do j=1,ncases                          ! Start of CALL SHOWER loop
                                             ! -------------------------
!--------------------------------------
!       Determine direction (isotropic)
!--------------------------------------
240     call randomset(w0)
        win=w0*(1.0-wimin)+wimin
        call randomset(phai0)
        phai=pi*(2.0*phai0-1.0)
        sinth=dsqrt(1.D0-win*win)
        uin=dcos(phai)*sinth
        vin=dsin(phai)*sinth
        dis=sposi/win
        xpf=dis*uin
        ypf=dis*vin
        if (dabs(xpf).gt.xhbeam.or.dabs(ypf).gt.yhbeam) go to 240
        if (sposi.gt.5.0) then
          disair=(sposi-5.0)/win
          xin=disair*uin
          yin=disair*vin
          zin=-5.D0
        else
          xin=0.D0
          yin=0.D0
          zin=-sposi
        end if

!       ------------------------------------
!       Get source region from cg input data
!       ------------------------------------

        if(irin.le.0.or.irin.gt.nreg) then
          call srzone(xin,yin,zin,iqin+2,0,irinn)
          if(irinn.le.0.or.irinn.ge.nreg) then
            write(6,fmt="(' Stopped in MAIN. irinn = ',i5)")irinn
            stop
          end if
          call rstnxt(iqin+2,0,irinn)
        else
          irinn=irin
        end if

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

        ekein=ekein
        wtin = 1.0

        wtsum = wtsum + wtin               ! Keep running sum of weights
        etot = ekein + iabs(iqin)*RM       ! Incident total energy (MeV)
        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

        latchi=0

!       ---------------------------------------------------
!       Print first NWRITE or NLINES, whichever comes first
!       ---------------------------------------------------
        if (ncount .le. nwrite .and. ilines .le. nlines) then
          ilines = ilines + 1
          write(6,250) etot,xin,yin,zin,uin,vin,win,iqin,irinn,idin
250       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(6,fmt="(' Stopped in MAIN.',
     1    ' (Incident kinetic energy + RM) > min(UE,UP+RM).')")
          stop
        end if

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

        uf(1)=0.0
        vf(1)=0.0
        wf(1)=0.0  ! Needed if lpolar(i)=1

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

!----------------------------------
!      Sum variable and its square.
!----------------------------------

        do kdet=1,ndet
          depeh(kdet)=depeh(kdet)+depe(kdet)
          depeh2(kdet)=depeh2(kdet)+depe(kdet)*depe(kdet)
          depe(kdet)=0.0
        end do

        faexps=faexps+faexp
        faexp2s=faexp2s+faexp*faexp
        faexp=0.0
        fexpss=fexpss+fexps
        fexps2s=fexps2s+fexps*fexps
        fexps=0.0

        fambdes=fambdes+fambde
        fambde2s=fambde2s+fambde*fambde
        fambde=0.d0
        sambdes=sambdes+sambde
        sambde2s=sambde2s+sambde*sambde
        sambde=0.d0

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

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

                                               ! -----------------------
      end do                                   ! End of CALL SHOWER loop
                                               ! -----------------------

!                     ========================
      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
      close(UNIT=39,status='keep')

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

!-----------------------------------------------------------------------
! Step 9:  Output-of-results
!-----------------------------------------------------------------------
!     ---------------------
!     Write out the results
!     ---------------------
      write(6,280) ncount,ncases,totke,totke/ncount
280   FORMAT(/,' Ncount=',I10,' (actual cases run)',/,
     *       ' Ncases=',I10,' (number of cases requested)',/,
     *       ' TotKE =',G15.5,' (total KE (MeV) in run)'/
     *       ' Average Kinetic enegy =',G15.5,'MeV'/)

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

!-----------------------------
!     Sampled source spectrum
!-----------------------------

      write(6,300) sposi 
300   FORMAT(/' Absorbed energy inside phantom for 1.253MeV photon'/
     *    ' Source position ',F10.1,' cm from phantom surface'/ 
     *    ' Within 1cm x 1 cm area after 5 cm air')

      write(6,310) ncases, xhbeam, yhbeam
310   FORMAT(1X,I8,' photons normally incident from front side'/
     *' Half width of beam is ',G15.5,'cm for X and ',G15.5,'cm for Y')

!----------------------------------------------
!     Calculate average dose and its deviation
!----------------------------------------------
!     Conversion from absorbed energy in MeV to absorbed dose in Gy
!     Ausgab scores absorbed energy in unit of MeV
!     Main routine converts this into Gy (J/kg).
!     vol=area*depth(=1.0): volume of region in cm^3
!     MeV/g:(absorbed energy in MeV)/(vol*density(=1.0))
!     1MeV=1.602E-13J,  1kg=1000g
!     1MeV/g=1.602E-13(J/MeV)*1000(g/kg)=1.602E-10 Gy

      area=1.D0*1.D0
      do kdet=1,ndet 
        vol=area*1.D0
        dose(kdet)=depeh(kdet)/ncases
        dose2(kdet)=depeh2(kdet)/ncases
        doseun(kdet)=dsqrt((dose2(kdet)-dose(kdet)*dose(kdet))/ncases)
        dose(kdet)=dose(kdet)*1.602E-10/vol
        doseun(kdet)=doseun(kdet)*1.602E-10/vol
        depths=kdet-1.0
        depthl=kdet
        write(6,320)depths,depthl,(media(ii,med(kdet+1)),ii=1,24),
     *  rhor(kdet+1),dose(kdet),doseun(kdet)
320     FORMAT(' At ',F4.1,'--',F4.1,'cm (',24A1,',rho:',F8.4,')=',
     *  G13.5,'+-',G13.5,'Gy/incident')
      end do
        

!-------------------------------------------------------------------------
!     Calculate average air absorbed dose in Gy (Air Gy) and its deviation
!-------------------------------------------------------------------------
!     Conversion from air absorbed dose in MeV cm2/g to that in Gy
!     Unit of mass energy absorption coefficient mu_en is cm2/g
!     Ausgab scores energy (MeV) times mu_en in unit of MeV cm2/g.
!     Main routine converts this into Gy (J/kg).
!     1MeV=1.602E-13J,  1kg=1000g
!     1MeV/g=1.602E-13(J/MeV)*1000(g/kg)=1.602E-10 Gy
!     Dividing by detector area (for example, 1cm2).

      faexpa=faexps/ncases
      faexp2s=faexp2s/ncases
      faexrr=dsqrt((faexp2s-faexpa*faexpa)/ncases) 
      faexpa=faexpa*1.602E-10/area
      faexrr=faexrr*1.602E-10/area
      fexpsa=fexpss/ncases
      fexps2s=fexps2s/ncases
      fexerr=dsqrt((fexps2s-fexpsa*fexpsa)/ncases) 
      fexpsa=fexpsa*1.602E-10/area
      fexerr=fexerr*1.602E-10/area
      if (faexpa.gt.0.0) then
        bsfa=fexpsa/faexpa
        bsferr=bsfa*dsqrt((faexrr/faexpa)**2.+(fexerr/fexpsa)**2.)
        write(6,330) faexpa,faexrr,fexpsa,fexerr,bsfa,bsferr
330     FORMAT(/' Exposure in free air (using mu_en) ',7X,'=', 
     *  G15.5,'+-',G15.5,' Gy/incident'/
     *  ' Exposure at phantom surface (using mu_en) =',G15.5,
     *  '+-',G15.5,' Gy/incident'/ ' Backscattering factor =',
     *  G15.5,'+-',G15.5)
      else
        write(6,340) faexpa,faexrr,fexpsa,fexerr
340     FORMAT(/' Exposure in free air (using mu_en) =', G15.5,'+-',
     *  G15.5,' Gy/incident'/
     *  ' Exposure at phantom surface (using mu_en) ='
     *  , G15.5,'+-',G15.5,'Gy/incident')
      end if

!-----------------------------------------------------------------
!     Calculate average ambient dose equivalent and its deviation
!-----------------------------------------------------------------
!     Conversion from ambient dose equivalent in MeV cm2/g to Sv
!     Ausgab scores absorbed energy of air in MeV cm2/g times
!     ratio of ambient dose equivalent, Sv, to air collision kerma in Gy.
!     Main routine obtaine ambient dose equivalent in Sv by
!     converting MeV cm2/g into Gy.

      fambdeq=fambdes/ncases
      fambde2s=fambde2s/ncases
      famberr=dsqrt((fambde2s-fambdeq*fambdeq)/ncases) 
      fambdeq=fambdeq*1.602E-10/area
      famberr=famberr*1.602E-10/area
      sambdeq=sambdes/ncases
      sambde2s=sambde2s/ncases
      samberr=dsqrt((sambde2s-sambdeq*sambdeq)/ncases) 
      sambdeq=sambdeq*1.602E-10/area
      samberr=samberr*1.602E-10/area
      write(6,350) fambdeq,famberr,sambdeq,samberr
350   FORMAT(/' Ambient dose equivalent in free air (using mu_en) ',
     *  7X,'=',G15.5,'+-',G15.5,' Sv/incident'/ 
     *  ' Ambient dose equivalent at phantom surface (using mu_en) =',
     *  G15.5,'+-',G15.5,' Sv/incident')

!     =============================
      call ecnsv1(1,nreg,totke)
!     =============================

!     ====================
      call counters_out(1)
!     ====================

!     -----------
!     Close files
!     -----------
      close(UNIT=1)
      close(UNIT=4)

      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 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)

      implicit none

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

      include 'include/egs5_epcont.f'    ! COMMONs required by EGS5 code
      include 'include/egs5_media.f'
      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
     * depe(20),faexp,fexps,fambde,sambde,maxpict,ndet
      real*8 depe,faexp,fexps,fambde,sambde
      integer maxpict,ndet

      integer                                                ! Arguments
     * iarg

      real*8                                           ! Local variables
     * cmod,dcon,edepwt,encoea,ekein,decoe,decon

      integer idet,ie,iql,irl

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

!  ---------------------------------
!  Keep track of how deep stack gets
!  ---------------------------------
      if (np.gt.MXSTACK) then
        write(6,100) np,MXSTACK
100     FORMAT(//' In AUSGAB, np=',I3,' >= maximum stack',
     *       ' allowed which is',I3/1X,79('*')//)
        stop
      end if

!     ------------------------
!     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
      end if

!-------------------------------------------------
!     Score data at detector region (region 2-21)
!--------------------------------------------------
      if (irl.ge.2.and.irl.le.nreg-3) then
        idet=irl-1
        if(idet.ge.1.and.idet.le.ndet) then
          depe(idet)=depe(idet)+edepwt/rhor(irl)
        end if
      end if

!--------------------------------
!     Check cross phantom surface
!--------------------------------
      if (abs(irl-irold).eq.1.and.iq(np).eq.0) then
        if((w(np).gt.0.0.and.irl.eq.2).or.(w(np).le.0.0.and.irl.eq.1))
     *  then
          if (dabs(w(np)).ge.0.0349) then
            cmod=dabs(w(np))
          else 
            cmod=0.0175
          end if
          ekein=e(np)
          dcon=encoea(ekein)            ! Absorbed energy in air
          decon=decoe(ekein)            ! Sv/Gy for ambient DE
          fexps=fexps+e(np)*dcon*wt(np)/cmod
          sambde=sambde+e(np)*dcon*decon*wt(np)/cmod
          if (w(np).lt.0.0) latch(np)=1
          if (w(np).gt.0.0.and.latch(np).eq.0) then
            faexp=faexp+e(np)*dcon*wt(np)/cmod
            fambde=fambde+e(np)*dcon*decon*wt(np)/cmod
          end if
        end if
      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

!     ----------------------------------------------------------------
!     Print out stack information (for limited number cases and lines)
!     ----------------------------------------------------------------
      if (ncount .le. nwrite .and. ilines .le. nlines) then
        ilines = ilines + 1
        write(6,110) e(np),x(np),y(np),z(np),u(np),v(np),w(np),
     *               iql,irl,iarg
 110    FORMAT(7G15.7,3I5)
      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(6,*) 'Stopped in howfar with ir(np) <=0'
        stop
      end if
c
      if(ir_np.gt.izonin) then
        write(6,*) '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(6,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(6,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---------------------
!-------------------------------encoea.f--------------------------------
! Version:   030831-1300
!-----------------------------------------------------------------------
!23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
!
!     double precision function encoea(energy) 
!     Function to evaluate the energy absorption coefficient of air.    
!      (Tables and Graphs oh photon mass attenuation coefficients and 
!      energy-absorption coefficients for photon energies 1 keV to 
!      20 MeV for elements Z=1 to 92 and some dosimetric materials, 
!      S. M. Seltzer and J. H. Hubbell 1995, Japanese Society of 
!      Radiological Technology)
!----------------------------------------------------------------------
      double precision function encoea(energy)

      real*8 hnu(38)/0.001,0.0015,0.002,0.003,0.0032029,0.0032029,
     *           0.004,0.005,0.006,0.008,0.01,0.015,0.02,0.03,0.04,
     *           0.05,0.06,0.08,0.10,0.15,0.2,0.3,0.4,0.5,0.6,0.8,1.0,
     *           1.25,1.5,2.0,3.0,4.0,5.0,6.0,8.0,10.0,15.0,20.0/

      real*8 enmu(38)/3599.,  1188.,  526.2,  161.4,  133.0,  146.0,  
     *      76.36, 39.31,  22.70,  9.446,  4.742,  1.334,  0.5389,
     *      0.1537,0.06833,0.04098,0.03041,0.02407,0.02325,0.02496,
     *      0.02672,0.02872,0.02949,0.02966,0.02953,0.02882,0.02789,
     *      0.02666,0.02547,0.02345,0.02057,0.01870,0.01740,0.01647,
     *      0.01525,0.01450,0.01353,0.01311/
                                                  
      real*8 energy,enm1,hnu1,ene0,slope

      integer i

      if (energy.gt.hnu(38)) then
        encoea=enmu(38)
        return
      end if
      if (energy.lt.hnu(1)) then
        encoea=enmu(1)
        return
      end if
      
      do i=1,38
        if(energy.ge.hnu(i).and.energy.lt.hnu(i+1)) then
          enm1=dlog(enmu(i+1))
          enm0=dlog(enmu(i))
          hnu1=dlog(hnu(i+1))
          hnu0=dlog(hnu(i))

          ene0=dlog(energy)
          slope=(enm1-enm0)/(hnu1-hnu0)
          encoea=exp(enm0+slope*(ene0-hnu0))
          return
        end if
        if(energy.eq.hnu(i+1)) then
          encoea=enmu(i+1)
          return
        end if
      end do

! If sort/interpolation cannot be made, indicate so by writing
! a comment and stopping here.
      write(6,100) energy
100   FORMAT(///,' *****STOPPED IN ENCOEA*****',/,' E=',G15.5,///)
      return
      end
                                                                       
!--------------------------last line of encoea.f------------------------
!------------------------------decoe.f--------------------------------
! Version:   100302-1000
!-----------------------------------------------------------------------
!23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
! ----------------------------------------------------------------------
! Function to evaluate the ratio of ambient dose equivalent to 
! air collision kerma (Sv/Gy). 
! Data tanken from ICRP pub 74 (1996).
! ----------------------------------------------------------------------
      double precision function decoe(energy)

      implicit none

      real*8 energy, slope
      integer i

      real*8 hnu(25)/
     *  0.01,0.015,0.02,0.03,0.04,0.05,0.06,0.08,
     *  0.10,0.15,0.2,0.3,0.4,0.5,0.6,0.8,1.0,1.5,2.0, 
     *  3.0,4.0,5.0,6.0,8.0,10.0/

      real*8 enmu(25)/0.008,0.26,0.61,1.10,1.47,1.67,1.74,1.72,
     *  1.65,1.49,1.40,1.31,1.26,1.23,1.21,1.19,1.17,1.15,1.14,
     *  1.13,1.12,1.11,1.11,1.11,1.10/

      if(energy.gt.hnu(25)) then
        decoe=enmu(25)
        return
      end if

      if (energy.lt.hnu(1)) then
        decoe=enmu(1)
        return
      end if

      do i=1,25
        if(energy.ge.hnu(i).and.energy.lt.hnu(i+1)) then
          slope=(dlog(enmu(i+1))-dlog(enmu(i)))/
     *          (dlog(hnu(i+1))-dlog(hnu(i)))
          decoe=dlog(enmu(i))+slope*(dlog(energy)-dlog(hnu(i)))
          decoe=exp(decoe)
          return
        end if
        if(energy.eq.hnu(i+1)) then
          decoe=enmu(i+1)
          return
        end if
      end do

! If sort/interpolation cannot be made, indicate so by writing
! a comment and stopping here.

      write(3,100) energy
100   format(///,' **** Stopped in decoe ****',/,' E=',G15.5,///)
      stop

      return
      end

!-------------------------last line of decoe.f------------------------
