!*********************************************************************** !***************************** KEK, High Energy Accelerator Research * !***************************** Organization * !*** u c n a i c g v ********* * !***************************** EGS5.0 USER CODE - 28 Jul 2012/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-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 ucnaicgv.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 corresponds to ucnai3cgp.mor for egs4. * ! Use Ranlux random number generator. * !*********************************************************************** ! * ! ----------------------- * ! cg Geometry (ucnaicgv) * ! ----------------------- * ! * ! * ! 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_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,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, * rr0,sigpe,sigte,sigph,sigspg,sigspe,sigspp,tefs,tef2s,wtin,wtsum, * xi0,yi0,zi0 real*8 * phs(50),ph2s(50),specs(3,50),spec2s(3,50) real ! Local variables * elow,eup,rdet,rtcov,rtgap,tcov,tdet,tgap 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,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 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,120) inseed 120 FORMAT(/,' inseed=',I12,5X, * ' (seed for generating unique sequences of Ranlux)') ! ============= call rluxinit ! Initialize the Ranlux random-number generator ! ============= !----------------------------------------------------------------------- ! Step 4: Determination-of-incident-particle-parameters !----------------------------------------------------------------------- ! Define initial variables for incident particle normally incident ! on the slab iqin=0 ! Incident particle charge - photons ekein=1.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,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=',G10.5,' MeV, pcut=',G10.5,' MeV') ! ----------------------------------------------- ! Print out energy information of K- and L-X-rays ! ----------------------------------------------- if (iedgfl(i) .ne. 0) then ! Output X-ray energy ner = nne(med(i)) do iiz=1,ner izn = zelem(med(i),iiz) ! Atomic number of this element write(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 depe=0.D0 pefs=0.D0 pef2s=0.D0 tefs=0.D0 tef2S=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=10000 ! Set maximum number for pict maxpict=50 tt=etime(tarray) tt0=tarray(1) !----------------------------------------------------------------------- ! Step 8: Shower-call !----------------------------------------------------------------------- ! Write batch number write(39,fmt="('0 1')") ! ======================== if(iwatch.gt.0) call swatch(-99,iwatch) ! ======================== ! ------------------------- do i=1,ncases ! Start of shower call-loop ! ------------------------- ! ---------------------- ! 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 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,irin,wtin) ! ========================================================== ! If some energy is deposited inside detector add pulse-height ! and efficiency. if (depe .gt. 0.D0) then 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,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=7.62 rdet=3.81 tcov=0.1 rtcov=0.1 tgap=0.5 rtgap=0.5 write(6,330) tdet,rdet,tcov,rtcov,tgap,rtgap 330 FORMAT(/' Detector length=',G15.5,' cm'/ * ' Detector radius=',G15.5,' cm'/ * ' Al cover thickness=',G10.2,' cm'/ * ' Al cover side thickness=',G10.2,' cm'/ * ' Front gap =',G10.2,' cm'/' Side gap =',G10.2,' cm'/) write(6,340) ekein 340 FORMAT(' Results for ',G15.5,'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,350) avpe,sigpe 350 FORMAT(' Peak efficiency =',G11.4,'+-',G9.2,' %') ! ---------------- ! Total efficiency ! ---------------- avte = tefs/ncount tef2s = tef2s/ncount sigte = dsqrt((tef2s-avte*avte)/ncount) avte = avte*100.0 sigte = sigte*100.0 write(6,360) avte,sigte 360 FORMAT(' Total efficiency =',G11.4,'+-',G9.2,' %') ! -------------------------- ! Pulse height distribution ! -------------------------- write(6,370) 370 FORMAT(/' Pulse height distribution ') do ie=1,50 elow=deltae*(ie-1) eup=deltae*ie avph = phs(ie)/ncount ph2s(ie)=ph2s(ie)/ncount sigph=dsqrt((ph2s(ie)-avph*avph)/ncount) avph = avph/deltae sigph= sigph/deltae write(6,380) eup,avph,sigph 380 FORMAT(' E (upper-edge --',G10.4,' MeV )=',G15.5,'+-',G15.5, * ' counts/MeV/incident') end do ! ---------------------------------------------------------- ! Particle spectrum. Incident particle spectrum to detector. ! ---------------------------------------------------------- write(6,400) 400 FORMAT(/' Particle spectrum crossing the detector plane'/ * 30X,'particles/MeV/source photon'/ * ' Upper energy',11X,' Gamma',18X,' Electron', * 14X,' Positron') do ie=1,50 elow=deltae*(ie-1) eup=deltae*ie ! ---------------------------------- ! 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/deltae sigspg= sigspg/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/deltae sigspe= sigspe/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/deltae sigspp= sigspp/deltae write(6,410) eup,avspg,sigspg,avspe,sigspe,avspp,sigspp 410 FORMAT(G10.5,' MeV--',3(G12.5,'+-',G12.5)) 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 * 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 particle information if it enters from outside ! ---------------------------------------------------- if (irl .ne. irold .and. iarg .eq. 0) then 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) + wt(np) 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------------------------ !-------------------------------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---------------------