program main

      implicit none

      include 'const.inc'
      integer i, j, n, nevts, counter, ranno
      parameter(nevts=1e3)
      parameter(n=19)
      real lm, HRNDM1, vect(6), const, xinit, yinit, plr,xtest
      double precision x0, y0, z0, th0, ph0
      double precision x(n), y(n), z(n)
      double precision bx(n), by(n), bz(n), t1
      double precision lmin, lmax, thindex
      double precision xmin, xmax, ymin, ymax, zmin, zmax
      double precision throtx, throty, throtz, en
c
c     -------------------------
c     fill the ntuple 'spin.rz' //// ntuple variables ....
      integer max_paw, iquest
      integer istat, icycle
      parameter(max_paw = 10000000) ! PAW common block size
      real quest, p
      common/quest/iquest(100)
      common/pawc/p(max_paw)
      real lam, energy, xi, yi, zi, thi, phi
      real trx, try, trz, xf, yf, zf
      common /spinblk/ i,lam, energy,
     &     xi, yi, zi, thi, phi,
     &     trx, try, trz, xf, yf, zf
      call hlimit(max_paw)
      iquest(10) = 6500000
      call hropen(1,'Myntuple','picoil.rz','NQ',8190,istat)
      call hbnt(10,'spinnt',' ')
c--------------------------------
      call hbname(10,'spin',i,'i:I,lam:R,energy:R,xi:R,yi:R,zi:R,'//
     &     'thi:R, phi:R,trx:R, try:R, trz:R, xf:R, yf:R, zf:R')
c-----------------------------------------------------------------------
      call hrget(100,'lam.hbook','')
c-----------------------------------------------------------------------
c
ccc      t1=time()
ccc      ranno=mod(t1,1576879653)
ccc      print *,"asking for seed"
ccc      read(*,*) ranno;
c
c      call rluxgo(3,ranno,0,0)  ! RANLUX luxury level 3
      call rluxgo(3,11,0,0)  ! RANLUX luxury level 3
c-----------------------------------------------------------------------
c
c      stop
      lmin = 2.0   ; lmax = 15.0
      xmin = -4.5  ; xmax = 4.5
      ymin = -2.0  ; ymax = 2.0
      zmin = -9.0  ; zmax = 9.0
c
      do 100 i = 1, nevts
c`
         lm = HRNDM1(100)
         if(lm.lt.1.) goto 100
         en = (1/(2*m_n))*(2*pi*h_par/lm)**2
c
         call ranlux(vect,4)
         x0 = xmin + (xmax-xmin)*vect(1)
         y0 = ymin + (ymax-ymin)*vect(2)
         z0 = zmin
         thindex = 0.003*lm
         th0  = thindex*vect(3)
         ph0  = 2.00*pi*vect(4)

         throtx = 0.
         throty = 0.
         throtz = 0.

         polx0 = 0.0
         poly0 = 1.0
         polz0 = 0.0

         do j=1,n

            z(j) = z0+(j-1)

            if(j.eq.1) then
               x(j) = x0
               y(j) = y0
               polx(j) = polx0
               poly(j) = poly0
               polz(j) = polz0
            else
               x(j) = x(j-1)+z(j)*dsin(th0)*dcos(th0)
               y(j) = y(j-1)+z(j)*dsin(th0)*dsin(th0)
               polx(j) = polx(j-1)+
               poly(j) = poly(j-1)+
               polz(j) = polz(j-1)+
            endif
c
c            call mag_b(x(j),y(j),z(j),bx(j),48,1)
c            call mag_b(x(j),y(j),z(j),by(j),48,2)
            bx(j)=0.
            by(j)=0.
            call mag_b(x(j),y(j),z(j),bz(j),60,3)

            Polx=Polx+

            
            throtx=throtx+(gy*0.0001*m_n*lm/h_blanck)*(bz(j)-by(j))
            throty=throty+(gy*0.0001*m_n*lm/h_blanck)*(bx(j)-bz(j))
            throtz=throtz+(gy*0.0001*m_n*lm/h_blanck)*(by(j)-bx(j))

         enddo

c         write(11,*) throtx,throty,throtz
c
c     Fill the ntuple variables
c
         lam=lm
         energy = en
         xi = x0
         yi = y0
         zi = z0
         thi = th0
         phi = ph0
         xf = x(n)
         yf = y(n)
         zf = z(n)
         trx = throtx
         try = throty
         trz = throtz
         call hfnt(10)

 100  continue
c
c     -------------------------
      call hrout(10,icycle,' ')
      call hrend('Myntuple')
c     -------------------------
c
      return
      end