/*
   7 June 2006 : gen_data_scat_par.c
      change name to gen100kev.c
*/
#include <stdio.h>
#include <math.h>

#define TINY 1.0e-4

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

main() {
   int    iseed1, iseed2, nline, count;
   double lenph, lenth, lenrr;
   double lambda, boxlength, tiny;
   double rrdot, thdot, phdot;
   double vrr, vth, vph, vvr;
   double muth, muph, sumv;
   double ph, th, th0, ph0;
   double vx, vy, vz, mu0, mu, phi;
   double r, rr, vv, v, mu1;
   double ran1();
   FILE *fp, *fopen();

   nline = 5000;
   iseed1 = -27659;
   iseed2 = -23439;
   lenph = 90*PI/180.;
   lenth = 90*PI/180.;
   lenrr = 10000.;
   v = 0.0461321482434411;
   tiny = 1.0e-9;

   r = 10.*PI/180.0;
   sumv = 0.0;
   ph0 = lenph/2.0;
   th0 = lenth/2.0;

   mu0 = cos(60.0*PI/180.0);
/*
   mu1 = cos(30.0*PI/180.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;) {
      ph = ran1(&iseed1)*lenph;
      th = ran1(&iseed1)*lenth;
      rr = ran1(&iseed1)*lenrr; 
      rr = 2.5;

/*      mu = ran1(&iseed2)*(1.-mu0)+mu0; */
      mu = ran1(&iseed2)*(1.-mu0)+mu0;
      phi = ran1(&iseed2)*TWOPI;


      if((((th-th0)*(th-th0))+((ph-ph0)*(ph-ph0))) <= r*r) { 
         vph = v*sqrt(1.0-mu*mu)*cos(phi);
         vth = v*sqrt(1.0-mu*mu)*sin(phi);
         vrr = v*mu;

         vvr = sqrt(vph*vph + vth*vth + vrr*vrr);

         count = count+1;
         fprintf(fp,"%5d %20.15lf %20.15f %20.15lf %20.15lf %20.15lf %20.15lf\n",count,ph,th,rr,vph,vth,vrr);
/*
         fprintf(fp,"%5d %15.10lf %15.10f %15.10lf %15.10lf %15.10lf %15.10lf\n",count,ph,th,rr,ph,th,rr);
*/
         sumv = sumv + sqrt(vph*vph + vth*vth + vrr*vrr);
      } 
   }
   fclose(fp);
   printf("v = %10.6lf = %10.6lf\n",vvr,sumv/count);
}


#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;
}
