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

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

main() {
   int    iseed, 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, mu, phi;
   double r, rr, vv, v;
   double ran1();
   FILE *fp, *fopen();

   nline = 1500;
   iseed = -27659;
   lenph = PI/2.;
   lenth = PI/2.;
   lenrr = 10000.;
/*
  LE ->  20keV   v = 0.0065291803950572;
  ME -> 100keV   v = 0.0145987576887395;
  HE ->   1MeV   v = 0.0461321482434411;
*/
   v = 0.0065291803950572;
   tiny = 1.0e-12;

   r = 0.3;
   sumv = 0.0;
   ph0 = PI/4.0;
   th0 = PI/4.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(&iseed)*lenph;
      th = ran1(&iseed)*lenth;
/*      rr = ran1(&iseed)*lenrr; */
      rr = 10.0;

      mu = ran1(&iseed);
      phi = ran1(&iseed)*TWOPI;
      vph = v*sqrt(1.0-mu*mu)*cos(phi);
      vth = v*sqrt(1.0-mu*mu)*sin(phi);
      vrr = v*mu;


      if((((th-th0)*(th-th0))+((ph-ph0)*(ph-ph0))) <= r*r) {
         muth = cos(th);
         muph = cos(ph);

/*
         vph = -vx*sqrt(1.-muph*muph) + vy*muph;
         vth = vx*muth*muph + vy*muth*sqrt(1.-muph*muph) - vz*sqrt(1.-muth*muth);
         vrr = vx*sqrt(1.-muth*muth)*muph + vy*sqrt(1.-muth*muth)*sqrt(1.-muph*muth) + vz*muth;
*/        
         phdot = vph/(rr*sqrt(1.0-muth*muth));
         thdot = vth/rr; 
         rrdot = vrr; 

         vx = phdot*rr*sqrt(1.0-muth*muth);
         vy = thdot*rr; 
         vz = rrdot; 

         vv = sqrt(vx*vx + vy*vy + vz*vz);
         vvr = sqrt(vph*vph + vth*vth + vrr*vrr);


         printf("vv = %10.6lf\n",vv);
         count = count+1;
         fprintf(fp,"%5d %25.15lf %25.15f %25.15lf %20.15lf %20.15lf %20.15lf\n",count,ph,th,rr,phdot,thdot,rrdot);
         sumv = sumv + sqrt(vph*vph + vth*vth + vrr*vrr);
      }
   }
   fclose(fp);
   printf("v = %10.6lf = %10.6lf = %10.6lf\n",vvr,vv,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;
}
