!*********************************************************************** !***************************** KEK, High Energy Accelerator Research * !***************************** Organization * !** u c n a i c g v _ v s p ** * !***************************** EGS5.0 USER CODE - 09 Oct 2015/1430 * !*********************************************************************** !* This is a general User Code based on the cg geometry scheme. * !*********************************************************************** ! * ! 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-6002 * ! 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 ucnaicgv_vsp.f User Code requires a cg-input file only * ! (e.g., ucnaicgv.data). * ! The following shows the geometry for ucnaicg.data. * ! Input data for CG geometry must be written at the top of data-input * ! file together with material assignment to each region. Cg-data can * ! be checked by CGview. * ! This user code same with ucnaicgv.f except to calculate volume * ! average spectrum inside NaI detector. * ! Use Ranlux random number generator. * !*********************************************************************** ! * ! ----------------------- * ! cg Geometry (ucnaicgv_vsp) * ! ----------------------- * ! * ! * ! R * ! ^ * ! | * ! +----+----+----+--------+------+--- * ! | * ! | Outer vacuum region * ! + +----+----+--------+------+-----+ R=9.41 * ! | | Air | * ! | | +--------------------+ + R=4.41 * ! | | | Al cover | | * ! | + + +-----------+---+ + R=4.31 * ! | | | | Gap | | | * ! + + + + +-------+ + + R=3.81 * ! | | | | | |Quartz | * ! | | | | | NaI | | | * ! 1.253 MeV | | | | | | | | * ! ============>--+----+---+-------+---+-----+--------> Z * ! photons -5.6 -0.6 -0.5 0.0 7.62 8.12 13.12 * ! -5.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 * depe,deltae,spec(3,50),maxpict real*8 depe,deltae,spec integer maxpict !**** real*8 ! Arguments real*8 totke real*8 rnnow,etot real*8 esumt real*8 ! Local variables * availke,avpe,avph,avspe,avspg,avspp,avte,desci2,pefs,pef2s, * adoses,adose2s,avdose,sigdose,rr0,sigpe,sigte,sigph,sigspg, * sigspe,sigspp,tefs,tef2s,wtin,wtsum,xi0,yi0,zi0,dvol,dweight, * elow,eup,emid,rdet,rtcov,rtgap,tcov,tdet,tgap real*8 * phs(50),ph2s(50),specs(3,50),spec2s(3,50) real * tarray(2),tt,tt0,tt1,cputime,etime integer * i,icases,idin,ie,ifti,ifto,ii,iiz,imed,ireg,isam, * izn,nlist,j,k,n,ner,ntype 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, they must be opened ! after getcg etc. 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 !----------------------------------------------------------------------- ! --------------------------------- ! Define media before calling PEGS5 ! --------------------------------- nmed=4 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 ! ============== medarr(1)='NAI ' medarr(2)='AL ' medarr(3)='QUARTZ ' medarr(4)='AIR-AT-NTP ' do j=1,nmed do i=1,24 media(i,j)=medarr(j)(i:i) end do end do chard(1) = 7.62d0 ! automatic step-size control chard(2) = 0.1d0 chard(3) = 0.5d0 chard(4) = 5.0d0 write(6,fmt="('chard =',5e12.5)") (chard(j),j=1,nmed) ! ----------------------------------- ! Run KEK PEGS5 before calling HATCH ! ----------------------------------- write(6,'(A/)') 'PEGS5-call comes next' ! ========== call pegs5 ! ========== !----------------------------------------------------------------------- ! Step 3: Pre-hatch-call-initialization !----------------------------------------------------------------------- write(6,*) '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(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,'(A)') '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) = 1 ! Switches for PE-angle sampling iedgfl(i) = 1 ! K & L-edge fluorescence iauger(i) = 0 ! K & L-Auger iraylr(i) = 0 ! 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,'(/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 !----------------------------------------------------------------------- ! Define initial variables for incident particle normally incident ! on the slab iqin=0 ! Incident particle charge - photons ekein=1.253 ! Incident particle kinetic energy xin=0.0 ! Source position yin=0.0 zin=-5.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 ! pdf data for many source deltae=0.05 ! Energy bin of response !----------------------------------------- ! 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(6,'(A,I5)') ' Stopped in MAIN. irin = ',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(6,'(/A)') ' 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') ! ========== call hatch ! ========== ! ------------------------------ ! Close files (after HATCH call) ! ------------------------------ close(UNIT=KMPI) close(UNIT=KMPO) ! ---------------------------------------------------------- ! Print various data associated with each media (not region) ! ---------------------------------------------------------- write(6,'(/A)') ' Quantities associated with each MEDIA:' do j=1,nmed write(6,'(/1X,24A1)') (media(i,j),i=1,24) write(6,'(5X,A,G15.7,A,G15.7,A)') * ' rho=',rhom(j),' g/cu.cm rlc=',rlcm(j),' cm' write(6,'(5X,A,G15.7,A,G15.7,A)') * ' ae=',ae(j),' MeV ue=',ue(j),' MeV' write(6,'(5X,A,G15.7,A,G15.7,A/)') * ' ap=',ap(j),' MeV up=',up(j),' MeV' end do ! ------------------------------------------------------- ! Print media and cutoff energies assigned to each region ! ------------------------------------------------------- do i=1,nreg if (med(i) .eq. 0) then write(6,'(A,I3,A)') ' medium(',i,')=vacuum' else write(6,'(A,I3,A,24A1,A,G10.5,A,G10.5,A)') * ' medium(',i,')=',(media(ii,med(i)),ii=1,24), * 'ecut=',ecut(i),' MeV, pcut=',pcut(i),' 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(6,'(A,I3)') ' X-ray information for Z=',izn write(6,'(A/4G15.5/4G15.5/2G15.5)') * ' K-X-ray energy in keV',(ekx(ii,izn),ii=1,10) write(6,'(A/4G15.5/4G15.5)') * ' L-1 X-ray in keV',(elx1(ii,izn),ii=1,8) write(6,'(A/4G15.5/G15.5)') * ' L-2 X-ray in keV',(elx2(ii,izn),ii=1,5) write(6,'(A/4G15.5/3G15.5)') * ' L-3 X-ray in keV',(elx3(ii,izn),ii=1,7) 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(6,'(/A/6X,6(A,14X),A,11X,A,3X,A,1X,A/)') * ' Energy/Coordinates/Direction cosines/etc.', * 'e','x','y','z','u','v','w','iq','ir','iarg' ! Energy bin width deltae=ekein / 50 ! Zero the variables depe=0.D0 pefs=0.D0 pef2s=0.D0 tefs=0.D0 tef2S=0.D0 adoses=0.d0 adose2s=0.d0 do j=1,50 phs(j)=0.D0 ph2s(j)=0.D0 do ntype=1,3 spec(ntype,j)=0.D0 specs(ntype,j)=0.D0 spec2s(ntype,j)=0.D0 end do end do ! Set histories ncases=100000 ! 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 ! ------------------------- ! ---------------------- ! Select incident energy ! ---------------------- 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 ! ---------------------- ! Select incident angle ! ---------------------- ! --------------------------------------------------- ! Print first NWRITE or NLINES, whichever comes first ! --------------------------------------------------- if (ncount .le. nwrite .and. ilines .le. nlines) then ilines = ilines + 1 write(6,'(7G15.7,3I5)') * etot,xin,yin,zin,uin,vin,win,iqin,irin,idin end if ! ----------------------------------------------------------- ! Compare maximum energy of material data and incident energy ! ----------------------------------------------------------- if(etot+(1-iabs(iqin))*RM.gt.emaxe) then write(6,'(A,A)') ' 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(6,'(A,3E12.5)') * 'Following source direction cosines are not normalized.', * 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. if(depe .gt. 0.D0) then adoses=adoses+depe adose2s=adose2s+depe*depe ie=depe/deltae + 1 if (ie .gt. 50) ie = 50 phs(ie)=phs(ie)+wtin ph2s(ie)=ph2s(ie)+wtin*wtin tefs=tefs + wtin tef2s=tef2s + wtin*wtin if(depe .ge. ekein*0.999) then pefs=pefs +wtin pef2s=pef2s +wtin*wtin end if depe = 0.D0 end if do ntype=1,3 do ie=1,50 specs(ntype,ie)=specs(ntype,ie)+spec(ntype,ie) spec2s(ntype,ie)=spec2s(ntype,ie)+ * spec(ntype,ie)*spec(ntype,ie) spec(ntype,ie)=0.D0 end do end do 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,'(A)') '9' ! Set end of batch for CG View tt=etime(tarray) tt1=tarray(1) cputime=tt1-tt0 write(6,'(/A,G15.5)') ' Elapsed Time (sec)=',cputime !----------------------------------------------------------------------- ! Step 9: Output-of-results !----------------------------------------------------------------------- write(6,'(/2(A,I10,A/),A,G15.5,A)') * ' Ncount=',ncount,' (actual cases run)', * ' Ncases=',ncases,' (number of cases requested)', * ' TotKE =',totke,' (total KE (MeV) in run)' if (totke .le. 0.D0) then write(6,'(//A,G15.5/A,G15.5/A,I10)') * ' Stopped in MAIN with TotKE=',totke, * ' AvailKE=',availke,' Ncount=',ncount stop end if ! Detector size. Must be changed for different size of detector. tdet=7.62 rdet=3.81 tcov=0.1 rtcov=0.1 tgap=0.5 rtgap=0.5 write(6,'(2(A,G15.5,A/),4(A,G10.2,A/))') * ' Detector length=',tdet,' cm',' Detector radius=',rdet,' cm', * ' Al cover thickness=',tcov,' cm', * ' Al cover side thickness=',rtcov,' cm', * ' Front gap =',tgap,' cm',' Side gap =',rtgap,' cm' write(6,'(A,G15.5,A/)') ' Results for ',ekein,'MeV photon' ! ----------------------------------- ! Calculate average and its deviation ! ----------------------------------- ! --------------- ! Peak efficiency ! --------------- avpe = pefs/ncount pef2s=pef2s/ncount sigpe=dsqrt((pef2s-avpe*avpe)/ncount) avpe = avpe*100.0 sigpe = sigpe*100.0 write(6,'(A,G11.4,A,G9.2,A)') * ' Peak efficiency =',avpe,'+-',sigpe,' %' ! ---------------- ! Total efficiency ! ---------------- avte = tefs/ncount tef2s = tef2s/ncount sigte = dsqrt((tef2s-avte*avte)/ncount) avte = avte*100.0 sigte = sigte*100.0 write(6,'(A,G11.4,A,G9.2,A)') * ' Total efficiency =', avte,'+-',sigte,' %' ! ------------------------- ! Absorbed dose of detector ! ------------------------- ! 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). ! MeV/g:(absorbed energy in Mev)/(weight) ! 1MeV=1.602E-13J, 1kg=1000g ! 1MeV/g=1.602E-13(J/MeV)*1000(g/kg)=1.602E-10 Gy ! Detector volume and weight dvol=pi*rdet*rdet*tdet dweight=dvol*rhom(1) avdose = adoses/ncount adose2s = adose2s/ncount sigdose = dsqrt((adose2s-avdose*avdose)/ncount) avdose = avdose*1.602E-10/dweight sigdose = sigdose*1.602E-10/dweight write(6,'(/A,G11.4,A,G9.2,A)') * ' Absorbed dose of detector =', avdose,'+-',sigdose, * ' Gy per incident' ! -------------------------- ! Pulse height distribution ! -------------------------- write(6,'(/A)') ' Pulse height distribution ' do ie=1,50 elow=deltae*(ie-1) eup=deltae*ie emid=(elow+eup)/2.d0 avph = phs(ie)/ncount ph2s(ie)=ph2s(ie)/ncount sigph=dsqrt((ph2s(ie)-avph*avph)/ncount) avph = avph/deltae sigph= sigph/deltae write(6,'(A,G10.4,A,G15.5,A,G15.5,A)') * ' E (middle --',emid,' MeV )=',avph,'+-',sigph, * ' counts/MeV/incident' end do ! -------------------------------------------------- ! Volume average particle spectrum inside detector. ! -------------------------------------------------- write(6,'(/A/30X,A/A,11X,A,18X,A,14X,A)') * ' Volume average particle spectrum inside the detector', * 'particles/MeV/cm2 per source photon',' Upper energy', * ' Gamma',' Electron',' Positron' do ie=1,50 elow=deltae*(ie-1) eup=deltae*ie emid=(elow+eup)/2.d0 ! Volume average spectrum inside detector can be calculated from ! track length distribution inside detector divided by its volume. ! ---------------------------------- ! Gamma spectrum per MeV per source ! ---------------------------------- avspg = specs(1,ie)/ncount spec2s(1,ie)=spec2s(1,ie)/ncount sigspg=dsqrt((spec2s(1,ie)-avspg*avspg)/ncount) avspg=avspg/dvol/deltae sigspg= sigspg/dvol/deltae ! ------------------------------------- ! Electron spectrum per MeV per source ! ------------------------------------- avspe = specs(2,ie)/ncount spec2s(2,ie)=spec2s(2,ie)/ncount sigspe=dsqrt((spec2s(2,ie)-avspe*avspe)/ncount) avspe= avspe/dvol/deltae sigspe= sigspe/dvol/deltae ! ------------------------------------ ! Positron spectrum per MeV per source ! ------------------------------------ avspp = specs(3,ie)/ncount spec2s(3,ie)=spec2s(3,ie)/ncount sigspp=dsqrt((spec2s(3,ie)-avspp*avspp)/ncount) avspp= avspp/dvol/deltae sigspp= sigspp/dvol/deltae write(6,'(G10.5,A,3(G12.5,A,G12.5))') * emid,' MeV --',avspg,' +-',sigspg,avspe,' +-', * sigspe,avspp,' +-',sigspp end do nlist=1 ! ============================= call ecnsv1(nlist,nreg,totke) call ntally(nlist,nreg) ! ============================= ! ==================== call counters_out(1) ! ==================== 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 energy deposition and tack length of each particle inside detector ! 3) Print out particle transport information ! 4) call plotxyz if ncount<=maxpict ! ---------------------------------------------------------------------- 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 * depe,deltae,spec(3,50),maxpict real*8 depe,deltae,spec integer maxpict integer ! Arguments * iarg real*8 ! Local variables * edepwt integer * ie,iql,irl,ntype ! ------------------------ ! 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 .ge. 5) return ! ---------------------------------------------- ! Score energy deposition inside NaI detector ! ---------------------------------------------- if (med(irl) .eq. 1) then depe = depe + edepwt ! ---------------------------------------------- ! Score track lenghth (ustep) inside detector ! ---------------------------------------------- if (iql .eq. 0) then ! photon ntype=1 ie = e(np)/deltae +1 if(ie .gt. 50) ie = 50 elseif (iql .eq. -1) then ! electron ntype=2 ie = (e(np) - RM)/deltae +1 if(ie .gt. 50) ie = 50 else ! positron ntype=3 ie = (e(np) - RM)/deltae +1 if(ie .gt. 50) ie = 50 end if spec(ntype,ie) = spec(ntype,ie) + ustep*wt(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,'(7G15.7,3I5)') e(np),x(np),y(np),z(np),u(np),v(np), * w(np),iql,irl,iarg 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(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---------------------