      program gen_data_scatter

      integer iseed1, nline, count, iseed
      real*8 lenx, leny, lenz
      real*8 lambda, boxlength, tiny, mu0
      real*8 x, y, z, vx, vy, vz, mu, phi
      real*8 x0, y0, r, r1, v, vv, sumv
      real*8 vmax, vmin

      real*8 PI
      parameter (PI = 3.14159265358979323846d0)

      real*8 TWOPI
      parameter (TWOPI = 6.28318530717958647692d0)

c  LE ->  20keV   v = 0.0065291803950572;
c  ME -> 100keV   v = 0.0145987576887395;
c  HE ->   1MeV   v = 0.0461321482434411;
c   vmax = 0.0461321482434411;
c   vmin = 0.0065291803950572;

      nline = 5000
      iseed1 = -27659
      lenx = PI/2.0
      leny = PI/2.0
      lenz = 10000.0
 
      tiny = 1.0e-9
      r =  0.2
      x0 = PI/4.0
      y0 = PI/4.0
      mu0 = cos(40.0*PI/180.)
      sumv = 0.0
      count = 0  
     
      open(10,file='data')
      do while (count .lt. nline) 
         x = ran1(iseed1)*lenx
         y = ran1(iseed1)*leny
         z = ran1(iseed1)*lenz
         z = 2.0

         mu = 1.0d0-(ran1(iseed1)*(1.0-mu0))
         phi = ran1(iseed1)*TWOPI

   
         v = ran1(iseed1)*(vmax-vmin)+vmin
         v = 0.0065291803950572

         vx = v*cos(phi)*sqrt(1.0-mu*mu)
         vy = v*sin(phi)*sqrt(1.0-mu*mu)
         vz = v*mu
         vv = sqrt(vx*vx + vy*vy + vz*vz)

         if ((x-x0)**2+(y-y0)**2 .le. r*r) then
            count = count+1
            write(10,100)count,x,y,z,vx,vy,vz
         endif
      enddo
 100  format(I5,6F12.5)
      close(10)
      end
      
      FUNCTION ran1(idum)
      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      double precision ran1,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     &       NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)

      INTEGER j,k,iv(NTAB),iy/0/
      SAVE iv,iy
      DATA iv /NTAB*0/, iy/0/

      if (idum.le.0.or.iy.eq.0) then
          idum=max(-idum,1)
          do j=NTAB+8,1,-1
             k=idum/IQ
             idum=IA*(idum-k*IQ)-IR*k
             if (idum.lt.0) idum=idum+IM
             if (j.le.NTAB) iv(j)=idum
          enddo
          iy=iv(1)
       endif
       k=idum/IQ
       idum=IA*(idum-k*IQ)-IR*k
       if(idum.lt.0) idum=idum+IM
       j=1+iy/NDIV
       iy=iv(j)
       iv(j)=idum
       ran1=min(AM*iy,RNMX)
       return
       END
              
      
      

