subroutine target(q, en0, x0, y0, z0, thx0, thy0, txin1, txin2,
     &     enp, x, y, z, thx, thy,xstot,
     &	 thsc, dxs, xst,enflg,thflg,thp,php,j,pla,plb) !bec

      implicit none

      include 'const.inc'

      integer i,j,loop		! bec
      real vec1(1), vec(3), qrand(1), HRNDM1
      double precision en0, x0, y0, z0, thx0, thy0
	double precision txin1, txin2, pla, plb
      double precision x, y, z, thx, thy, ep, enp
      double precision rfx, rfy, rfz
      double precision dxs, 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, sinphp, tharg, pharg
	logical enflg,thflg
c 
	    enflg = .false.
	    thflg = .false.

      if((x0.ge.txin1.and.x0.le.txin2).and.
     &   (y0.ge.tary1.and.y0.le.tary2)) 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,en0,k0)
         ep  = en0- (h_par*wp)

	if (ep.le.0.d0) write(*,*) ' ep < 0:',ep
c
c		en0<1 can lead to unphysical results: ep<0 or unphysical scattering
c		 angle depending on what q was chosen.  Now, in main_cwn.f the code
c		 goes back and randomly chooses q again and again until physical results.
c          In addition, if q<0.1, wq subroutine decreases omega from table value until
c          the angle can be calculated (momentum conservation)  - bec
c 
c        if(en0.le.1.) then    
c            xdis = x0 + tarz*dtan(thx0)
c            ydis = y0 + tarz*dtan(thy0)
c            if( (xdis.ge.txin1.and.xdis.le.txin2) .and.
c     &           (ydis.ge.tary1.and.ydis.le.tary2) ) then
c               x   = xdis
c               y   = ydis
c               z   = z0 + tarz
c               thx = thx0
c               thy = thy0
c               enp = en0
c            else
c               x   = -1.d6
c               y   = -1.d6
c               z   = -1.d6
c               thx = -1.d6
c               thy = -1.d6
c               enp = -1.d6
c            endif
c         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) then
			th_test = 1.0d0
			if (th_test.gt.1.05) then
				write(*,*) '----------ooo----------'
				write(*,*) 'k0:',k0,'   kp:',kp,'   q:',q
			endif
		  endif
c	write(6,*) ' target thsc '
c	write(6,*) th_test, k0, kp, ep
            thsc = dacos(th_test)
            call ranlux(vec,3)
            phsc = 2.0*pi*vec(1)
cc            php  = datan(dtan(phsc)*dcos(thx0)/dcos(thy0))
cc      by definition cos(thx0)=cos(thy0) 
cc		  php = phsc    ! bec
cc	if (abs(dsin(thsc)*dcos(phsc)/(dcos(php)*dcos(thx0)*dcos(th0)))
cc     &	.gt.1.0) write(*,*) 'thsc:',thsc,'  phsc:',phsc,'  php:',php
cc     &  ,'  arg:',dsin(thsc)*dcos(phsc)/(dcos(php)*dcos(thx0)*dcos(th0))
cc            thp  = dasin(dsin(thsc)*dcos(phsc)/(dcos(php)*dcos(thx0)*
cc     &           dcos(th0)))
c
c	Using transformation matrix from x,y,z to X,Y,Z (lab) where k-vect is along 
c        z-axis.  Rotate by phi about Z, then theta about y, and chi=0 about z (no need for this rot).
c        This rotation creates CS, x,y,z relative to X,Y,Z, where theta and phi locate z axis
c        in X,Y,Z (lab) frame.  We'll take k-vector along z, so 
c        theta_o and phi_o equal theta and phi of rotated frame.  kp (scattered vec) is 
c        theta_sc from z (=k) with azimuthal angle phi_sc.  Use following to rotate back to X,Y,Z, giving 
c        theta_p and phi_p in X,Y,Z.


	tharg = -dsin(th0)*dsin(thsc)*dcos(phsc)+dcos(th0)*dcos(thsc)

	if (tharg.gt.1.0d0) then
		write(*,*) ' theta_p error :  cos(th_p) =',tharg
		tharg=1.0d0
	else if (tharg.lt.-1.0d0) then
		write(*,*) ' theta_p error :  cos(th_p) =',tharg
		tharg=-1.0d0
	endif
c		write(6,*) ' target tharg '
	thp = dacos(tharg)

	pharg =  (dsin(thsc)*dcos(phsc)*dcos(ph0)*dcos(th0) -
     &      dsin(thsc)*dsin(phsc)*dsin(ph0) + 
     &	  dcos(thsc)*dcos(ph0)*dsin(th0) )/dsin(thp)	

	if (pharg.gt.1.0d0) then
		write(*,*) ' phi_p error :  cos(ph_p) =',pharg
		pharg=1.0d0
	else if (pharg.lt.-1.0d0) then
		write(*,*) ' phi_p error :  cos(ph_p) =',pharg
		pharg=-1.0d0
	end if

