#include <stdio.h>
#include <math.h>

#define PI (4.0*atan(1.0))
#define TWOPI (8.0*atan(1.0))

main() {
   int    iseed1, nline, count, iseed;
   double lenx, leny, lenz;
   double lambda, boxlength, tiny, mu0;
   double x, y, z, vx, vy, vz, mu, phi;
   double x0, y0, r, r1, v, vv, sumv;
   double vmax, vmin;
   double ran1();
   FILE *fp, *fopen();

   nline = 5000;
   iseed1 = -27659;
   lenx = PI/2.0;
   leny = PI/2.0;
   lenz = 10000.0;
/*
  LE ->  20keV   v = 0.0065291803950572;
  ME -> 100keV   v = 0.0145987576887395;
  HE ->   1MeV   v = 0.0461321482434411;
*/
   vmax = 0.0461321482434411;
   vmin = 0.0065291803950572;
   tiny = 1.0e-9;

   r = 0.2;
   x0 = PI/2.0;
   y0 = PI/2.0;
   mu0 = cos(40.0*PI/180.);
   sumv = 0.0;
/*
   lambda = 1.0;
   boxlength = lenz*lambda;
*/
   if((fp=fopen("data","w"))==NULL) {
      printf("Error to open 'data' file for writting.\n");
      exit(1);
   }
   for(count=0;count<nline;) {
      x = ran1(&iseed1)*lenx;
      y = ran1(&iseed1)*leny;
      z = ran1(&iseed1)*lenz;  
      z = 2.0; 

/*      mu = ran1(&iseed1)*2.0-1.0; */
      mu = 1.0-(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)*(x-x0)+(y-y0)*(y-y0))<= r*r) {
            count = count+1;
            fprintf(fp,"%5d %19.15lf %19.15lf %19.15lf %19.15lf %19.15lf %19.15lf\n",count,x,y,z,vx,vy,vz);
            sumv = sumv+sqrt(vx*vx+vy*vy+vz*vz);
/*            printf("%5d %15.8lf %20.18lf\n",count,mu,sumv); */
      }
   }
   fclose(fp);
   printf("count = %d\n",count);
   printf("v = %20.15lf = %20.15lf\n",sumv/(double)count,v);
}


#define ia 16807
#define im 2147483647
#define am (1.0/im)
#define iq 127773
#define ir 2836 
#define ntab 32 
#define ndiv (1+(im-1)/ntab)
#define eps 1.2e-7
#define rnmx (1.0-eps)

double ran1(idum)
int *idum;
{
   int j, k, imax(); 
   static int i=0, iy, iv[ntab+1];
   double  dmax(), dmin(), temp;

   if(i==0) {
      i+=1;   
      iy = 0; 
      for(j=1;j<=ntab;j++) {
         iv[j] = 0;
      }
   } 
   if((*idum<=0)||(iy<=0)) {
      *idum = imax(-*idum,1); 
      for(j=ntab+8;j>=1;j-=1) {
         k = *idum/iq;
         *idum = ia*(*idum-k*iq)-ir*k;
         if(*idum<0) *idum += im;
         if(j<=ntab) iv[j] = *idum;
      }
      iy = iv[1];
   }
   k = *idum/iq;
   *idum = ia*(*idum-k*iq)-ir*k;
   if(*idum<0) *idum = *idum+im;
   j = 1+iy/ndiv;
   iy = iv[j];
   iv[j] = *idum;
   temp = dmin(am*iy,rnmx);
   return temp;
}

int imax(x1,x2)
int x1, x2;
{
   if(x1>=x2) return x1;
   if(x1<x2) return x2;
}

double dmin(x1,x2)
double x1, x2;
{
   if(x1>=x2) return x2;
   if(x1<x2) return x1;
}


