!*********************************************************************** !***************************** KEK, High Energy Accelerator Research * !***************************** Organization * !* u c i o n c h _ c g v ***** * !***************************** EGS5.0 USER CODE - 30 May 2021/1430 * ! correct score condition * !*********************************************************************** !* 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-5451 * ! Fax: +81-29-864-4051 * ! * !*********************************************************************** !*********************************************************************** ! The ucionch_cgv.f User Code requires a cg-input file only * ! (e.g., ucionch_cgv.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. * ! Use Ranlux random number generator. * !*********************************************************************** ! * ! ------------------------- * ! cg Geometry (ucionch_cgv) * ! ------------------------- * ! * ! * ! R * ! ^ * ! | * ! +----+----+----+--------+------+--- * ! | * ! | Outer vacuum region * ! + +----+----+--------+------+-----+ R=10.0 * ! | | Air | * ! | | +--------------------+ + R=5.2 * ! | | | Case | | * ! | + + +-----------+ + + R=5.0 * ! | | | | 1 atm | | | * ! | | | | Air | | | * ! | | | | | | | * ! | | | | | | | * ! 0.662 MeV | | | | | | | * ! ==============>+----+-----------+---+-----+--------> Z * ! photons -10.0 -0.2 -0.0 0.0 10.0 10.2 15.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 * dose,dosep,deltae,spec(50),maxpict real*8 dose,dosep,deltae,spec integer maxpict !**** real*8 ! Arguments real*8 totke real*8 rnnow,etot real*8 esumt real*8 ! Local variables * availke,doses,dose2s,doseps,dosep2s,avdose,sigdose,avspe,sigspe, * rr0,wtin,wtsum, xi0,yi0,zi0,vol,weight real*8 * specs(50),spec2s(50) real ! Local variables * elow,eup,rdet,rcase,tcase,tdet 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=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 ! ============== medarr(1)='20 degree 1 atm air ' medarr(2)='POLYETHYLENE ' 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) = 0.2d0 write(6,fmt="('chard =',5e12.5)") (chard(j),j=1,nmed) ! ----------------------------------- ! Run KEK PEGS5 before calling HATCH ! ----------------------------------- write(6,100) 100 FORMAT('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,110) 110 FORMAT('CEND') !-------------------------------- ! Get nreg from cg input data !-------------------------------- nreg=izonin ! Read material for each refion from egs5job.data read(4,*) (med(i),i=1,nreg) ! Set option except vacuum region do i=1,nreg 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,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=0.662 ! 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,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(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 ! ------------------------------------------------------- ! Print media and cutoff energies assigned to each region ! ------------------------------------------------------- do i=1,nreg if (med(i) .eq. 0) then write(6,200) i 200 FORMAT(' medium(',I3,')=vacuum') else write(6,210) i,(media(ii,med(i)),ii=1,24),ecut(i),pcut(i) 210 FORMAT(' medium(',I3,')=',24A1, * 'ecut=',G12.5,' MeV, pcut=',G12.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(6,220) izn 220 FORMAT(' X-ray information for Z=',I3) write(6,230) (ekx(ii,izn),ii=1,10) 230 FORMAT(' K-X-ray energy in keV',/, * 4G15.5,/,4G15.5,/,2G15.5) write(6,240) (elx1(ii,izn),ii=1,8) 240 FORMAT(' L-1 X-ray in keV',/,4G15.5,/,4G15.5) write(6,250) (elx2(ii,izn),ii=1,5) 250 FORMAT(' L-2 X-ray in keV',/,5G15.5) write(6,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(6,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',/) ! Energy bin width deltae=ekein / 50 ! Zero the variables dose=0.D0 doses=0.D0 dose2s=0.D0 dosep=0.D0 doseps=0.D0 dosep2s=0.D0 do j=1,50 spec(j)=0.D0 end do ! Set histories ncases=10000000 ! 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,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(6,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(6,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) ! ========================================================== ! ---------------------------- ! Sum variable and its square. ! ---------------------------- doses=doses+dose dose2s=dose2s+dose*dose dose=0.d0 doseps=doseps+dosep dosep2s=dosep2s+dosep*dosep dosep=0.d0 do ie=1,50 specs(ie)=specs(ie)+spec(ie) spec2s(ie)=spec2s(ie)+spec(ie)*spec(ie) spec(ie)=0.D0 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,fmt="('9')") ! Set end of batch for CG View tt=etime(tarray) tt1=tarray(1) cputime=tt1-tt0 write(6,300) cputime 300 format(' Elapsed Time (sec)=',G15.5) !----------------------------------------------------------------------- ! Step 9: Output-of-results !----------------------------------------------------------------------- write(6,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(6,320) totke,availke,ncount 320 FORMAT(//,' Stopped in MAIN with TotKE=',G15.5,/, * ' AvailKE=',G15.5, /,' Ncount=',I10) stop end if tdet=10.0 rdet=5.0 tcase=0.2 rcase=0.2 write(6,330) tdet,rdet,tcase,rcase 330 FORMAT(/' Detector length=',G15.5,' cm'/ * ' Detector radius=',G15.5,' cm'/ * ' Polyetylene case thickness=',G10.2,' cm'/ * ' Polyetylene case side thickness=',G10.2,' cm'/) write(6,340) ekein 340 FORMAT(' Results for ',G15.5,'MeV photon'/) write(6,350) (media(i,2),i=1,24) 350 FORMAT(' Cover material is ',24A1/) write(6,360) (media(i,1),i=1,24) 360 FORMAT(' Inside chamber is ',24A1/) write(6,365) zin 365 FORMAT(' Cs-137 exists at z=',F8.2,' cm, x=y=0.0'/) ! ----------------------------------- ! Calculate average and its deviation ! ----------------------------------- ! Absorbed energy (MeV) devided by weight of air inside chamber (dweight) ! MeV/g ! 1MeV=1.602E-13J, 1kg=1000g ! 1MeV/g=1.602E-13(J/MeV)*1000(g/kg)=1.602E-10 Gy tdet=10.0 rdet=5.0 vol=pi*rdet*rdet*tdet ! volume of air cm3 weight=vol*rhom(med(1)) ! weight of air inside chamber g ! rhom(1) is a density of air (medium 1) write(6,370) vol,weight 370 FORMAT(' Volume of chamber : ',G15.5,' cm3'/ * ' weight of chamber : ',G15.5,' g') ! --------------- ! Absorbed dose ! --------------- write(6,380) 380 FORMAT(/' Result from absorbed energy in air inside chamber') avdose = doses/ncount dose2s=dose2s/ncount sigdose=dsqrt((dose2s-avdose*avdose)/ncount) avdose = avdose/weight*1.602E-10 sigdose = sigdose/weight*1.602E-10 write(6,390) avdose,sigdose 390 FORMAT(' Air absorbed dose =',G15.4,'+-',G15.4, * ' Gy/source photon') write(6,400) 400 FORMAT(/' Result from absorbed energy in air'/ * 'using track length of photon inside chamber and', * ' mass energy absorption coefficient of air') ! Track lenght (cm) gives photon fluence deviding by volume ! 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 avdose = doseps/ncount dosep2s=dosep2s/ncount sigdose=dsqrt((dosep2s-avdose*avdose)/ncount) avdose = avdose/vol*1.602E-10 sigdose = sigdose/vol*1.602E-10 write(6,410) avdose,sigdose 410 FORMAT(' Air absorbed dose =',G15.4,'+-',G15.4, * ' Gy/source photon') ! ---------------------------------------- ! Average photon spectrum inside chamber. ! ---------------------------------------- write(6,420) 420 FORMAT(/' Average photon spectrum inside chamber'/) do ie=1,50 elow=deltae*(ie-1) eup=deltae*ie avspe = specs(ie)/ncount spec2s(ie)=spec2s(ie)/ncount sigspe=dsqrt((spec2s(ie)-avspe*avspe)/ncount) avspe=avspe/vol/deltae sigspe= sigspe/vol/deltae write(6,430) eup,avspe,sigspe 430 FORMAT(' Upper energy', G15.5,' MeV--',G12.5,'+-',G12.5, * ' photons/cm2/MeV per source photon') 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 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 * dose,dosep,deltae,spec(50),maxpict real*8 dose,dosep,deltae,spec integer maxpict integer ! Arguments * iarg real*8 ! Local variables * edepwt,dcon,decon,encoea 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 .ge. 5) return ! ---------------------------------------------- ! Score energy deposition inside NaI detector ! ---------------------------------------------- if (irl .eq. 1) then ! inside ion chamber dose = dose + edepwt ! -------------------------------------------------------- ! Score particle information if it proceeds inside chamber ! -------------------------------------------------------- if (iql .eq. 0) then ! photon dcon=encoea(e(np)) dosep=dosep+wt(np)*e(np)*dcon*ustep ie = e(np)/deltae +1 if(ie .gt. 50) ie = 50 spec(ie) = spec(ie) + wt(np)*e(np)*dcon*ustep end if 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,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------------------------ !-------------------------------encoea.f-------------------------------- ! Version: 030831-1300 ! Reference: SLAC-265 (p.19-20, Appendix 2) !----------------------------------------------------------------------- !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 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 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=alog(enmu(i+1)) enm0=alog(enmu(i)) hnu1=alog(hnu(i+1)) hnu0=alog(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------------------------ !-------------------------------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---------------------