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

      implicit none

      include 'const.inc'

      integer i
      double precision en0, lam
      double precision x0, y0, z0, thx0, thy0, x, y, z, thx, thy
      double precision gdx1, gdx2, gdx3, gdx4
      double precision thtest, testx, xdis, ydis, iz0
      double precision ix1, ix2, iy1, iy2, iz, gl, gin
      parameter	(ix1 = -2.80,   ! cm
     +     ix2 = 2.80,          ! cm
c     +     iy1 = -2.30,         ! cm  
c     +     iy2 = 2.30)          ! cm
     +     iy1 = -2.80,         ! cm  
     +     iy2 = 2.80)          ! cm
c
      iz = gl                   ! cm
      lam = 2.*pi*h_par/sqrt(2.*m_n*en0)
      thtest = gin*lam
c
      if((y0.ge.iy1.and.y0.le.iy2).and.(x0.ge.ix1.and.x0.le.ix2) ) then
         xdis = x0 + iz*dtan(thx0)
         ydis = y0 + iz*dtan(thy0)
         if(xdis.lt.ix1.or.xdis.gt.ix2) then
            if(abs(thx0).le.thtest) then
               if(thx0.ge.0.) then
                  thx = -1.*thx0
                  x = 2.*ix2 - iz*dtan(thx0) - x0
                  if(x.lt.ix1.or.x.gt.ix2) then
                     thx = thx0
                     x = 2.*ix1 - 2.*ix2 + iz*dtan(thx0) + x0
                  endif
               else
                  thx = -1.*thx0
                  x   = 2.*ix1 - iz*dtan(thx0) - x0
                  if(x.lt.ix1.or.x.gt.ix2) then
                     thx = thx0
                     x = 2.*ix2 - 2.*ix1 + iz*dtan(thx0) + x0
                  endif
               endif
            else
               x=-1.d6
               thx=-1.d6
            endif
         else
            x=xdis
            thx=thx0
         endif
         if(ydis.lt.iy1.or.ydis.gt.iy2) then
            if(abs(thy0).le.thtest) then
               if(thy0.ge.0.) then
                  thy = -1.*thy0
                  y   = 2.*iy2 - iz*dtan(thy0) - y0
                  if(y.lt.iy1.or.y.gt.iy2) then
                     thy = thy0
                     y = 2.*iy1 - 2.*iy2 + iz*dtan(thy0) + y0
                  endif
               else
                  thy = -1.*thy0
                  y   = 2.*iy1 - iz*dtan(thy0) - y0
                  if(y.lt.iy1.or.y.gt.iy2) then
                     thy = thy0
                     y = 2.*iy2 - 2.*iy1 + iz*dtan(thy0) + y0
                  endif
               endif
            else
               y=-1.d6
               thy=-1.d6
            endif
         else
            y=ydis
            thy=thy0
         endif
         z=z0+iz
      else
         x=-1.d6
         y=-1.d6
         z=-1.d6
         thx=-1.d6
         thy=-1.d6
      endif

      return
      end