subroutine scatt(q, een0, x0, y0, z0, thx0, thy0, & txin1,txin2,tyin1,tyin2,zlm, & xsc, ysc, zsc, scattFlag, & enp, x, y, z, thx, thy,xstot,thsc,lam,rfail,scwgt) implicit none include 'const.inc' integer hws, scattFlag real vec(4), HRNDM1 double precision een0, k0, x0, y0, z0, thx0, thy0 double precision enp, kp, x, y, z, thx, thy double precision xsc, ysc, zsc, scwgt 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 txin1, txin2, tyin1, tyin2, pla, plb double precision rfx, rfy, rfz, zx, zy, xlm, ylm, zlm ! zlm set in target as tarz doesn't change 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 (gtindexx.lt.1.0E-5) then rflx = 0.0 else if (dabs(thx0)/lam.lt.0.00173) then rflx = 1.0 c rflx = 1.0-17.34*dabs(thx0)/lam else if ((dabs(thx0)/lam.ge.0.00173).and. & (dabs(thx0)/lam.le.gtindexx)) then c rflx = 1-(0.20/(gtindexx-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.gtindexx) then rflx = 0.d0 end if cccc y if (gtindexy.lt.1.0E-5) then rfly = 0.0 else if (dabs(thy0)/lam.lt.0.00173) then rfly = 1.0 c rfly = 1.0-17.34*dabs(thy0)/lam else if ((dabs(thy0)/lam.ge.0.00173).and. & (dabs(thy0)/lam.le.gtindexy)) then c rfly = 1-(0.20/(gtindexy-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.gtindexy) then rfly = 0.d0 end if end if ccccc k0 = rsqrt(2*m_n*een0)/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 (hegastar) then ! if using He gas target c if (j.eq.1) write(6,*) ' !! He gas target run !! ' c if (j.eq.1) write(3,*) ' !! He gas target run !! ' c Get sigma(een0) een0 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 (een0.lt.6.0) then xstot = 0.0899/(een0/1000.0)**(0.474) else xstot = 0.460/(een0/1000.0)**(0.173) endif denat = ndenhegas norm = 1.0 else if (airtar) then c Get sigma(een0) een0 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 (een0.lt.2.2) then xstot = 1.0742/(een0/1000.0)**(0.47) else xstot = 6.7799/(een0/1000.0)**(0.162) endif denat = ndenairtar norm = 1.0 else call sommers(T, een0, xstot) call nden(T,norm) call heden(T,hedens) denat = hedens*avog/he4mol c bec affects mltp scat xstot = 5.d0*xstot ! increasing scattering likelihood - factor here and below xstot = xstot ! increasing scattering likelihood - factor here and below end if call ranlux(vec,4) dreact = -1.0*log(vec(1))/(denat*xstot*norm) c bec affects mltp scat if (hegastar) then ! if using He gas target c bec affects mltp scat scwgt=1.0d0 c bec affects mltp scat else if (airtar) then c bec affects mltp scat scwgt=1.0d0 c bec affects mltp scat else cc bec affects mltp scat scwgt = (1.0d0-dexp(-denat*xstot*norm/5.0d0*dreact))/ ! prob. of scat with 5*xstot cc bec affects mltp scat & (1.0d0-dexp(-denat*xstot*norm*dreact)) c bec affects mltp scat dreact = dreact/5.0d0 ! increasing scattering likelihood c bec affects mltp scat scwgt = 1.0d0/5.0d0 c bec affects mltp scat end if c distance at which scattering will happen from the target upstream rfx = dreact*dsin(th0)*dcos(ph0) rfy = dreact*dsin(th0)*dsin(ph0) rfz = dreact*dcos(th0) xsc = x0 + rfx ! entering position + dist. to target scatter ysc = y0 + rfy zsc = z0 + rfz ccccccccccc Scatterings ccccccccccc if( (xsc.ge.txin1.and.xsc.le.txin2) .and. & (ysc.ge.tyin1.and.ysc.le.tyin2) .and. & (zsc.ge.z0.and.zsc.le.zlm) ) then ! scat in target since scatter point in tar scattFlag = 1 c Get w-histogram id associated with q hws = q*1.0d3+700.0d0 call hrget(hws,'hbooks/sqw.hbook','') wp=HRNDM1(hws) call hdelet(hws) k0 = rsqrt(2*m_n*een0)/h_par ep = een0 - (h_par*wp) if (ep.lt.0.01d0) then write(6,*) ' Scatt tar, ep<0.01 ',ep goto 1000 endif kp = rsqrt(2*m_n*ep)/h_par th_test = (k0**2 + kp**2 - q**2)/(2.*k0*kp) if(th_test.gt.1.0) th_test = 1.0 if(th_test.lt.-1.0) th_test = -1.0 ! see below thsc = dacos(th_test) phsc = 2.0*pi*vec(2) c*** Want isotropic scattering angle if NOT liquid He scattering if (hegastar.or.airtar) then ! if using He or Air gas target - isotropic 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 = een0*(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 else thsc = dacos(th_test) ! liquid He, thsc from k, k', q end if if (ep.le.0.01d0) then write(*,*) ' ep < 0.01:',ep goto 1000 end if 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 tar (scat in tar) ydis = ysc + (zlm-z0-rfz)*dtan(thyp) ! y loc at end of tar if( (xdis.ge.txin1.and.xdis.le.txin2) .and. & (ydis.ge.tyin1.and.ydis.le.tyin2) ) then ! at zlm inside xlm,ylm after scattering x = xdis y = ydis z = zlm thx = thxp thy = thyp enp = ep else if(xdis.ge.txin1.and.xdis.le.txin2) then ! at zlm, inside xlm if(ydis.lt.tyin1) then ! and outside ylm- y = tyin1 ! return to target with current position at wall z = (tyin1-ysc)/dtan(thyp) + zsc x = xsc + (tyin1-ysc)*dtan(thxp)/dtan(thyp) else if(ydis.gt.tyin2) then ! and outside ylm+ y = tyin2 z = (tyin2-ysc)/dtan(thyp) + zsc x = xsc + (tyin2-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.tyin1.and.ydis.le.tyin2) then ! at zlm, inside ylm if(xdis.lt.txin1) then ! and outside xlm- x = txin1 z = (txin1-xsc)/dtan(thxp) + zsc y = ysc + (txin1-xsc)*dtan(thyp)/dtan(thxp) else if(xdis.gt.txin2) then ! and outside xlm- x = txin2 z = (txin2-xsc)/dtan(thxp) + zsc y = ysc + (txin2-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 target if(xdis.lt.txin1) then xlm = txin1 ! set limit to tar wall else if(xdis.gt.txin2) then xlm = txin2 endif if(ydis.lt.tyin1) then ylm = tyin1 else if(ydis.gt.tyin2) then ylm = tyin2 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 ! scater outside target, i.e., scat reaction length is scattFlag = 0 ! such that xcs, ysc, or zsc (scatter point) is outside target 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.txin1.and.xt.le.txin2) .and. & (yt.ge.tyin1.and.yt.le.tyin2) ) then ! inside x-y walls at end of tar x = xt y = yt z = zlm thx = thx0 thy = thy0 enp = een0 else if(xt.ge.txin1.and.xt.le.txin2) then ! inside x-walls at end of tar if(yt.lt.tyin1) then ! y < neg y-wall y = tyin1 z = (tyin1-y0)/dtan(thy0) + z0 ! position at wall hit x = x0 + (tyin1-y0)*dtan(thx0)/dtan(thy0) ! x0+z*tan(thx) else if(yt.gt.tyin2) then ! y > pos y-wall y = tyin2 z = (tyin2-y0)/dtan(thy0) + z0 x = x0 + (tyin2-y0)*dtan(thx0)/dtan(thy0) endif thx = thx0 thy = -1*thy0 ! reflect in y enp = een0 if (vec(4).gt.rfly) goto 1000 ! apply Ry to neutron else if(yt.ge.tyin1.and.yt.le.tyin2) then ! inside y-walls at end of tar if(xt.lt.txin1) then x = txin1 z = (txin1-x0)/dtan(thx0) + z0 y = y0 + (txin1-x0)*dtan(thy0)/dtan(thx0) else if(xt.gt.txin2) then x = txin2 z = (txin2-x0)/dtan(thx0) + z0 y = y0 + (txin2-x0)*dtan(thy0)/dtan(thx0) endif thx = -1*thx0 ! reflect in x thy = thy0 enp = een0 if (vec(4).gt.rflx) goto 1000 ! apply Rx to neutron else ! outside x-y walls at end of tar if(xt.lt.txin1) then xlm = txin1 else if(xt.gt.txin2) then xlm = txin2 endif if(yt.lt.tyin1) then ylm = tyin1 else if(yt.gt.tyin2) then ylm = tyin2 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 if (vec(4).gt.rflx) goto 1000 ! apply Rx to neutron thy = thy0 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 = een0 endif endif cccccccccc goto 1001 1000 rfail = .true. ! neutron is lost through wall c write(6,*) ' scatt ',rflx,rfly c write(6,*) thx, thy 1001 return end