subroutine beampipe(beamflg, en0, x0, y0, z0, 
     &      thx0, thy0, xin1, xin2, yin1, yin2, zlm, enp, x, y, z,
     &      thx, thy, scflag, nbx, nby, pl, totrot) 
  
      implicit none

      include 'const.inc'

      integer scf, scflag, nbx, nby, beamflg
      integer i, j, n, k
      parameter(n=50)
      double precision q, en0, x0, y0, z0, thx0, thy0
      double precision lam, thtestx, thtesty, zin1 
      double precision xin1, xin2, yin1, yin2, zlm, pl, totrot, rotang
      double precision x, y, z, thx, thy, ep, enp
      double precision xi(n), yi(n), zi(n), thxi(n), thyi(n), eni(n)
      double precision xs(n), ys(n), zs(n)
      double precision xo(n), yo(n), zo(n), thxo(n), thyo(n), eno(n)
	double precision xsig,thscat,pl_sc,pl_o,pl_ns
	double precision, external :: rsqrt,ratan
	logical rfailflg

c 
      lam = 2.*pi*h_par/rsqrt(2.*m_n*en0) ! Ang
c
	rfailflg = .false.
	xsig = -1.d6
      thscat  = -1.d6
      scf = 0
      scflag = 0
      totrot = 0.0d0
	zin1 = z0
      nbx=0
	nby=0

	if(beamflg.eq.1) then
	  thtestx = ginx*lam
	  thtesty = giny*lam
        zlm = z0 + gl      

	else if(beamflg.eq.2) then
	  thtestx = ginx*lam
	  thtesty = giny*lam
        zlm = z0 + iz0

	else if(beamflg.eq.3) then
	  thtestx = gotx*lam
	  thtesty = goty*lam
        zlm = z0 + oz
	 
	end if   

      if((x0.ge.xin1.and.x0.le.xin2).and.
     &     (y0.ge.yin1.and.y0.le.yin2)) then
         xi(1) = x0
         yi(1) = y0
         zi(1) = z0
         thxi(1) = thx0
         thyi(1) = thy0
         eni(1) = en0
         pl = 0

         do j=1, n		! loop through multiple scatterings
            if (j.ge.2) then
               xi(j) = xo(j-1)	! starting position now at last ending position
               yi(j) = yo(j-1)
               zi(j) = zo(j-1)
               thxi(j) = thxo(j-1)
               thyi(j) = thyo(j-1)
               eni(j) = eno(j-1)
            endif


C  load current position into scatb.f subroutine, check for scattering, do 
C    reflections and return new positions, angles, energy
            call scatb(beamflg, eni(j), 
     &           xi(j), yi(j), zi(j),thxi(j), thyi(j), 
     &           xin1, xin2, yin1, yin2, zlm,xs(j), ys(j), zs(j), scf,
     &           eno(j), xo(j), yo(j), zo(j), thxo(j), thyo(j), 
     &           xsig,thscat,lam,rfailflg)

		if (rfailflg) goto 110		! neutron is lost through wall
		if (thxo(j).eq.-1*thxi(j)) nbx = nbx + 1
 		if (thyo(j).eq.-1*thyi(j)) nby = nby + 1
		scflag = scflag + scf  !  
c         scf is used in scatb.f to indicate if scatter position is inside
c             target (1) or outside target (0)

c      update rotation  -- need begining and ending positions and energy
		if (scf.eq.1) then
        pl_sc =rsqrt((xs(j)-xi(j))**2+(ys(j)-yi(j))**2+(zs(j)-zi(j))**2) 
	  rotang = gy*mag_b_beam*pl_sc*0.1/rsqrt(2.*eni(j)/m_n)
	  pl = pl_sc + pl
		   totrot = totrot + rotang
	  pl_o = rsqrt((xo(j)-xs(j))**2+(yo(j)-ys(j))**2+(zo(j)-zs(j))**2)
	  rotang = gy*mag_b_beam*pl_o*0.1/rsqrt(2.*eno(j)/m_n)
	  pl = pl_o + pl
		   totrot = totrot + rotang
		else
			if (eno(j).ne.eni(j)) write(6,*) ' no scat but Ei.ne.Eo '
	  pl_ns=rsqrt((xo(j)-xi(j))**2+(yo(j)-yi(j))**2+(zo(j)-zi(j))**2) 
              rotang = gy*mag_b_beam*pl_o*0.1/rsqrt(2.*eno(j)/m_n)
	  pl=pl_ns + pl
			totrot = totrot + rotang
		end if

c		if (j.gt.1) then
c		write(6,*) ' empty ', j	
c		write(6,*) xo(j),yo(j),zo(j)
c		write(6,*) scf, rotang, totrot
c	    write(6,*) ' pl ',pl
c		end if

cccc		end rotation update

c   below is somewhat redundant because of rflx,y in scate.f, but if you set 
c       rflx,y=1 in scatb.f, then this is needed.  So leave in.	
	   if ((thxo(j).eq.-1*thxi(j)).and.(abs(thxo(j)).gt.thtestx)) then
            x = -1.d6
            y = -1.d6
            z = -1.d6
            thx = -1.d6
            thy = -1.d6
            enp = -1.d6         
	      exit
         endif
	   if ((thyo(j).eq.-1*thyi(j)).and.(abs(thyo(j)).gt.thtesty)) then
            x = -1.d6
            y = -1.d6
            z = -1.d6
            thx = -1.d6
            thy = -1.d6
            enp = -1.d6         
	      exit
         endif
c
           
	if (j.gt.49) write(6,*) ' j,en, xsig ',j,eni(j), xsig
	if (j.gt.49) write(6,*) ' x,y,z ',xo(j),yo(j),zo(j)


            if(zo(j)-zlm.ge.0) then
               x = xo(j)
               y = yo(j)
               z = zo(j)
               thx = thxo(j)
               thy = thyo(j)
               enp = eno(j)
               exit	!  Jump out of loop
            endif
         enddo

      else
         x = -1.d6
         y = -1.d6
         z = -1.d6
         thx = -1.d6
         thy = -1.d6
         enp = -1.d6
      endif

	goto 111

 110   continue
	   x = -1.d6
         y = -1.d6
         z = -1.d6
         thx = -1.d6
         thy = -1.d6
         enp = -1.d6

 111      return
      end