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
      double precision ix1, ix2, iy1, iy2, iz, gindex
      parameter	(ix1 = -2.85,   ! cm 
c      parameter	(ix1 = -2.3,   ! cm 
     +     ix2 = 2.85,          ! cm
c     +     ix2 = 2.3,          ! cm
     +     iy1 = -2.25,         ! cm
     +     iy2 = 2.25,          ! cm
     +     iz  = 28.0,          ! cm
c     +     iz  = 140.0,          ! cm
     +     gindex = 0.003)      ! (Pyrex glass+ Ti + Ga) rad/Ang
c     +     gindex = 0.001)      ! (Pyrex glass+ Ti + Ga) rad/Ang

c     The real x axis of the supermirror is -3.0 to 3.0 but I subtracted
c     the area of the plates

      lam = 2.*pi*h_par/sqrt(2.*m_n*en0)
      thtest = gindex*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
            thx = thx0
            x = xdis
         endif

         if( (ydis.lt.iy1.or.ydis.gt.iy2) ) then ! cond 4
            if(abs(thy0).le.thtest) then ! cond 5
               if(thy0.ge.0.) then ! cond 6
                  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             ! cond 6
                  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            ! 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 + iz
      else
         x = -1.d6
         y = -1.d6
         z = -1.d6
         thx = -1.d6
         thy = -1.d6
      endif

      return
      end