subroutine alwindow(en0, x0, y0, z0, thx0, thy0, txin1, txin2,
     &     tyin1, tyin2, tzin, enp, x, y, z, thx, thy, scatflg)

      implicit none

      include 'const.inc'

      double precision en0, x0, y0, z0, thx0, thy0, rx0, ry0
	double precision txin1, txin2, tyin1, tyin2, tzin
      double precision x, y, z, thx, thy, ep, enp
	double precision thp, php, sinthp, sinphp
      double precision xdis, ydis
	double precision rfx, rfy, rfz, rdist(3), dreact
	double precision thsc, phsc, thsclab, th0, ph0, tharg, pharg
	double precision sig, thxp, thyp, rand
	real vec(3)
	logical scatflg
c     x0, etc. are current location from main
c     txin1 etc. are constraints of this Al window from main
c      enp,x,y etc are results of passing through this Al window


	    scatflg = .false.

      if((x0.ge.txin1.and.x0.le.txin2).and.
     &   (y0.ge.tyin1.and.y0.le.tyin2)) then   ! we're inside Al

        call ranlux(vec,3)
c	write(6,*) ' rand vector 1, 2',vec(1),vec(2)
c	write(6,*) ' rand vector 3',vec(3)

c       Get polar LAB angles from thx0, thy0
         if(thx0.ge.0..and.thy0.ge.0.) then
            ph0 = datan(dtan(thy0)/dtan(thx0))
         elseif(thx0.lt.0..and.thy0.ge.0.) then
            ph0 = pi + datan(dtan(thy0)/dtan(thx0))
         elseif(thx0.lt.0..and.thy0.lt.0.) then
            ph0 = pi + datan(dtan(thy0)/dtan(thx0))
         elseif(thx0.ge.0..and.thy0.lt.0.) then
            ph0 = 2.*pi + datan(dtan(thy0)/dtan(thx0))
         endif

         th0      = datan(sqrt(dtan(thx0)**2 + dtan(thy0)**2))

c       choose scattering angle -- isotropic in CM - z-axis along
c            neutron momentum
	
		thsc  = dacos(1.d0-2.0*vec(1)) ! isotropic (0,pi) bec
		phsc  = 2*pi*vec(2)

c	 New kinetic energy from LAB scattering angle but still in CS with z along k

        thsclab = datan(dsin(thsc)/(dcos(thsc)+(m_n/m_al)))
	if (thsclab.lt.0.) thsclab = pi + thsclab

	  ep = en0*(m_n/(m_n+m_al))**2*( dcos(thsclab)+
     &       dsqrt( (m_al/m_n)**2-(dsin(thsclab))**2 ) )**2
	if (ep.le.0.d0) write(*,*) ' ep < 0:',ep


c
c      Transform angles back from CS with z along k, to z along beam axis

	tharg = -dsin(th0)*dsin(thsclab)*dcos(phsc)+dcos(th0)*dcos(thsclab)

	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)

	pharg =  (dsin(thsclab)*dcos(phsc)*dcos(ph0)*dcos(th0) -
     &      dsin(thsclab)*dsin(phsc)*dsin(ph0) + 
     &	  dcos(thsclab)*dcos(ph0)*dsin(th0) )/dsin(thp)	

	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(thsclab)*dcos(phsc)*dsin(ph0)*dcos(th0) +
     &      dsin(thsclab)*dsin(phsc)*dcos(ph0) + 
     &      dcos(thsclab)*dsin(ph0)*dsin(th0) )/dsin(thp) 
      if (sinphp.lt.0) php = 2.0d0*pi - php

c
c	write(6,*) ' thsc ', thsc,',  thsclab ',thsclab
c	write(6,*)  ' en0 ', en0,', ep ',ep
c	write(6,*) '    thp ',thp,',  php ',php

c       Get sigma(en0)  en0 in (meV) and find reaction length
c      !!This sigma(E) for Al is fit of intepreted data from ENDF --
c            not much actual data!!!!!

	if (en0.lt.2.2) then
	   sig = 0.1378/(en0/1000.0)**(0.456)    
	 else
	   sig = 1.0602/(en0/1000.0)**(0.129)    
	endif	   

c         sig = 2.6
c
c	write(6,*) 'En ',en0,' sig ',sig

        dreact = -1.0*log(1-vec(3))/(sig*ndenAl)  ! n*sig in 1/cm if sig in barn

c	write(6,*) ' nsig ',sig*ndenAl, ' dreact ',dreact  

            thxp = datan(dtan(thp)*dcos(php))
            thyp = datan(dtan(thp)*dsin(php))
            rfx = dreact*dsin(th0)*dcos(ph0)  ! before scat dist.
            rfy = dreact*dsin(th0)*dsin(ph0)
            rfz = dreact*dcos(th0)
c
c	write(6,*) 'x,y,z',x0,y0,z0
c	write(6,*) 'dx,dy,dz',rfx,rfy,rfz

            rdist(1) = x0 + rfx   ! distance to scatter point
            rdist(2) = y0 + rfy
            rdist(3) = z0 + rfz

c      ****Inner Loop  to determine if scattered, if in target, etc..***	 
       if(  (rdist(1).ge.txin1.and.rdist(1).le.txin2) .and.
     &      (rdist(2).ge.tyin1.and.rdist(2).le.tyin2) .and.
     &      (rdist(3).ge.z0.and.rdist(3).le.(z0+tzin)) ) then  !scat in tar
c       if(  (rdist(3).ge.z0.and.rdist(3).le.(z0+tzin)) ) then  !scat in tar
             
		   xdis  = rdist(1) + (tzin-rfz)*dtan(thxp)		! x loc at end of tar
             ydis  = rdist(2) + (tzin-rfz)*dtan(thyp)		! y loc at end of tar
		   scatflg=.true.
         if( (xdis.ge.txin1.and.xdis.le.txin2) .and.
     &        (ydis.ge.tyin1.and.ydis.le.tyin2) ) then   ! still in Al
            x = xdis
            y = ydis
            z = z0 + tzin
            thx = thxp
            thy = thyp
            enp = ep
         else			!  starts in Al but leaves
c			write(6,*) ' Al scatter - leaves target '
            x = -1.d6
            y = -1.d6
            z = -1.d6
            thx = -1.d6
            thy = -1.d6
            enp = -1.d6
         endif
	 else     ! doesn't scatter

           xdis = x0 + tzin*dtan(thx0)
           ydis = y0 + tzin*dtan(thy0)

         if( (xdis.ge.txin1.and.xdis.le.txin2) .and.
     &        (ydis.ge.tyin1.and.ydis.le.tyin2) ) then   ! still in Al
c		write(6,*) ' Al NO scatter - still in Al '
            x = xdis
            y = ydis
            z = z0 + tzin
            thx = thx0
            thy = thy0
            enp = en0
         else			!  starts in Al but leaves
c			write(6,*) ' Al NO scatter - leaves Al '
            x = -1.d6
            y = -1.d6
            z = -1.d6
            thx = -1.d6
            thy = -1.d6
            enp = -1.d6
         endif
	 endif
      else				!  never made it into Al window
c			write(6,*) ' Al Never in Al '
         x = -1.d6
         y = -1.d6
         z = -1.d6
         thx = -1.d6
         thy = -1.d6
         enp = -1.d6
      endif

      return
      end