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

      lam = 2.*pi*h_par/sqrt(2.*m_n*en0)
      thtest = gsm*lam
c
      if((y0.ge.sy1.and.y0.le.sy2).and.(x0.ge.sx1.and.x0.le.sx2)) then
         xdis = x0 + sz*dtan(thx0)
         ydis = y0 + sz*dtan(thy0)
         if( (xdis.lt.sx1.or.xdis.gt.sx2) ) then
            if(abs(thx0).le.thtest) then
               if(thx0.ge.0.) then
                  thx = -1.*thx0
                  x   = 2.*sx2 - sz*dtan(thx0) - x0
                  if(x.lt.sx1.or.x.gt.sx2) then
                     thx = thx0
                     x = 2.*sx1 - 2.*sx2 + sz*dtan(thx0) + x0
                  endif
               else
                  thx = -1.*thx0
                  x   = 2.*sx1 - sz*dtan(thx0) - x0
                  if(x.lt.sx1.or.x.gt.sx2) then
                     thx = thx0
                     x = 2.*sx2 - 2.*sx1 + sz*dtan(thx0) + x0
                  endif
               endif
            else
               x   = -1.d6
               thx = -1.d6
            endif
         else
            thx = thx0
            x = xdis
         endif

         if( (ydis.lt.sy1.or.ydis.gt.sy2) ) then ! cond 4
            if(abs(thy0).le.thtest) then ! cond 5
               if(thy0.ge.0.) then ! cond 6
                  thy = -1.*thy0
                  y   = 2.*sy2 - sz*dtan(thy0) - y0
                  if(y.lt.sy1.or.y.gt.sy2) then
                     thy = thy0
                     y = 2.*sy1 - 2.*sy2 + sz*dtan(thy0) + y0
                  endif
               else             ! cond 6
                  thy = -1.*thy0
                  y   = 2.*sy1 - sz*dtan(thy0) - y0
                  if(y.lt.sy1.or.y.gt.sy2) then
                     thy = thy0
                     y = 2.*sy2 - 2.*sy1 + sz*dtan(thy0) + y0
                  endif
               endif            ! cond 6
            else                ! cond 5
               y = -1.d6
               thy = -1.d6
            endif               ! cond 5
         else                   ! cond 4
            thy = thy0
            y = ydis
         endif                  ! cond 4
         z = z0 + sz
      else
         x = -1.d6
         y = -1.d6
         z = -1.d6
         thx = -1.d6
         thy = -1.d6
      endif

      return
      end