program sigma-q


	implicit none

      include 'const.inc'


	integer loop
	double precision ep, wp, kp, k1, thp,php,qlast
      integer i, j, nevts,counter,lim,Event,cin,cot,trig  !bec
	integer Nenscat1,Ntar1,Nenscat2,Ntar2,Ntar12,NEdet,Nthscat1,Nthscat2
	integer encnt, Nsdet,Nthdet
      parameter(nevts=3e7)
      real xlam, HRNDM1, vect(5), lm(3), const, xinit, yinit, plr,xtest
      character*28 outfname
      double precision th0, ph0, thindex
      double precision en0, x0, y0, z0, thx0, thy0
      double precision en1, x1, y1, z1, thx1, thy1
      double precision en2, x2, y2, z2, thx2, thy2
      double precision x3, y3, z3, thx3, thy3
      double precision x4, y4, z4, thx4, thy4
      double precision x5, y5, z5, thx5, thy5
      double precision x6, y6, z6, thx6, thy6
      double precision x7, y7, z7, thx7, thy7
      double precision x8, y8, z8, thx8, thy8
      double precision x9, y9, z9, thx9, thy9
      double precision x10, y10, z10, thx10, thy10
      double precision time1, gl, gin, got
	logical tar2flg,enflg1,enflg2,thflg1,thflg2
      integer ranno
      integer seedc
      integer n1, n2, n3
C      parameter(ranno=1)
      double precision T 
      parameter	(gl = 111.0,     ! Guide length (cm)
     +          gin = 0.001095,    ! Input critical index (rad/Ang)
     +           T  = 2.10)     ! helium temp. (K)
c     +           T  = 4.25)     ! helium temp. (K)
c     +     got = 0.001)        ! Output critical index (rad/Ang)

      double precision r2d, lambda, sig, intlen, thscat, sigsom, sigsq
      double precision pl1_ns, th1_ns, pl2_ns, th2_ns, pl_ns, th_ns
      double precision scd1, scd2, pl1a_sc, pl1b_sc, th1a_sc, th1b_sc
      double precision pl2a_sc, pl2b_sc, th2a_sc, th2b_sc
      double precision pl1_sc, th1_sc, pl2_sc, th2_sc, pl3, th3, phi
	double precision pltot,thtot,pl_del
      CHARACTER day*2,month*2,hour*2,minute*2, filename*20
c-----------------------------------
      integer qmin, qmax
      double precision k0, qdummy, q
c-----------------------------------
c
c     -------------------------
      integer max_paw, iquest
      integer istat, icycle, datetime(8)
      parameter(max_paw = 1000000) ! PAW common block size
      real quest, p
	character date*12, time*12, zone*12
      common/pawc/p(max_paw)


      do 101 i = 1, nevts
c
         Event = 1



123        q=HRNDM1(301)
  
         call hf1(305,real(q),1.)

            if(q.le.0.562803) then

      rc  = (-4.7092+7.1047*T-1.9944*T**2+0.21769*T**3)/sqrt(3.)
*
      call dsg0(T,ds0dw)

      call heden(T,hedens)
      denat = hedens*avog/he4mol

       xstot = 4.*pi*ds0dw/(1.+(rc*q)**2)/denat  

            else
               call sommers(T, en0, dxs)
			  xstot = dxs   ! bec
            endif