c	write(6,*) ' target pharg '
		php = dacos(pharg)
	sinphp =     (dsin(thsc)*dcos(phsc)*dsin(ph0)*dcos(th0) +
     &      dsin(thsc)*dsin(phsc)*dcos(ph0) + 
     &      dcos(thsc)*dsin(ph0)*dsin(th0) )/dsin(thp) 
      if (sinphp.lt.0) php = 2.0d0*pi - php

c     Cross section .........................
            if(q.le.0.562803) then
               call intq(k0, q, thsc, xst)
               xstot = xst
            else
c               domga = phsc*(1.-dcos(thsc))/(2.*pi)  !  don't need -- bec
               call sommers(T, en0, dxs)
c               xstot = domga*dxs/3.33   !! why /3.33??  bec
			  xstot = dxs   ! bec
c		We want probabilty of scattering, found from interaction length (dreact) found below.  
c         This depends on the total cross section.  The angular dependence is taken care of 
c		in choosing q from S(q).  To try and include some angular dependence here as well, 
c		would be redundant.  (see dsg0.f for same modification) -- bec
            endif
c     ............................................
c
c	the following lines are to turn off the cross section if you want to eliminate scattering effect
c		-- comment them out if you want scattering !!
c
c		if (j.eq.1) print *, " !Non-scattering run! (see target.f)"
c		xstot = 1.d-6
c
c	
		if (hegastar) then    ! if using He gas target
		  if (j.eq.1) write(6,*) ' !! He gas target run !! '
		  if (j.eq.1) write(3,*) ' !! He gas target run !! '
c             Get sigma(en0)  en0 in (meV) and find reaction length
c             !!This sigma(E) for air gap is fit of intepreted data from ENDF --
c               for n-4He  but not much actual data!!!!!
c
	        if (en0.lt.6.0) then
	           xstot = 0.0899/(en0/1000.0)**(0.474)    
	         else
	           xstot = 0.460/(en0/1000.0)**(0.173)    
	        endif	   
			denat = ndenhegas
			norm = 1.0
	    else if (airtar) then
	
c       Get sigma(en0)  en0 in (meV) and find reaction length
c      !!This sigma(E) for air gap is fit of intepreted data from ENDF --
c            for n-14N  , 14N cross section much bigger than 16O  
c            but not much actual data!!!!!
			 if (en0.lt.2.2) then
	           xstot = 1.0742/(en0/1000.0)**(0.47)   
	         else
	           xstot = 6.7799/(en0/1000.0)**(0.162)    
	        endif	   
			denat = ndenair
			norm = 1.0
		else
	      call heden(T,hedens) 
            denat = hedens*avog/he4mol
            call nden(T,norm)
		end if

c	write(6,*) en0, xstot, (denat*xstot*norm)

            dreact = -1.0*log(1-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.txin1.and.rdist(1).le.txin2) .and.
     &           (rdist(2).ge.tary1.and.rdist(2).le.tary2) .and.
     &           (rdist(3).ge.z0.and.rdist(3).le.(z0+tarz)) ) then  !scat in tar
	         xdis  = rdist(1) + (tarz-rfz)*dtan(thxp)		! x loc at end of tar
               ydis  = rdist(2) + (tarz-rfz)*dtan(thyp)		! y loc at end of tar
			pla = dsqrt(rfx**2+rfy**2+rfz**2)  ! path length before scat
			plb = dsqrt(((tarz-rfz)*dtan(thxp))**2+
     &              ((tarz-rfz)*dtan(thyp))**2+(tarz-rfz)**2)   ! path length after scat
		     if (ep.ne.en0) enflg=.true.
		     if (thsc.gt.0.d0) thflg=.true.
               if( (xdis.ge.txin1.and.xdis.le.txin2) .and.
     &              (ydis.ge.tary1.and.ydis.le.tary2) ) then
                  x = xdis
                  y = ydis
                  z = z0 + tarz
                  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 + tarz*dtan(thx0)
               ydis = y0 + tarz*dtan(thy0)
               if( (xdis.ge.txin1.and.xdis.le.txin2) .and.
     &              (ydis.ge.tary1.and.ydis.le.tary2) ) then
                  x = xdis
                  y = ydis
                  z = z0 + tarz
			pla = 0.0
			plb = 0.0
                  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
c         endif
      else
         x = -1.d6
         y = -1.d6
         z = -1.d6
         thx = -1.d6
         thy = -1.d6
         enp = -1.d6
      endif

      return
      end