/*
   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,rr0;
   double ran1();
   FILE *fp, *fopen(),fp_in;

   rr0 = 50.0;
/*
   nline = 5000;
*/
   nline = 20000;   

   iseed1 = -27659;
   iseed2 = -23439;
   lenph = 90*PI/180.;
   lenth = 90*PI/180.;
   lenrr = 10000.;
/*
  LE ->  20keV   v = 0.0065291803950572;
  ME -> 100keV   v = 0.0145987576887395;
  HE ->   1MeV   v = 0.0461321482434411;
*/
/*  10Mev   voc = 0.144844008620218;  */
/* 100Mev   voc = 0.428195496560756;  */
/*   1Gev   voc = 0.875025655404499;  */
   
   v = 0.428195496560756;
   tiny = 1.0e-9;

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

   mu0 = cos(60.0*PI/180.0);
/*
   mu1 = cos(30.0*PI/180.0);
   lambda = 1.0;
   boxlength = lenz*lambda;
*/
   
/*
   input b_slab
*/   
   if ((fp_in = fopen("bx1d.dat","r")) == NULL ) {
      printf("File could not be opened\n");
      exit();
   } 
   else {   
      printf("input bxslab\n");
      for (i=1;i<=NZ;i++) fscanf(fp_in,"%lf\n",&bx1d[i]);
   }         
   fclose(fp_in);   
   if ((fp_in = fopen("by1d.dat","r")) == NULL ) {
      printf("File could not be opened\n");
      exit();
   } 
   else {   
      printf("input byslab\n");
      for (i=1;i<=NZ;i++) fscanf(fp_in,"%lf\n",&by1d[i]);
   }         
   fclose(fp_in);   
/*
     input b2d 
*/
   if ((fp_in = fopen("bx2d.dat","r")) == NULL ) {
      printf("File could not be opened\n");
      exit();
   } 
   else {   
      printf("input bx2D\n");
      for (k=1;k<=(NX*NY);k++) fscanf(fp_in,"%lf\n",&bx2d[k-(k/NX)*NX][k/NX+1]);
   }         
   fclose(fp_in);   
   if ((fp_in = fopen("by2d.dat","r")) == NULL ) {
      printf("File could not be opened\n");
      exit();
   } 
   else {   
      printf("input by2D\n");
      for (k=1;k<=(NX*NY);k++) fscanf(fp_in,"%lf\n",&by2d[k-(k/NX)*NX][k/NX+1]);
   }         
   fclose(fp_in);   
   
   /* check input slab */
   for (i=1;i<=NZ;i++) {
      printf("bxslab[%d]=%le\n ",i,bx1d[i]);
      printf("byslab[%d]=%le\n ",i,by1d[i]);
   
   }

   /* check input a2d */
    for(j=1;j<=NY;j++) {
      for(i=1;i<=NX;i++) { 
         printf("bx[%d][%d]=%le\n ", i,j,bx2d[i][j]); 
         printf("by[%d][%d]=%le\n ", i,j,by2d[i][j]); 
      }
   }   

   
   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;
*/
      ph = 2.*ran1(&iseed1)*lenph - lenph;
      th = 2.*ran1(&iseed1)*lenth;
      rr = ran1(&iseed1)*lenrr; 
/*
      rr = 2.5;
*/
      rr = 5.0;

/*      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;
*/
         
         zz = rr/delz - rr;   
         xx = (ph+x0)/delx;  /* x0 = lenx/2. */
         yy = (y0-th)/dely;  /* y0 = pi/2.  */
/*  
    interpolate bslab
*/
         k = z;
         v = z-(double) k;
         omv = 1.d0-v;
   	 k = k+1;
         kp = k+1;
         bxslab = bx1d[kp]*v + bx1d[k]*omv;
         byslab = bxslab; 
/*	 
   interpolate b2D
*/
    /* b_\phi_2d */   	 
         i = xx;
         j = yy;
         s = xx - (double) i;
         u = yy - (double) j;
         oms = 1.0-s;
         omu = 1.0-u;
	 
         i = i+1;
         j = j+1;
         ip = i+1;
         jp = j+1;

         if(i > NX)  i=i-NX;
         if(ip > NX) ip=ip-NX;
         if(j > NY)  j=j-NY;
         if(jp > NY) jp=jp-NY;

         a1 = bx2d[ip][jp];
         a2 = bx2d[i][jp];
         a3 = bx2d[ip][j];
         a4 = bx2d[i][j];
         bx2d2 = a1*s*u+a2*oms*u + a3*s*omu+a4*oms*omu;
     	 
    /* b_theta  */
         a1 = by2d[ip][jp];
         a2 = by2d[i][jp];
         a3 = by2d[ip][j];
         a4 = by2d[i][j];
         by2d2 = a1*s*u+a2*oms*u + a3*s*omu+a4*oms*omu;

    /* magnetic field at that position  */	  
         bph = bxslab*rr0*rr0/(rr*rr)+bx2d2d*rr0*rr0/(rr*rr);
         bth = byslab*rr0*rr0/(rr*rr)-by2d2d*rr0*rr0/(sin(th)*rr*rr);
         brr = b0*rr0*rr0/(rr*rr);

	 bb = sqrt(bph*bph + bth*bth + brr*brr);
		 
	 bphhat = bph/bb;
	 bthhat = bth/bb;
	 brrhat = br/bb;
	 bphth = sqrt(bphhat*bphhat + bthhat*bthhat);

     /* define phip and thetap */
	 thetap = atan2(bphth,brrhat);
         phip = atan2(bthhat,bphhat);  

     /* v||b */     
         vph = v*cos(phip)*sin(thetap);
         vth = v*sin(phip)*sin(thetap);
	 vrr = v*cos(thetap);
	 
         vvr = sqrt(vph*vph + vth*vth + vrr*vrr);
	 
     /* check if v||b  */
         cosvb = (vph*bph + vth*bth + vrr*brr)/(vrr*bb);	 
	 printf("cosvb = %lf\n",cosvb);

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