#include "peg_pic.h"


void ccvel(float *, float *, float *, float *, float *,  float *, float, float, float);

/***************************************************************/
/* Monte Carlo scheme for binary Coulomb collisions ddd  */ 

void bincc() {
	int          isp, s, s_1, nnp, ncl, i, im, ii, ncl1;
	float        ppcc, m_0l, m_1l;
	static int   init_flag=1, *ldav, *Pdav, *col;
	static float *ppl;

	if (init_flag) {
		Pdav= (int *)malloc(MAXLEN*sizeof(int));
		col = (int *)malloc(nc*sizeof(int));
		ppl = (float *)malloc(nc*sizeof(float));
		
		if (ccc>1) {
			ldav= (int *)malloc(nc*sizeof(int));
		}

		init_flag = 0;
	}

	/*  collision with the particles of same species  */

	for (isp=0; isp<ccc; isp++)	{
		if (ccc>1) {
			for (s=0; s<nc; s++) {
				for (i=1; i<=MAXD; i++) {
					Ldav[isp][s][i] = 0;
				}
			}
			
			for (ncl=0; ncl<nc; ncl++) {
				ldav[ncl]=0;
			}
		}
		
		for (ncl=0; ncl<nc; ncl++) {
			col[ncl]=0;
		}
		
		for (i=1; i<=np[isp]; i++) {
			Pdav[i]=0;
		}
		
		for (ii=0; ii<nc; ii++) { 
			ncl1 = ii+1;
			ppl[ii] = pp[isp]*(srho[isp][ii]+srho[isp][ncl1]); 
		}

		for (nnp=np[isp]-1; nnp>=0; nnp--) { 
			s = 1+(nnp + 1)*frand(); /* random particle */ 


			if (Pdav[s]) {
				im = s;
				s = Pdav[s];

				i = nnp+1;                 /* if "s" will be chosen once more, then the real particle will be "nnp" or "Pdav_(nnp+1)-1"  */

				if (Pdav[i]) {
					Pdav[im] = Pdav[i];
				} else {
					Pdav[im] = i;
				}
			} else {
				i = nnp+1;                 /* if "s" will be chosen once more, then
                         the real particle will be "nnp" or "Pdav_(nnp+1)-1"  */
				if (Pdav[i]) {
					Pdav[s] = Pdav[i];
				} else {
					Pdav[s] = i;
				}
			}

			s_1 = s - 1;

			if (x[isp][s_1] < xnc) {                    /***   for i_t=2 stuff ***/
				ncl = x[isp][s_1];

				if (ccc > 1) {
					ldav[ncl]++;
					Ldav[isp][ncl][ldav[ncl]] = s;
				}

				/** if there is no saved particle in the given cell then   
                                    particle is saved with the address ncl **/
				if (!col[ncl]) {
					col[ncl] = s; 
				} else {   /*  in the same cell there is particle for binary collision  */
					ii = col[ncl] -1 ;
					col[ncl]=0;

					/****  collision between the "ii" and "s_1" particles  *****/
					ccvel(&vx[isp][s_1], &vy[isp][s_1], &vz[isp][s_1], &vx[isp][ii], &vy[isp][ii], &vz[isp][ii], 0.5, 0.5, ppl[ncl]);
					/*************** Cheking of energy conservation *********************
          chekk[isp] += (vx[isp][s_1]-vx[isp][ii]-0.5*dux)*dux +
                        (vy[isp][s_1]-vy[isp][ii]-0.5*duy)*duy +
                        (vz[isp][s_1]-vz[isp][ii]-0.5*duz)*duz;

      printf("delta_kin= %22.21f\n", chekk[0]);  
					 ********************************************************************/

				}     /* end of else  */
			}  /* end of for if (x>=0) */
		}   /* end of for (nnp...   */
	}    /* end of for (isp=0... */


	/***** Collisions between different species, max. 2 species (0--1)  ******/

	if (ccc > 1) {
		m_0l= m_0;
		m_1l= m_1;

		for (ncl=0; ncl < nc; ncl++) {
			ncl1 = ncl+1;
			ppcc = pp_01*max((srho[0][ncl]+srho[0][ncl1]), (srho[1][ncl]+srho[1][ncl1]));

			im = 1; 

			while ((Ldav[0][ncl][im]) && (Ldav[1][ncl][im])) {
				i = Ldav[0][ncl][im] - 1;
				s = Ldav[1][ncl][im] - 1;

				/* collision of the electron "i" and the ion "s"*/
				ccvel(&vx[1][s], &vy[1][s], &vz[1][s], &vx[0][i], &vy[0][i], &vz[0][i], m_0l, m_1l, ppcc);

				im++;
			}   /* end of while (...         */
		}    /* end of for (ncl=0   ...  */
	}     /*  end of if (ccc > 0)  */ 
}   /* end of bincc() */ 

/*****************************************************/ 

void ccvel(float *vx1, float *vy1, float *vz1, float *vx2, float *vy2, float *vz2, float m1, float m2, float coef) {
	float mag;
	double  u, ux, uy, uz, up, c_th, col_th, s_th, sc, ss, dux, duy, duz;

	ux = *vx1 - *vx2;
	uy = *vy1 - *vy2;
	uz = *vz1 - *vz2;
	up = sqrt(ux*ux + uy*uy);
	u = sqrt(ux*ux + uy*uy + uz*uz);
	mag=frand()*TWOPI;

	col_th = sqrt(coef*fabs(log(frand()))/(u*u*u + DBL_MIN));

	if ( col_th>1.570796327) {
		col_th = frand()*1.570796327;
	}
	
	c_th = cos(col_th)-1.;

	if (up)	{
		s_th = sin(col_th)/up;
		ss = s_th*sin(mag);
		sc = s_th*cos(mag);

		dux = ux*uz*sc - uy*u*ss + ux*c_th;
		duy = uy*uz*sc + ux*u*ss + uy*c_th;
		duz = -up*up*sc + uz*c_th;
	} else {
		s_th = sin(col_th)*uz;

		dux = s_th*cos(mag);
		duy = s_th*sin(mag);
		duz = uz*c_th;
	}

	*vx1 += m1*dux;
	*vy1 += m1*duy;
	*vz1 += m1*duz;

	*vx2 -= m2*dux;
	*vy2 -= m2*duy;
	*vz2 -= m2*duz;
}
/* END of ccvel() */

/*******************************************************/
