subroutine scatb(beamflg, en0, x0, y0, z0, & thx0, thy0, xin1,xin2,yin1,yin2, zlm, & xsc, ysc, zsc, scattFlag, enp, x, y, z, & thx, thy, xstot, thsc, lam, rfail) implicit none include 'const.inc' integer hws, scattFlag, beamflg real vec(4), HRNDM1 double precision en0, k0, x0, y0, z0, thx0, thy0 double precision enp, kp, x, y, z, thx, thy double precision xsc, ysc, zsc double precision xstot, norm, hedens, denat, dreact double precision q, ep, wp, th_test, thsc, phsc double precision tharg, pharg, thp, php, sinphp double precision th0, ph0, thxp, thyp double precision xin1, xin2, yin1, yin2, pla, plb double precision rfx, rfy, rfz, zx, zy, xlm, ylm, zlm double precision xdis, ydis, xdis0, ydis0, rdist(3) double precision xt, yt double precision thscCM, rflx, rfly, lam double precision, external :: rsqrt,ratan logical rfail rfail = .false. ccccc reflectivity if (reflectflg) then rflx=1.0 ; rfly=1.0 else if (dabs(thx0)/lam.lt.0.00173) then rflx = 1.0 c rflx = 1.0-17.34*dabs(thx0)/lam else if ((beamflg.eq.1).or.(beamflg.eq.2)) then if((dabs(thx0)/lam.ge.0.00173).and. & (dabs(thx0)/lam.le.ginx)) then c rflx = 1-(0.20/(ginx-0.00173))*(dabs(thx0)/lam-0.00173) rflx = 1-75.38*(dabs(thx0)/lam-0.00173) if (rflx.lt.0.0d0) rflx = 0.0d0 else if (dabs(thx)/lam.ge.ginx) then rflx = 0.d0 end if else if (beamflg.eq.3) then if((dabs(thx0)/lam.ge.0.00173).and. & (dabs(thx0)/lam.le.gotx)) then c rflx = 1-(0.20/(ginx-0.00173))*(dabs(thx0)/lam-0.00173) rflx = 1-75.38*(dabs(thx0)/lam-0.00173) if (rflx.lt.0.0d0) rflx = 0.0d0 else if (dabs(thx)/lam.ge.gotx) then rflx = 0.d0 end if end if cccc y if (dabs(thy0)/lam.lt.0.00173) then rfly = 1.0 c rfly = 1.0-17.34*dabs(thy0)/lam else if ((beamflg.eq.1).or.(beamflg.eq.2)) then if((dabs(thy0)/lam.ge.0.00173).and. & (dabs(thy0)/lam.le.giny)) then c rfly = 1-(0.20/(giny-0.00173))*(dabs(thy0)/lam-0.00173) rfly = 1-75.38*(dabs(thy0)/lam-0.00173) if (rfly.lt.0.0d0) rfly = 0.0d0 else if (dabs(thy)/lam.ge.giny) then rfly = 0.d0 end if else if (beamflg.eq.3) then if((dabs(thy0)/lam.ge.0.00173).and. & (dabs(thy0)/lam.le.goty)) then c rfly = 1-(0.20/(giny-0.00173))*(dabs(thy0)/lam-0.00173) rfly = 1-75.38*(dabs(thy0)/lam-0.00173) if (rfly.lt.0.0d0) rfly = 0.0d0 else if (dabs(thy)/lam.ge.goty) then rfly = 0.d0 end if end if ccccc k0 = rsqrt(2*m_n*en0)/h_par th0 = ratan(rsqrt(dtan(thx0)**2 + dtan(thy0)**2)) if(thx0.ge.0..and.thy0.ge.0.) then ph0 = ratan(dtan(thy0)/dtan(thx0)) elseif(thx0.lt.0..and.thy0.ge.0.) then ph0 = pi + ratan(dtan(thy0)/dtan(thx0)) elseif(thx0.lt.0..and.thy0.lt.0.) then ph0 = pi + ratan(dtan(thy0)/dtan(thx0)) elseif(thx0.ge.0..and.thy0.lt.0.) then ph0 = 2.*pi + ratan(dtan(thy0)/dtan(thx0)) endif if (hebeam) then ! He gas in beam c Get sigma(en0) en0 in (meV) and find reaction length c !!This sigma(E) for air gap is fit of intepreted data from ENDF -- c for n-4He but not much actual data!!!!! c if (en0.lt.6.0) then xstot = 0.0899/(en0/1000.0)**(0.474) else xstot = 0.460/(en0/1000.0)**(0.173) endif denat = 0.0000269 !room temperature helium norm = 1.0 else !air in beam c Get sigma(en0) en0 in (meV) and find reaction length c !!This sigma(E) for air gap is fit of intepreted data from ENDF -- c for n-14N , 14N cross section much bigger than 16O c but not much actual data!!!!! if (en0.lt.2.2) then xstot = 1.0742/(en0/1000.0)**(0.47) else xstot = 6.7799/(en0/1000.0)**(0.162) endif denat = 0.000054 !room temperature norm = 1.0 end if call ranlux(vec,4) dreact = -1.0*log(vec(1))/(denat*xstot*norm) c distance at which scattering will happen from the beginning of beam section rfx = dreact*dsin(th0)*dcos(ph0) rfy = dreact*dsin(th0)*dsin(ph0) rfz = dreact*dcos(th0) xsc = x0 + rfx ! entering position + dist. to scatter ysc = y0 + rfy zsc = z0 + rfz ccccccccccc Scatterings ccccccccccc if( (xsc.ge.xin1.and.xsc.le.xin2) .and. & (ysc.ge.yin1.and.ysc.le.yin2) .and. & (zsc.ge.z0.and.zsc.le.zlm) ) then ! scat inside beampipe scattFlag = 1 c*** Want isotropic scattering angle for air and He scat thscCM = dacos(1.d0-2.d0*vec(3)) ! isotropic Thsc in CM frame thsc = ratan(dsin(thscCM)/(dcos(thscCM)+(m_n/m_he))) if (thsc.lt.0.) thsc = pi + thsc ep = en0*(m_n/(m_n+m_he))**2*( dcos(thsc)+ ! Lab Thsc, but with z along k & rsqrt( (m_he/m_n)**2-(dsin(thsc))**2 ) )**2 if (ep.le.0.01d0) then write(*,*) ' ep < 0.01:',ep goto 1000 end if phsc = 2.0*pi*vec(2) c*** transform from z along k to z along beam transport axis tharg = -dsin(th0)*dsin(thsc)*dcos(phsc)+dcos(th0)*dcos(thsc) if (tharg.gt.1.0d0) then write(*,*) ' theta_p error : cos(th_p) =',tharg tharg=1.0d0 else if (tharg.lt.-1.0d0) then write(*,*) ' theta_p error : cos(th_p) =',tharg tharg=-1.0d0 endif thp = dacos(tharg) if (thp.lt.1.0d-9) then pharg = 1.0 else pharg = (dsin(thsc)*dcos(phsc)*dcos(ph0)*dcos(th0) - & dsin(thsc)*dsin(phsc)*dsin(ph0) + & dcos(thsc)*dcos(ph0)*dsin(th0) )/dsin(thp) end if if (pharg.gt.1.0d0) then write(*,*) ' phi_p error : cos(ph_p) =',pharg pharg=1.0d0 else if (pharg.lt.-1.0d0) then write(*,*) ' phi_p error : cos(ph_p) =',pharg pharg=-1.0d0 end if php = dacos(pharg) sinphp = (dsin(thsc)*dcos(phsc)*dsin(ph0)*dcos(th0) + & dsin(thsc)*dsin(phsc)*dcos(ph0) + & dcos(thsc)*dsin(ph0)*dsin(th0) )/dsin(thp) if (sinphp.lt.0) php = 2.0d0*pi - php thxp = ratan(dtan(thp)*dcos(php)) thyp = ratan(dtan(thp)*dsin(php)) xdis = xsc + (zlm-z0-rfz)*dtan(thxp) ! x loc at end of pipe (scat in pipe) ydis = ysc + (zlm-z0-rfz)*dtan(thyp) ! y loc at end of pipe if( (xdis.ge.xin1.and.xdis.le.xin2) .and. & (ydis.ge.yin1.and.ydis.le.yin2) ) then ! at zlm inside xlm,ylm x = xdis y = ydis z = zlm thx = thxp thy = thyp enp = ep else if(xdis.ge.xin1.and.xdis.le.xin2) then ! at zlm, inside xlm if(ydis.lt.yin1) then ! and outside ylm- y = yin1 ! return to target with current position at wall z = (yin1-ysc)/dtan(thyp) + zsc x = xsc + (yin1-ysc)*dtan(thxp)/dtan(thyp) else if(ydis.gt.yin2) then ! and outside ylm+ y = yin2 z = (yin2-ysc)/dtan(thyp) + zsc x = xsc + (yin2-ysc)*dtan(thxp)/dtan(thyp) endif thx = thxp ! so reflect in y thy = -1*thyp enp = ep if (vec(4).gt.rfly) goto 1000 ! apply Ry to neutron else if(ydis.ge.yin1.and.ydis.le.yin2) then ! at zlm, inside ylm if(xdis.lt.xin1) then ! and outside xlm- x = xin1 z = (xin1-xsc)/dtan(thxp) + zsc y = ysc + (xin1-xsc)*dtan(thyp)/dtan(thxp) else if(xdis.gt.xin2) then ! and outside xlm- x = xin2 z = (xin2-xsc)/dtan(thxp) + zsc y = ysc + (xin2-xsc)*dtan(thyp)/dtan(thxp) endif thx = -1*thxp ! so reflect in x thy = thyp enp = ep if (vec(4).gt.rflx) goto 1000 ! apply Rx to neutron else ! at zlm outside pipe x-y if(xdis.lt.xin1) then xlm = xin1 ! set limit to pipe wall else if(xdis.gt.xin2) then xlm = xin2 endif if(ydis.lt.yin1) then ylm = yin1 else if(ydis.gt.yin2) then ylm = yin2 endif zx = (xlm-xsc)/dtan(thxp) ! delta z btw wall hit and scat zy = (ylm-ysc)/dtan(thyp) if(zx.le.zy) then x = xlm z = zx + zsc y = ysc + zx*dtan(thyp) thx = -1*thxp ! reflect in x thy = thyp if (vec(4).gt.rflx) goto 1000 ! apply Rx to neutron else y = ylm z = zy + zsc x = xsc + zy*dtan(thxp) thx = thxp thy = -1*thyp ! reflect in y if (vec(4).gt.rfly) goto 1000 ! apply Ry to neutron endif enp = ep endif else ! scatter outside pipe, i.e., scat reaction length is scattFlag = 0 ! such that xcs, ysc, or zsc (scatter point) is outside pipe xsc = x0 ! so xlm, ylm are defined as reflecting wall ysc = y0 ! and the angle of reflecting direction is inverted zsc = z0 ! xt = (zlm-z0)*dtan(thx0)+x0 ! x-position end of target yt = (zlm-z0)*dtan(thy0)+y0 if( (xt.ge.xin1.and.xt.le.xin2) .and. & (yt.ge.yin1.and.yt.le.yin2) ) then ! inside x-y walls at end of pipe x = xt y = yt z = zlm thx = thx0 thy = thy0 enp = en0 else if(xt.ge.xin1.and.xt.le.xin2) then ! inside x-walls at end of pipe if(yt.lt.yin1) then ! y < neg y-wall y = yin1 z = (yin1-y0)/dtan(thy0) + z0 ! position at wall hit x = x0 + (yin1-y0)*dtan(thx0)/dtan(thy0) ! x0+z*tan(thx) else if(yt.gt.yin2) then ! y > pos y-wall y = yin2 z = (yin2-y0)/dtan(thy0) + z0 x = x0 + (yin2-y0)*dtan(thx0)/dtan(thy0) endif thx = thx0 thy = -1*thy0 ! reflect in y enp = en0 if (vec(4).gt.rfly) goto 1000 ! apply Ry to neutron else if(yt.ge.yin1.and.yt.le.yin2) then ! inside y-walls at end of pipe if(xt.lt.xin1) then x = xin1 z = (xin1-x0)/dtan(thx0) + z0 y = y0 + (xin1-x0)*dtan(thy0)/dtan(thx0) else if(xt.gt.xin2) then x = xin2 z = (xin2-x0)/dtan(thx0) + z0 y = y0 + (xin2-x0)*dtan(thy0)/dtan(thx0) endif thx = -1*thx0 ! reflect in x thy = thy0 enp = en0 if (vec(4).gt.rflx) goto 1000 ! apply Rx to neutron else ! outside x-y walls at end of pipe if(xt.lt.xin1) then xlm = xin1 else if(xt.gt.xin2) then xlm = xin2 endif if(yt.lt.yin1) then ylm = yin1 else if(yt.gt.yin2) then ylm = yin2 endif zx = (xlm-x0)/dtan(thx0) ! delta z btw zo and z at wall hit zy = (ylm-y0)/dtan(thy0) if(zx.le.zy) then ! position at first x or y wall hit x = xlm z = zx + z0 y = y0 + zx*dtan(thy0) thx = -1*thx0 ! reflect in x thy = thy0 if (vec(4).gt.rflx) goto 1000 ! apply Rx to neutron else y = ylm z = zy + z0 x = x0 + zy*dtan(thx0) thx = thx0 thy = -1*thy0 ! reflect in y if (vec(4).gt.rfly) goto 1000 ! apply Ry to neutron endif enp = en0 endif endif ccccccccccc goto 1001 1000 rfail = .true. ! neutron is lost through wall 1001 return end