subroutine guide(en0,x0,y0,z0,thx0,thy0,x,y,z,thx,thy,nxr,nyr)

      implicit none

      include 'const.inc'

      integer i, nxr, nyr
      double precision en0, lam
      double precision x0, y0, z0, thx0, thy0, x, y, z, thx, thy
      double precision gxt, gxtp, gyt, gytp, iz, xdis, ydis
      double precision thxtest, thytest, rflx,rfly
	double precision, external :: rsqrt,ratan
	real vec(2)
c
	  call ranlux(vec,2)

	nxr=0;nyr=0
      iz = gl                   ! cm
      lam = 2.*pi*h_par/rsqrt(2.*m_n*en0) ! Ang
      thxtest = ginx*lam
      thytest = giny*lam
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 ((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
cccc   y
		  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.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
	end if
ccccc
 
c   non-bounce final location in x and y
         xdis = x0 + iz*dtan(thx0)
         ydis = y0 + iz*dtan(thy0)

      
      if (sptmflg.eq..true.) then          !septum -rcm
      if((y0.ge.gy1.and.y0.le.gy2).and.((x0.ge.gx1.and.x0.le.gx2)
     &	.or.(x0.ge.gx3.and.x0.le.gx4))) then
c
! 1st do case of neutron making it to the end of the guide regardless of
!  it's entering angle, i.e., despite large angle never hits wall
!  this is the largest number of neutrons for NSRII setup

	if( (((x0.gt.gx1.and.x0.lt.gx2).and.(xdis.gt.gx1.and.xdis.lt.gx2))
     &.or.((x0.gt.gx3.and.x0.lt.gx4).and.(xdis.gt.gx3.and.xdis.lt.gx4)))
     & .and.(ydis.gt.gy1.and.ydis.lt.gy2) ) then
		x=xdis
		y=ydis   
	    z = z0 + iz
!	 now do bounces, including multiple bounces
!      x,y separately since might only bounce on one wall, i.e., might
!      bounce with ok angle in x and make it to the end despite a large
!      angle in y
  
        else 
c     x-axis
	   if (x0.ge.gx1.and.x0.le.gx2) then
		if (abs(thx0).le.thxtest) then
	      if(thx0.ge.0) then
              gxt = gx2
              gxtp = gx1
            else
              gxt = gx1
              gxtp = gx2
            endif
            nxr = (iz*dtan(thx0)+x0-gxtp)/(gxt-gxtp) ! number of x-reflections
            if(nxr.eq.1) x = 2*gxt-x0-iz*dtan(thx0)    !corrected x-position -rcm
	      if(nxr.eq.2) x = x0+iz*dtan(thx0)-2*(gxt-gxtp)
	      if(nxr.eq.3) x = 2*gxt-x0-iz*dtan(thx0)+2*(gxt-gxtp)
            thx = (-1)**nxr*thx0
ccccc   reflectivity
				rflx = rflx**nxr
c	write(6,*) thxtest,thx0,nxr,rflx
ccccc
	     else
		  x=xdis
	     end if
	   else
           if (abs(thx0).le.thxtest) then      !septum -rcm
	      if(thx0.ge.0) then
              gxt = gx4
              gxtp = gx3
            else
              gxt = gx3
              gxtp = gx4
            endif
            nxr = (iz*dtan(thx0)+x0-gxtp)/(gxt-gxtp) ! number of x-reflections
            if(nxr.eq.1) x = 2*gxt-x0-iz*dtan(thx0)    !corrected x-position -rcm
	      if(nxr.eq.2) x = x0+iz*dtan(thx0)-2*(gxt-gxtp)
	      if(nxr.eq.3) x = 2*gxt-x0-iz*dtan(thx0)+2*(gxt-gxtp)
            thx = (-1)**nxr*thx0

ccccc   reflectivity
				rflx = rflx**nxr
c	write(6,*) thxtest,thx0,nxr,rflx
ccccc
	     else
		  x=xdis
	     end if
	    end if
c     y-axis
	    if (abs(thy0).le.thytest) then
	      if(thy0.ge.0) then
               gyt = gy2
               gytp = gy1
           else
              gyt = gy1
              gytp = gy2
           endif
           nyr = (iz*dtan(thy0)+y0-gytp)/(gyt-gytp) ! number of y-reflections
           if(nyr.eq.1) y = 2*gyt-y0-iz*dtan(thy0)    !corrected y-position -rcm
	     if(nyr.eq.2) y = y0+iz*dtan(thy0)-2*(gyt-gytp)
	     if(nyr.eq.3) y = 2*gyt-y0-iz*dtan(thy0)+2*(gyt-gytp)
           thy = (-1)**nyr*thy0
ccccc   reflectivity
			rfly = rfly**nyr
c	write(6,*) thytest,thy0,nyr,rfly
ccccc
	    else
	     y=ydis
	    end if

c         write(35,*) x0, thx0,y0,thy0

	    if ( ((x.gt.gx1.and.x.lt.gx2).or.(x.gt.gx3.and.x.lt.gx4))
     &        .and.(y.gt.gy1.and.y.lt.gy2) ) then
c     z-axis
		  if (vec(1).le.rflx*rfly) then  !  apply Rx*Ry to neutron
               z=z0+iz
		  else
		   nxr=0
		   nyr=0
             x = -1.d6
             y = -1.d6
             z = -1.d6
             thx = -1.d6
             thy = -1.d6	
	      end if
	     else
		   nxr=0
		   nyr=0
             x = -1.d6
             y = -1.d6
             z = -1.d6
             thx = -1.d6
             thy = -1.d6
	     end if
	  end if


      else
	   nxr=0
	   nyr=0
         x = -1.d6
         y = -1.d6
         z = -1.d6
         thx = -1.d6
         thy = -1.d6
      endif

      end if
	
************This is for no septum	-rcm*********************
	if (sptmflg.eq..false.) then 
	      if((y0.ge.gy1.and.y0.le.gy2).and.(x0.ge.gx1.and.x0.le.gx4)) then
c
! 1st do case of neutron making it to the end of the guide regardless of
!  it's entering angle, i.e., despite large angle never hits wall
!  this is the largest number of neutrons for NSRII setup
	  if ( (xdis.gt.gx1.and.xdis.lt.gx4).and.
     &       (ydis.gt.gy1.and.ydis.lt.gy2) ) then
		x=xdis
		y=ydis   
	    z = z0 + iz
!	 now do bounces, including multiple bounces
!      x,y separately since might only bounce on one wall, i.e., might
!      bounce with ok angle in x and make it to the end despite a large
!      angle in y
  
        else 
c     x-axis
	    if (abs(thx0).le.thxtest) then
	      if(thx0.ge.0) then
              gxt = gx4
              gxtp = gx1
            else
              gxt = gx1
              gxtp = gx4
            endif
            nxr = (iz*dtan(thx0)+x0-gxtp)/(gxt-gxtp) ! number of x-reflections
            if(nxr.eq.1) x = 2*gxt-x0-iz*dtan(thx0)    !corrected x-position -rcm
	      if(nxr.eq.2) x = x0+iz*dtan(thx0)-2*(gxt-gxtp)
	      if(nxr.eq.3) x = 2*gxt-x0-iz*dtan(thx0)+2*(gxt-gxtp)
            thx = (-1)**nxr*thx0
ccccc   reflectivity
				rflx = rflx**nxr
c	write(6,*) thxtest,thx0,nxr,rflx
ccccc
	     else
		  x=xdis
	     end if

c     y-axis
	    if (abs(thy0).le.thytest) then
	      if(thy0.ge.0) then
               gyt = gy2
               gytp = gy1
           else
              gyt = gy1
              gytp = gy2
           endif
           nyr = (iz*dtan(thy0)+y0-gytp)/(gyt-gytp) ! number of y-reflections
           if(nyr.eq.1) y = 2*gyt-y0-iz*dtan(thy0)    !corrected y-position -rcm
	     if(nyr.eq.2) y = y0+iz*dtan(thy0)-2*(gyt-gytp)
	     if(nyr.eq.3) y = 2*gyt-y0-iz*dtan(thy0)+2*(gyt-gytp)
           thy = (-1)**nyr*thy0
ccccc   reflectivity
			rfly = rfly**nyr
c	write(6,*) thytest,thy0,nyr,rfly
ccccc
	    else
	     y=ydis
	    end if

c         write(35,*) x0, thx0,y0,thy0

	    if ((x.gt.gx1.and.x.lt.gx4).and.
     &        (y.gt.gy1.and.y.lt.gy2) ) then
c     z-axis
		  if (vec(1).le.rflx*rfly) then  !  apply Rx*Ry to neutron
               z=z0+iz
		  else
		   nxr=0
		   nyr=0
             x = -1.d6
             y = -1.d6
             z = -1.d6
             thx = -1.d6
             thy = -1.d6	
	      end if
	     else
		   nxr=0
		   nyr=0
             x = -1.d6
             y = -1.d6
             z = -1.d6
             thx = -1.d6
             thy = -1.d6
	     end if
	  end if


      else
	   nxr=0
	   nyr=0
         x = -1.d6
         y = -1.d6
         z = -1.d6
         thx = -1.d6
         thy = -1.d6
      endif
      
      end if

c	if(((x0*x).lt.0).and.(x.gt.-1d5)) then   !checks for crossovers
c      write(6,*) ' x0 ',x0,' thx0 ',thx0
c	write(6,*) ' x1 ',x,' nxr ',nxr
c      end if

	return
      end