subroutine wq(val,eval,en,k)

      implicit none

      include 'const.inc'
 
      integer i,n
      parameter(n=37)
      double precision q(n), w(n), en, k, step, tmp
      double precision val,eval,slop,y_int

      data q/ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
     &     1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0,
     &     2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0,
     &     3.1, 3.2, 3.3, 3.4, 3.5, 3.6/             ! (1/A)
      data w/  0.000000, 0.235735, 0.484567, 0.739947, 0.969134,
     &     1.198321, 1.407863, 1.538827, 1.656695, 1.722177, 1.774562,
     &     1.807303, 1.800755, 1.768014, 1.695984, 1.597761, 1.466797,
     &     1.342381, 1.211417, 1.139387, 1.172128, 1.309640, 1.525730,
     &     1.774562, 2.029942, 2.154358, 2.226388, 2.265677, 2.291870,
     &     2.318063, 2.337707, 2.357352, 2.376997, 2.396641, 2.403189,
     &     2.409738, 2.416286/ ! 1/psec

c	elastic scattering -- no energy change
c	eval = 0.d0
c	return


      do i=1,n-1
         if(val.lt.q(1).or.val.gt.q(n)) then
            write(*,*) i, val,q(1),q(n),
     +           ' The value is out of range ... wq code'
            if(val.lt.q(1)) eval=w(1)
            if(val.gt.q(n)) eval = w(n)
         elseif(val.eq.q(i)) then
c		write(*,*) ' val eq qi '
            eval = w(i)
         elseif(val.eq.q(i+1)) then
c		write(*,*) ' val eq qi+1 '
            eval = w(i+1)
         elseif(val.gt.q(i).and.val.lt.q(i+1)) then
            slop = (w(i+1)-w(i))/(q(i+1)-q(i))
            y_int = w(i) - slop*q(i)
            eval = slop*val+y_int
c	write(*,*) val,slop,y_int,eval
         endif
      enddo

c
c	The following finds an omega value when q<0.1 that will conserve both momentum
c         and energy.  In the original code all neutrons with E<1meV were let pass
c         through the target with no interaction to avoid potential problems when
c         calculating theta_scatter.  Rather than doing that, main_cwn.f now looks for another
c         q value that will convserve p and E.  In addition, the following takes the omega 
c		value given from the interpolation of the table above, checks if theta_scat can 
c		be calculated, and if not reduces omega by 1/20th of its value and checks again.  -- bec
c


	if (val.lt.0.0D0) then
		write(*,*) ' wq:  q<0 ', q
		goto 102
	endif
	if (val.lt.0.1d0) then
		i = 0
		step = eval/20
 100		i=i+1
c		if (i.gt.1) print *, ' wq loop ', i
		if (eval.lt.0.d0) then
			eval = 0.d0
			print *, ' setting wp=0'
			goto 102
		endif
		if ((en-h_par*eval).lt.0.d0) goto 101
		tmp = 2.d0*m_n*(en-h_par*eval)/h_par**2-
     &		2.d0*k*dsqrt(2.d0*m_n*(en-h_par*eval))/h_par
		if (tmp.lt.val**2-k**2) goto 102
 101		eval = eval - step
		goto 100
	end if

 102      return
      end