subroutine target(q, T, en0, x0, y0, z0, thx0, thy0, tx1, tx2,
     &     ty1, ty2, tz, enp, x, y, z, thx, thy)

      implicit none

      include 'const.inc'

      integer i
      real vec1(1), vec(3), qrand(1), HRNDM1
      double precision en0, x0, y0, z0, thx0, thy0, tx1, tx2
      double precision ty1, ty2, tz, x, y, z, thx, thy, ep, enp
      double precision rfx, rfy, rfz
      double precision dxs, T, hedens, cs, qp, eqp, th_test
      double precision thsc, phsc, thxsc, thysc
      double precision henumb, dreact, domga
      double precision th0, ph0, domega, rdist(3)
      double precision k0, kp, thp, php, thxp, thyp, xdis, ydis
      double precision qp1, qp2, csthsc1, csthsc2
      integer qmin, qmax
      double precision q, eq, wmin, wmax, wtest, xst, xstot
      double precision denat, rc, rg, dsig0, ratio, thxtst, thytst
      double precision c1,c2,c3,con1,con2,con3,mm, thp1, thp2
      double precision norm, wp
c
      if((x0.ge.tx1.and.x0.le.tx2).and.(y0.ge.ty1.and.y0.le.ty2)) then

         if(thx0.ge.0..and.thy0.ge.0.) then
            ph0 = datan(dtan(thy0)/dtan(thx0))
         elseif(thx0.lt.0..and.thy0.ge.0.) then
            ph0 = pi + datan(dtan(thy0)/dtan(thx0))
         elseif(thx0.lt.0..and.thy0.lt.0.) then
            ph0 = pi + datan(dtan(thy0)/dtan(thx0))
         elseif(thx0.ge.0..and.thy0.lt.0.) then
            ph0 = 2.*pi + datan(dtan(thy0)/dtan(thx0))
         endif

         th0      = datan(sqrt(dtan(thx0)**2 + dtan(thy0)**2))

         k0 = sqrt(2*m_n*en0)/h_par
         call wq(q,wp)
         ep  = en0- (h_par*wp)
         if(en0.le.1.) then
            xdis = x0 + tz*dtan(thx0)
            ydis = y0 + tz*dtan(thy0)
            if( (xdis.ge.tx1.and.xdis.le.tx2) .and.
     &           (ydis.ge.ty1.and.ydis.le.ty2) ) then
               x   = xdis
               y   = ydis
               z   = z0 + tz
               thx = thx0
               thy = thy0
               enp = en0
            else
               x   = -1.d6
               y   = -1.d6
               z   = -1.d6
               thx = -1.d6
               thy = -1.d6
               enp = -1.d6
            endif
         else
            kp   = sqrt(2*m_n*ep)/h_par
            th_test=(k0**2 + kp**2 - q**2)/(2.*k0*kp)
            if(th_test.gt.1.0) write(*,*) '----------ooo----------'
            thsc = dacos(th_test)
            call ranlux(vec,3)
            phsc = 2.0*pi*vec(1)
            php  = datan(dtan(phsc)*dcos(thx0)/dcos(thy0))
            thp  = dasin(dsin(thsc)*dcos(phsc)/(dcos(php)*dcos(thx0)*
     &           dcos(th0)))
c     Cross section .........................
            if(q.le.0.562803) then
               call intq(T, k0, q, thsc, xst)
               xstot = xst
            else
               domga = phsc*(1.-dcos(thsc))/(2.*pi)
               call sommers(T, en0, dxs)
               xstot = domga*dxs/3.33
            endif
c     ............................................
            call heden(T,hedens)
            denat = hedens*avog/he4mol
            call nden(T,norm)


            dreact = -1.0*log(vec(2))/(denat*xstot*norm)
 
            thxp = datan(dtan(thp)*dcos(php))
            thyp = datan(dtan(thp)*dsin(php))
            rfx = dreact*dsin(th0)*dcos(ph0)
            rfy = dreact*dsin(th0)*dsin(ph0)
            rfz = dreact*dcos(th0)

            rdist(1) = x0 + rfx
            rdist(2) = y0 + rfy
            rdist(3) = z0 + rfz
            if(  (rdist(1).ge.tx1.and.rdist(1).le.tx2) .and.
     &           (rdist(2).ge.ty1.and.rdist(2).le.ty2) .and.
     &           (rdist(3).ge.z0.and.rdist(3).le.(z0+tz)) ) then
               xdis  = rdist(1) + (tz-rfz)*dtan(thxp)
               ydis  = rdist(2) + (tz-rfz)*dtan(thyp)
               if( (xdis.ge.tx1.and.xdis.le.tx2) .and.
     &              (ydis.ge.ty1.and.ydis.le.ty2) ) then
                  x = xdis
                  y = ydis
                  z = z0 + tz
                  thx = thxp
                  thy = thyp
                  enp = ep
               else
                  x = -1.d6
                  y = -1.d6
                  z = -1.d6
                  thx = -1.d6
                  thy = -1.d6
                  enp = -1.d6
               endif
            else
               xdis = x0 + tz*dtan(thx0)
               ydis = y0 + tz*dtan(thy0)
               if( (xdis.ge.tx1.and.xdis.le.tx2) .and.
     &              (ydis.ge.ty1.and.ydis.le.ty2) ) then
                  x = xdis
                  y = ydis
                  z = z0 + tz
                  thx = thx0
                  thy = thy0
                  enp = en0
               else
                  x = -1.d6
                  y = -1.d6
                  z = -1.d6
                  thx = -1.d6
                  thy = -1.d6
                  enp = -1.d6
               endif
            endif
         endif
      else
         x = -1.d6
         y = -1.d6
         z = -1.d6
         thx = -1.d6
         thy = -1.d6
         enp = -1.d6
      endif

      return
      end