!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 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 22 Oct 2015 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 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 integer mxbin,ne parameter (mxbin = 1024) common/totals/depe,deltae,spec(3,mxbin),maxpict,ne real*8 depe,deltae,spec,specs(3,mxbin),spec2s(3,mxbin) real*8 totke,rnnow,etot,esumt,wtin real*8 elow,eup,avsp(3),sigsp(3),rdet,tdet,dvol,emaxhisto real*8 wtsum,pefs,pef2s,tefs,tef2s,adoses,adose2s real*8 phs(mxbin),ph2s(mxbin),avpe,avph,avte,avdose,dweight real*8 sigpe,sigph,sigte,sigdose 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 open(90,FILE='spc.dat' ,STATUS='unknown') ! spc-gamma output open(91,FILE='de.dat' ,STATUS='unknown') ! dE output (***day-4) ! (2) Initialization call counters_out(0) call block_set !ooo (3) User setting parameters-1 ncases = 10000 ! number of calculations maxpict = 150 ! number of incident radiations in CG chard(1) = 1d0 ! set character dimensions chard(1) = 10d0 ! set character dimensions !ooo (4) User setting parameters-2 (link PEGS5-material to EGS5 med) nmed=2 medarr(1)='AL ' medarr(2)='AIR ' !ooo (5) Initialize tally parameters ne = 50 ! number of energy-steps for output histogram (max mxbin) emaxhisto = 1.5 ! max energy of the output histogram (MeV) deltae = emaxhisto /real(ne) rdet = 3.81 ! calculate "dvol" (Volume of the tally) (cm^3) tdet = 7.62 dvol = pi*rdet*rdet*tdet !ooo (6) Initialize Track-length tally do j=1,ne do ntype=1,3 spec (ntype,j) = 0.D0 specs (ntype,j) = 0.D0 spec2s(ntype,j) = 0.D0 end do end do !ooo (7) Initialize Energy-deposit tally ncount = 0 wtsum = 0.D0 depe = 0.D0 pefs = 0.D0 pef2s = 0.D0 tefs = 0.D0 tef2s = 0.D0 adoses = 0.d0 adose2s= 0.d0 do j=1,ne phs(j) = 0.D0 ph2s(j) = 0.D0 end do ! (8) Some treatments call geneout(1) call media0(medarr) ! set media call geneout(2) call pegs5 npreci = 3 call cgcontrol(1,4,39) !ooo (9) 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 ! (10) 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) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~ (11) 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 !ooo (12) Track-length tally do ntype=1,3 do ie=1,ne 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 !ooo (13) Energy-deposition tally if(depe .gt. 0.D0) then adoses = adoses + depe adose2s = adose2s + depe*depe ie = depe/deltae + 1 if (ie .gt. ne) ie = ne 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 ! peak definition pefs = pefs + wtin pef2s = pef2s + wtin*wtin end if depe = 0.D0 end if ncount = ncount + 1 if(iwatch.gt.0) call swatch(-1,iwatch) end do !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ MAIN LOOP END !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! (14) 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) !ooo (15) output tallies !ooo track-length write(90,*)"'Energy spectrum (track length)'" write(90,*)'x:E (MeV)' write(90,*)'y:Fluence per source (/cm^2/MeV)' write(90,*)'p:ylog nofr noms' write(90,*)'h:x n y1(photon),lh0b d1 y2(electron),lh0r d2 $ y3(positron),lh0rrr d3' do ie = 1,ne elow = deltae*(ie-1) eup = deltae*ie avsp(1) = specs(1,ie)/real(ncount) ! gamma spectrum spec2s(1,ie)= spec2s(1,ie)/real(ncount) sigsp(1) = dsqrt((spec2s(1,ie)-avsp(1)*avsp(1))/real(ncount)) avsp(1) = avsp(1) /dvol/deltae sigsp(1) = sigsp(1)/dvol/deltae avsp(2) = specs(2,ie) /real(ncount) ! e- spectrum spec2s(2,ie)= spec2s(2,ie)/real(ncount) sigsp(2) = dsqrt((spec2s(2,ie)-avsp(2)*avsp(2))/real(ncount)) avsp(2) = avsp(2) /dvol/deltae sigsp(2) = sigsp(2)/dvol/deltae avsp(3) = specs(3,ie) /real(ncount) ! e+ spectrum spec2s(3,ie)= spec2s(3,ie)/real(ncount) sigsp(3) = dsqrt((spec2s(3,ie)-avsp(3)*avsp(3))/real(ncount)) avsp(3) = avsp(3) /dvol/deltae sigsp(3) = sigsp(3)/dvol/deltae c$$$ elow = deltae*(ie-1) c$$$ eup = deltae*ie c$$$ avspg = specs(1,ie)/real(ncount) ! gamma spectrum c$$$ spec2s(1,ie) = spec2s(1,ie)/real(ncount) c$$$ sigspg = dsqrt((spec2s(1,ie)-avspg*avspg)/real(ncount)) c$$$ avspg = avspg /dvol/deltae c$$$ sigspg = sigspg/dvol/deltae c$$$ avspe = specs(2,ie) /real(ncount) ! e- spectrum c$$$ spec2s(2,ie) = spec2s(2,ie)/real(ncount) c$$$ sigspe = dsqrt((spec2s(2,ie)-avspe*avspe)/real(ncount)) c$$$ avspe = avspe /dvol/deltae c$$$ sigspe = sigspe/dvol/deltae c$$$ avspp = specs(3,ie) /real(ncount) ! e+ spectrum c$$$ spec2s(3,ie) = spec2s(3,ie)/real(ncount) c$$$ sigspp = dsqrt((spec2s(3,ie)-avspp*avspp)/real(ncount)) c$$$ avspp = avspp /dvol/deltae c$$$ sigspp = sigspp/dvol/deltae write(90,'(2f10.3,6(1pe12.5))') $ elow,eup,avsp(1),sigsp(1),avsp(2),sigsp(2),avsp(3),sigsp(3) end do !ooo Peak efficiency avpe = pefs /real(ncount) pef2s = pef2s/real(ncount) sigpe = dsqrt((pef2s-avpe*avpe)/real(ncount)) avpe = avpe *100.0 sigpe = sigpe*100.0 write(6,*)'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' write(6,'(A,G11.4,A,G9.2,A)') $ ' Peak efficiency =',avpe,'+-',sigpe,' %' !ooo Total efficiency avte = tefs /real(ncount) tef2s = tef2s/real(ncount) sigte = dsqrt((tef2s-avte*avte)/real(ncount)) avte = avte *100.0 sigte = sigte*100.0 write(6,'(A,G11.4,A,G9.2,A)') * ' Total efficiency =', avte,'+-',sigte,' %' !ooo 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 dweight = dvol*rhom(1) avdose = adoses /real(ncount) adose2s = adose2s/real(ncount) sigdose = dsqrt( (adose2s-avdose*avdose)/real(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' write(6,*)'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' !ooo enegy deposition tally write(91,*)"'Energy deposition'" write(91,*)'x:E (MeV)' write(91,*)'y:Events' write(91,*)'p:ylog nofr noms' write(91,*)'h:x n y,lh0cc d' do ie = 1,ne elow = deltae*(ie-1) eup = deltae*ie avph = phs(ie) / real(ncount) ph2s(ie) = ph2s(ie) / real(ncount) sigph = dsqrt((ph2s(ie)-avph*avph)/real(ncount)) ! avph = avph/deltae ! in case of (/MeV) output ! sigph = sigph/deltae ! in case of (/MeV) output write(91,'(2f10.3,g12.5,g12.5)')elow,eup,avph,sigph end do write(6,*)'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' write(6,*)' EGS5 sucessfully finished! ' write(6,*)'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' end !ooo (16) Source description subroutine source(wtin) implicit none include 'auxcommons/instuf.f' include 'include/randomm.f' real*8 pi,wtin,phi,thet,r4,r5 real*8 r1,r2,r3,a0,b0,r0,z0 parameter (pi=3.14159265358979d0) iqin = 0 ! incident particle charge ekein = 1.0 ! incident particle kinetic energy ! cylinder source r0 = 2.5 ! radius of the source cylinder (cm) a0 = 5 ! height of the source cylinder (cm) z0 = -10 ! source position (distance to detector) b0 = max(r0,a0) ! pick up longer of radius or height 1 call randomset(r1) call randomset(r2) call randomset(r3) call randomset(r4) call randomset(r5) xin = (2*r1 - 1)*b0 yin = (2*r2 - 1)*b0 zin = (2*r3 - 1)*b0 + z0 if( zin .lt. (z0-a0/2) .or. zin .gt. (z0+a0/2)) goto 1 if( ( xin**2.0 + yin**2.0) .gt. r0**2.0 ) goto 1 phi = r4 *pi ! phi : 0 ~ pi thet = (2*r5 - 1)*pi ! theta : -pi ~ pi uin = sin(phi)*cos(thet) vin = sin(phi)*sin(thet) win = cos(phi) wtin = 1.0 ! weight if(abs(uin*uin+vin*vin+win*win-1.0).gt.1.e-6) then write(6,'(A,3E12.5)')'The src-direction are not normalized.', * uin,vin,win stop end if 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' integer mxbin,ne parameter (mxbin = 1024) common/totals/depe,deltae,spec(3,mxbin),maxpict,ne real*8 depe,deltae,spec,edepwt integer iarg,maxpict,ie,iql,irl,ntype irl = ir(np) iql = iq(np) !ooo set deposit energy edepwt = edep*wt(np) 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 if (iwatch .gt. 0) call swatch(iarg,iwatch) if(iarg .ge. 5) return !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !ooo Tally (for single-irl ver.) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if ( irl.eq.1 ) then !ooo Energy deposition tally depe = depe + edepwt !ooo Track length tally if ( iql.eq.0 ) then ! photon ntype = 1 ie = e(np)/deltae +1 if( ie.gt.mxbin ) ie = mxbin elseif ( iql.eq.-1 ) then ! electron ntype = 2 ie = (e(np) - RM)/deltae +1 if( ie.gt.mxbin ) ie = mxbin else ! positron ntype = 3 ie = (e(np) - RM)/deltae +1 if( ie.gt.mxbin) ie = mxbin endif spec(ntype,ie) = spec(ntype,ie) + ustep*wt(np) endif !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Tally end !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 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