subroutine outcoil(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

c
      lam = 2.*pi*h_par/sqrt(2.*m_n*en0)
      thtest = got*lam
c
      if((y0.ge.iy1.and.y0.le.iy2).and.((x0.ge.ox1.and.x0.le.ox2)
     &     .or.(x0.ge.ox3.and.x0.le.ox4))) then

         xdis = x0 + oz*dtan(thx0)
         ydis = y0 + oz*dtan(thy0)

         if((x0.ge.ox3.and.x0.le.ox4).and.
     &        (xdis.lt.ox3.or.xdis.gt.ox4)) then
            if(abs(thx0).le.thtest) then
               if(thx0.ge.0) then
                  thx = -1.*thx0
                  x = 2.*ox4 - oz*dtan(thx0) - x0
                  if(x.lt.ox3) then
                     thx = thx0
                     x = 2.*ox3 - 2.*ox4 + oz*dtan(thx0) + x0
                  endif
               else
                  thx = -1.*thx0
                  x = 2.*ox3 - oz*dtan(thx0) - x0
                  if(x.gt.ox4) then
                     thx = thx0
                     x   = 2.*ox4 - 2.*ox3 + oz*dtan(thx0) + x0
                  endif
               endif
            else
               x   = -1.d6
               thx = -1.d6
            endif
         elseif((x0.ge.ox1.and.x0.le.ox2).and.
     &           (xdis.lt.ox1.or.xdis.gt.ox2)) then
            if(abs(thx0).le.thtest) then
               if(thx0.ge.0) then
                  thx = -1.*thx0
                  x = 2.*ox2 - oz*dtan(thx0) - x0
                  if(x.lt.ox1) then
                     thx = thx0
                     x = 2.*ox1 - 2.*ox2 + oz*dtan(thx0) + x0
                  endif
               else
                  thx = -1.*thx0
                  x = 2.*ox1 - oz*dtan(thx0) - x0
                  if(x.gt.ox2) then
                     thx = thx0
                     x   = 2.*ox2 - 2.*ox1 + oz*dtan(thx0) + x0
                  endif
               endif
            else
               x   = -1.d6
               thx = -1.d6
            endif
         else
            thx = thx0
            x = xdis
         endif
         if( (y0.ge.iy1.and.y0.le.iy2) .and.
     &        (ydis.lt.iy1.or.ydis.gt.iy2) ) then
            if(abs(thy0).le.thtest) then
                  thy = -1.*thy0
               if(thy0.ge.0.) then
                  y   = 2.*iy2 - oz*dtan(thy0) - y0
               else
                  y   = 2.*iy1 - oz*dtan(thy0) - y0
               endif
            else
               y   = -1.d6
               thy = -1.d6
            endif
         else
            thy = thy0
            y = ydis
         endif
         z = z0 + oz
      else
         x = -1.d6
         y = -1.d6
         z = -1.d6
         thx = -1.d6
         thy = -1.d6
      endif

      return
      end