subroutine guide(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, iz
c
      iz = gl                   ! cm
      lam = 2.*pi*h_par/sqrt(2.*m_n*en0)
      thtest = gin*lam
c
      if((y0.ge.gy1.and.y0.le.gy2).and.(x0.ge.gx1.and.x0.le.gx2) ) then
         xdis = x0 + iz*dtan(thx0)
         ydis = y0 + iz*dtan(thy0)
         if(xdis.lt.gx1.or.xdis.gt.gx2) then
            if(abs(thx0).le.thtest) then
               if(thx0.ge.0.) then
                  thx = -1.*thx0
                  x = 2.*gx2 - iz*dtan(thx0) - x0
                  if(x.lt.gx1.or.x.gt.gx2) then
                     thx = thx0
                     x = 2.*gx1 - 2.*gx2 + iz*dtan(thx0) + x0
                  endif
               else
                  thx = -1.*thx0
                  x   = 2.*gx1 - iz*dtan(thx0) - x0
                  if(x.lt.gx1.or.x.gt.gx2) then
                     thx = thx0
                     x = 2.*gx2 - 2.*gx1 + iz*dtan(thx0) + x0
                  endif
               endif
            else
               x=-1.d6
               thx=-1.d6
            endif
         else
            x=xdis
            thx=thx0
         endif
         if(ydis.lt.gy1.or.ydis.gt.gy2) then
            if(abs(thy0).le.thtest) then
               if(thy0.ge.0.) then
                  thy = -1.*thy0
                  y   = 2.*gy2 - iz*dtan(thy0) - y0
                  if(y.lt.gy1.or.y.gt.gy2) then
                     thy = thy0
                     y = 2.*gy1 - 2.*gy2 + iz*dtan(thy0) + y0
                  endif
               else
                  thy = -1.*thy0
                  y   = 2.*gy1 - iz*dtan(thy0) - y0
                  if(y.lt.gy1.or.y.gt.gy2) then
                     thy = thy0
                     y = 2.*gy2 - 2.*gy1 + 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