program omegatest


	implicit none
      include 'const.inc'



	double precision q, wp, en0, ep, k0, kp
	double precision xlam, HRNDM1, p



	integer max_paw, iquest,qmin, qmax
      integer i, loop, encnt, nevts
      parameter(max_paw = 1000000) ! PAW common block size
	parameter(nevts = 1e5)
      common/pawc/p(max_paw)

      call hlimit(max_paw)

      call hrget(100,'lam.hbook','')
c      call hrget(102,'xyspot.hbook','')
      call hrget(151,'gradx.hbook','')
      call hrget(152,'gradx.hbook','')
      call hrget(301,'sq.hbook','')


	call hbook1(305,'q sim(1/Ang)',100,0,4.0,0.)      ! bec test q choices
      call hbook1(306,'w sim(1/psec)',100,0,4.0,0.)      ! bec test q choices
      call hbook1(307,'Eo sim(meV)',3000,0.,40.0,0.)      ! bec test q choices
      call hbook1(308,'Ep sim(meV)',3000,0.,40.0,0.)      ! bec test q choices
      call hbook1(309,'ko sim(1/Ang)',100,0,4.0,0.)      ! bec test q choices
      call hbook1(310,'kp sim(1/Ang)',100,0,4.0,0.)      ! bec test q choices

c-----------------------------------------------------------------------


	do i = 1, nevts

         xlam = HRNDM1(100)
c
c     The energy calculated from the neutron's wavelength
c
         if(xlam.lt.1.) goto 100
         en0 = (1/(2*m_n))*(2*pi*h_par/xlam)**2
            call hf1(20,real(en0),1.)
	   if (en0.le.1.0) encnt =encnt+1
c
         k0 = sqrt(2*m_n*en0)/h_par
         qmin = 0
c         qmax = int(k0*25.)     ! for sq_1.hbook
         qmax = int(2.0*k0*25.)     ! bec  qmax = 2*k  from p-consv. 
         if((2.0*k0).gt.3.6) qmax=int(3.6*25.)
 		if (i.gt.1) call hdelet(302)
c	write(*,*) ' k0, e0 ',k0, en0
c	write(*,*) ' qmin, qmax ',qmin,qmax

         call hcopyr(301,302,"",qmin,qmax,0.,0.,"")

c
c         qdummy=HRNDM1(301)  ! bec 
c   HRNDM1 replaces spec with integrated spec!, so next pass through the loop
c   hcopyr is copying the integrated spec into 302, not the original spec!!
c   without the qdummy=HRNDM1(301) line, the simulated q-spec (305) looks
c   much more reasonable!  bec
		loop=0
 123        q=HRNDM1(302)  
c  		 q = 3.6/nevts*i 
c	write(*,*) '   q = ', q
		loop=loop+1
       	call wq(q,wp,en0,k0)  !  wq modified to prevent E,p consv. probs
		ep  = en0- (h_par*wp)
		if (ep.le.0.) goto 123

c		if (ep.le.0.) then
c		  ep = 1.
c		  else
		  kp=sqrt(2*m_n*ep)/h_par
c		endif

		if (dabs((k0**2+kp**2-q**2)/(2.*k0*kp)).le.1.d0) goto 124
      goto 123
 124	continue


         call hf1(305,real(q),1.)
         call hf1(306,real(wp),1.)
         call hf1(307,real(en0),1.)
         call hf1(308,real(ep),1.)
         call hf1(309,real(k0),1.)
         call hf1(310,real(kp),1.)



 100	continue

	enddo 

      call hrput(0,'omegatest.o','N')
	print *, 'omegatest.o'

	end
***