subroutine incoil(en0,x0,y0,z0,thx0,thy0,x,y,z,thx,thy)


      implicit none

      include 'const.inc'

      integer i
      double precision en0, en, 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

      iz = iz0                  ! cm
      lam = 2.*pi*h_par/sqrt(2.*m_n*en0)
      thtest = gin*lam
      en=en0
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   !  refl. so return new angle as -old angle
                  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  ! reflected twice so new ang = old ang
                     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