subroutine scate(en0, x0, y0, z0, thx0, thy0,
     &     txin1,txin2,tyin1,tyin2,zlm,
     &     xsc, ysc, zsc, scattFlag,
     &     enp, x, y, z, thx, thy,xstot,thsc,lam,rfail)

      implicit none

      include 'const.inc'

      integer hws, scattFlag
      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 txin1, txin2, tyin1, tyin2, 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 (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*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 (airtar) then  ! air in empty and target cells
	
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 = ndenairtar
			norm = 1.0
		else		! normally He gas in empty cell
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 = ndenhegas
			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 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*** 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 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
            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 = en0
         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 = en0    
		  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 = en0
		  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
               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