!*********************************************************************** !***************************** KEK, High Energy Accelerator Research * !** u c n a i c g v _ base ** Organization * !***************************** EGS5.0 USER CODE - 09 Oct 2015/1430 * !*********************************************************************** !* This is a general User Code based on the cg geometry scheme. * !*********************************************************************** ! Original "ucnaicgv_vsp.f" by H. Hirayama and Y. Naito ! the simplified version by H. Iwase !*********************************************************************** implicit none include 'include/egs5_h.f' 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' include 'auxcommons/aux_h.f' include 'auxcommons/edata.f' include 'auxcommons/etaly1.f' include 'auxcommons/instuf.f' include 'auxcommons/lines.f' include 'auxcommons/nfac.f' include 'auxcommons/watch.f' include 'auxcommons/geom_common.f' ! geom-common file common/totals/depe,deltae,spec(3,50),maxpict real*8 depe,deltae,spec real*8 totke,rnnow,etot,esumt,wtin integer maxpict integer i,icases,idin,ie,ifti,ifto,ii,iiz,imed,ireg,isam, * izn,nlist,j,k,n,ner,ntype,idum character*24 medarr(MXMED) ! (1) Open files open(6, FILE='egs5job.out',STATUS='unknown') ! General output open(4, FILE='egs5job.inp',STATUS='old') ! CG input open(39,FILE='egs5job.pic',STATUS='unknown') ! CG output ! (2) Initialization call counters_out(0) call block_set !#####(3) User setting parameters-1 (link PEGS5-material to EGS5 med) nmed=1 medarr(1)='AL ' !#####(4) User setting parameters-2 ncases = 1000 ! number of calculations maxpict = 200 ! number of incident radiations in CG chard(1) = 7.62d0 ! set character dimensions chard(2) = 0.1d0 chard(3) = 0.5d0 chard(4) = 5.0d0 ! (5) Some treatments call geneout(1) call media0(medarr) ! set media call geneout(2) call pegs5 npreci = 3 call cgcontrol(1,4,39) !#####(6) User setting parameters-3 (EGS options) 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 ! (7) Some treatments call rluxinit0 emaxe = 0.D0 ! Autoset ecuts when emaxe=0 call hatch0 call geneout(3) call geneout(4) call cgcontrol(4,idum,idum) call cgcontrol(5,idum,idum) if(iwatch.gt.0) call swatch(-99,iwatch) !~~~~ (8) MAIN LOOP ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ MAIN LOOP do i=1,ncases call source(wtin) ! source generating here irin = 0 ! set 0 for Autoset call cgcontrol(2,idum,idum) ! Autoset irin etot = ekein + iabs(iqin)*RM call shower (iqin,etot,xin,yin,zin,uin,vin,win,irin,wtin) ! main calc ncount = ncount + 1 if(iwatch.gt.0) call swatch(-1,iwatch) end do !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ MAIN LOOP END ! (9) Some treatments after calculation 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) call cgcontrol(6,idum,idum) call counters_out(1) write(6,*)'----- EGS5 sucessfully finished! -----' end !#####(10) Source description subroutine source(wtin) implicit none include 'auxcommons/instuf.f' include 'include/randomm.f' real*8 wtin,phi,thet,r1,r2,pi parameter (pi=3.14159265358979d0) iqin = 0 ! incident particle charge ekein = 1.253 ! incident particle kinetic energy xin = 0.0 ! source position of x yin = 0.0 ! source position of ya zin = -7.0 ! source position of z call randomset(r1) call randomset(r2) phi = r1*pi thet = (2*r2 - 1)*pi uin = sin(phi)*cos(thet) vin = sin(phi)*sin(thet) win = cos(phi) wtin = 1.0 ! weight end !-------------------------------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 irl = ir(np) iql = iq(np) if (iwatch .gt. 0) call swatch(iarg,iwatch) if(iarg .ge. 5) return 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 !----------------------------------------------------------------------- ! Required (geometry) subroutine for use with the EGS5 Code System ! ---------------------------------------------------------------------- ! This is a CG-HOWFAR. ! ---------------------------------------------------------------------- subroutine howfar implicit none 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 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 ir_np = ir(np) iq_np = iq(np) + 2 if(ir_np.le.0) then write(6,*) 'Stopped in howfar with ir(np) <=0' stop end if if(ir_np.gt.izonin) then write(6,*) 'Stopped in howfar with ir(np) > izonin' stop end if if(ir_np.EQ.izonin) then idisc=1 return end if tval=1.d+30 itvalm=0 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) 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--------------------- subroutine geneout(iout) implicit none !------------------------------------------------------------------ ! iout and general outputs ! 1 Media information (general output) ! 2 Character dimensions (general output) ! 3 Media information (general output) ! 4 X-ray information (general output) !------------------------------------------------------------------ integer iout,i,j,ii,iiz,izn,ner include 'include/egs5_h.f' 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' include 'auxcommons/aux_h.f' include 'auxcommons/edata.f' include 'auxcommons/etaly1.f' include 'auxcommons/instuf.f' include 'auxcommons/lines.f' include 'auxcommons/nfac.f' include 'auxcommons/watch.f' include 'auxcommons/geom_common.f' ! geom-common file if( iout.eq.1 )then 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 elseif( iout.eq.2)then write(6,fmt="('chard =',5e12.5)") (chard(j),j=1,nmed) write(6,'(A/)') 'PEGS5-call comes next' elseif( iout.eq.3 )then 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 elseif( iout.eq.4 )then 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' 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 !---------------- endif !---------------- end subroutine cgcontrol(icnt,i1,i2) implicit none !------------------------------------------------------------------ ! icnt and treatment ! 1 call geomgt and write dummy geom ! 2 ! 3 ! 4 region and medium information (CG output) !------------------------------------------------------------------ integer i,j,icnt,ifti,ifto,i1,i2 include 'include/egs5_h.f' 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' include 'auxcommons/aux_h.f' include 'auxcommons/edata.f' include 'auxcommons/etaly1.f' include 'auxcommons/instuf.f' include 'auxcommons/lines.f' include 'auxcommons/nfac.f' include 'auxcommons/watch.f' include 'auxcommons/geom_common.f' ! geom-common file if( icnt.eq.1 )then ifti = i1 ifto = i2 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' nreg = izonin read(4,*) (med(i),i=1,nreg) elseif( icnt.eq.2 )then 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 elseif( icnt.eq.4 )then write(39,fmt="('MSTA')") write(39,fmt="(i4)") nreg write(39,fmt="(15i4)") (med(i),i=1,nreg) write(39,fmt="('MEND')") elseif( icnt.eq.5 )then write(39,fmt="('0 1')") elseif( icnt.eq.6 )then write(39,'(A)') '9' ! Set end of batch for CG View endif end subroutine media0(medarr) implicit none include 'include/egs5_h.f' include 'include/egs5_media.f' character*24 medarr(MXMED) integer i,j do j=1,nmed do i=1,24 media(i,j)=medarr(j)(i:i) end do end do end subroutine rluxinit0 include 'include/randomm.f' 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 end subroutine hatch0 include 'include/egs5_h.f' include 'include/egs5_misc.f' write(6,'(/A)') ' Call hatch to get cross-section data' open(UNIT=KMPI,FILE='pgs5job.pegs5dat',STATUS='old') open(UNIT=KMPO,FILE='egs5job.dummy',STATUS='unknown') call hatch close(UNIT=KMPI) close(UNIT=KMPO) end