#include "peg_pic.h"


/***************************************************************/
/* Monte Carlo scheme for electron and ion-neutrals collisions */
/* First calculates maximum number of charged particles to be  */
/* collided per step. Then each time step chooses this number  */
/* of particles (randomly) and for each of these particles     */
/* checks probability for different type of collisions and if  */
/* probability is larger than the random number, then collides */

void mcc()
{
	static int init_flag=1;     
	static float ecol_extra[3], icol_extra[3][3];
	static float max_sigmav_e[3], max_sigmav_i[3][3];
	static float col_prob_e[3], col_prob_i[3][3], Nmax[3];

	int   N, i, ii, k, nnp, j, s, index, p, h;
	float random, dum, vel;
	float engy, rengy, del;
	float temp, vneutx, vneuty, vneutz;
	float vxtemp, vytemp, vztemp;
	float tempsigma1;
	double tempx;  

	if (init_flag) {
		k = nc2;
		s = nc3 + 1;

		for (j=0; j< Nnum; j++) {
			/* Max. neutral density */ 
			if (sourcD==3) {
				// No neutr. inside the source
				for (i=k; i<= s; i++) {
					Nneut[j][i] = 0.;
				}
			}

			Nmax[j] = Nneut[j][0];
			
			for (i=1; i< ng; i++) {
				Nmax[j] = max(Nmax[j], Nneut[j][i]);
			}
		}

		/* Calculating the null collision prob.  */

		for (j=0; j< Nnum; j++) {
			if(ecollisional) {
				/* Electrons */
				max_sigmav_e[j] = 0.0;

				for (i=0; i< 10*NEMAX; i++)	{
					engy = 0.1*i;
					max_sigmav_e[j]=max(max_sigmav_e[j],sqrt(2.0*1.6022e-19*engy/m[ecolsp])*(sigma1(j, engy)+sigma2(j, engy, 0)+sigma3(j, engy)));
				}
				
				col_prob_e[j] = 1 -exp(-max_sigmav_e[j]*Nmax[j]*dt); 
				ecol_extra[j]= 0.5;
			}
		}  /* end of for (j=0; j< Nnum; j++) loop */

		for (j=0; j< Nnum; j++) {
			/* Ions */
			for (k=0; k<Nnum; k++) {
				if(icollisional[j][k]) {
					max_sigmav_i[j][k] = 0.0;
					
					for (i=0; i<10*NEMAX; i++) {
						engy = 0.1*i;
						max_sigmav_i[j][k]=max(max_sigmav_i[j][k],sqrt(2.0*1.6022e-19*engy/m[icolsp[j]]/nimiu[j][k])*(sigma4(j, k, engy)+sigma5(j, k, engy)));
					}
					
					col_prob_i[j][k] = 1 -exp(-max_sigmav_i[j][k]*Nmax[k]*dt); 
					icol_extra[j][k]= 0.5;
					// printf(" Emax= %g, sigm_max= %g, nn= %g \n", engy, max_sigmav_i[j][k], Nmax[k]);
				}
			}  /* end of for (k=0; k< Nnum; k++) loop */
		}  /* end of for (j=0; j< nsp; j++) loop */

		init_flag = 0;
	}

	/* Step for calculation of temperature used in effective ionization routine */

	if ((!slow) && (efi)) {
		if (tstep%efi) {
			efiFL= 0;
		} else {
			efiFL= 1;
		}
	}

	/*****************************************************************/
	/***************** Collisions ************************************/

	for (ii=0; ii<Nnum; ii++) {
		/********** Electron collisions  ********************/
		if(ecollisional) {
			ecol_extra[ii] += np[ecolsp]*col_prob_e[ii];
			N = ecol_extra[ii];
			ecol_extra[ii] -= N;

			nnp = np[ecolsp];
			
			for(j=0; j< N; j++) {
				index= nnp*frand();

				if (x[ecolsp][index] < xnc) {
					nnp--;

					if (i_t[ecolsp]==2) {
						c_omt[ecolsp][index] = c_omt[ecolsp][nnp]; 
						c_omt[ecolsp][nnp] =0.;
					}

					tempx = x[ecolsp][nnp];              
					x[ecolsp][nnp] = x[ecolsp][index];
					x[ecolsp][index] = tempx;

					temp = vx[ecolsp][nnp];
					vx[ecolsp][nnp] = vx[ecolsp][index];
					vx[ecolsp][index] = temp;

					temp = vy[ecolsp][nnp];
					vy[ecolsp][nnp] = vy[ecolsp][index];
					vy[ecolsp][index] = temp;

					temp = vz[ecolsp][nnp];
					vz[ecolsp][nnp] = vz[ecolsp][index];
					vz[ecolsp][index] = temp;
				} else {  
					j--;
				}
			}

			for(j=nnp; j<nnp+N; j++) {
				/* determine if a collision occurs */
				s = x[ecolsp][j];
				dum = Ep(ecolsp, j);
				engy= Escale[ecolsp]*dum;
				vel = sqrt(dum)+1e-12;
				// sigma_total = Nmax[ii]*max_sigmav_e[ii]/vel/vscale; 
				tempsigma1 = Nneut[ii][s]/(Nmax[ii]*max_sigmav_e[ii]/vel/vscale);
				random= frand();
				// printf("en= %f, sigma= %e\n", engy, sigma2(ii, engy, 0));

				/*********************************************************************/
				/* determine the type of collision and calculate new velocities, etc */
				/*******************************/
				/* if the collision is elastic */

				if ((random -= tempsigma1*sigma1(ii, engy))<= 0.0) {
					/* first normalize vel */
					vx[ecolsp][j] /= vel;
					vy[ecolsp][j] /= vel;
					vz[ecolsp][j] /= vel;

					/* scatter the electron */
					newvel(engy,vel,&vx[ecolsp][j],&vy[ecolsp][j],&vz[ecolsp][j],ii,1,sctane[0][ii]);
				} else if ((random -= tempsigma1*sigma2(ii, engy, 0)) <= 0.0) {
					/* if the collision is excitation */
					if (slow) {
						del= x[ecolsp][j] - s;
						extrate[ii][s]  += (!s) ? 2*(1-del) : 1-del;
						extrate[ii][s+1]+= (s== nc-1) ? 2*del : del;
						if (s<=nc05) {
							Halhs_temp[ii] +=1.;
						} else {
							Harhs_temp[ii] +=1.;
						}
					} else if (Ef_flag) {
						if ( x[ecolsp][j]<=nc05) {
							Halhs_temp[ii] +=1.;
						} else {
							Harhs_temp[ii] +=1.;
						}
					}

					/* Multy-step ionization collisions with neutrals of the type 1 */

					if ((ii==0) && (efi)) {
						if (efsigma(Tx[ecolsp][s],srho[ecolsp][s]) > frand()) {
							goto EFF1;
						}
					}

					/* first normalize vel */
					vx[ecolsp][j] /= vel;
					vy[ecolsp][j] /= vel;
					vz[ecolsp][j] /= vel;

					if (Nexel[ii]) {
						for (k=1; k<=Nexel[ii]; k++) {
							if ((random -= tempsigma1*sigma2(ii, engy, k))<=0.0) {  
								break;
							}
						}
						
						engy -= extengy0[k-1][ii];
					} else { 
						engy -= extengy0[0][ii];
					}

					vel   = sqrt(fabs(engy)/Escale[ecolsp]);

					/* scatter the electron */
					newvel(engy,vel,&vx[ecolsp][j],&vy[ecolsp][j],&vz[ecolsp][j],ii,0,sctane[1][ii]);
				} else if((random - tempsigma1*sigma3(ii, engy)) <= 0.0) {
					/***************************************************************************/
					/* if the collision is ionization, add the released electron and ion first */
					EFF1:
					/* first normalize vel */
					vx[ecolsp][j] /= vel;
					vy[ecolsp][j] /= vel;
					vz[ecolsp][j] /= vel;

					/* subtract the ion. energy and partition the remaining energy */
					Escat(&engy, &rengy, ii);

					/* scatter the created electron */

					vel = sqrt(fabs(rengy)/Escale[ecolsp]);
					k = np[ecolsp];
					vx[ecolsp][k] = vx[ecolsp][j];
					vy[ecolsp][k] = vy[ecolsp][j];
					vz[ecolsp][k] = vz[ecolsp][j];
					newvel(rengy,vel,&vx[ecolsp][k],&vy[ecolsp][k],&vz[ecolsp][k],ii,0,sctane[2][ii]);
					x[ecolsp][k] = x[ecolsp][j];

					if (i_t[ecolsp]==2) {
						c_omt[ecolsp][k]=0;
					}

					/* assign velocities to the created ion */

					k = np[icolsp[ii]];
					fmaxwell(&vx[icolsp[ii]][k],&vy[icolsp[ii]][k],&vz[icolsp[ii]][k],Tneut[ii][s]); 
					x[icolsp[ii]][k] = x[ecolsp][j];
					
					if (i_t[icolsp[ii]]==2) {
						c_omt[icolsp[ii]][k]=0; 
					}

					/* finally scatter the incident electron */

					vel = sqrt(fabs(engy)/Escale[ecolsp]);
					newvel(engy,vel,&vx[ecolsp][j],&vy[ecolsp][j],&vz[ecolsp][j],ii,0,sctane[2][ii]);

					np[icolsp[ii]]++;
					np[ecolsp]++;

					if (slow) {
						del= x[icolsp[ii]][k] - s;
						ionrate[ii][s]  += (!s) ? 2*(1-del) : 1-del;
						ionrate[ii][s+1]+= (s== nc-1) ? 2*del : del;
					}
				}
			}
		}        /* end of if(ecollisional) loop */


		/******* Ion collisions  ************/
		/* "ii" now is the ion, "k" is the neutral  */

		i= icolsp[ii];

		for (k=0; k<Nnum; k++) {
			if(icollisional[ii][k]) {
				icol_extra[ii][k] +=  np[i]*col_prob_i[ii][k];
				N = icol_extra[ii][k];
				icol_extra[ii][k] -= N;

				nnp = np[i];
				
				for(j=0; j< N; j++) {
					index= nnp*frand();
					
					if (x[i][index] < xnc) {
						nnp--;

						if (i_t[i]==2) {
							c_omt[i][index] = c_omt[i][nnp];
							c_omt[i][nnp] =0.;
						}

						tempx = x[i][nnp];     
						x[i][nnp] = x[i][index];
						x[i][index] = tempx;

						temp = vx[i][nnp];
						vx[i][nnp] = vx[i][index];
						vx[i][index] = temp;

						temp = vy[i][nnp];
						vy[i][nnp] = vy[i][index];
						vy[i][index] = temp;

						temp = vz[i][nnp];
						vz[i][nnp] = vz[i][index];
						vz[i][index] = temp;
					} else {
						j--;
					}
				}

				for(j=nnp; j<nnp+N; j++) {
					s = x[i][j];              
					fmaxwell(&vneutx, &vneuty, &vneutz, Tneut[k][s]);  /* Neutral velocity */

					vxtemp = vx[i][j] - vneutx;                        /* Relative velocity */
					vytemp = vy[i][j] - vneuty;
					vztemp = vz[i][j] - vneutz;

					dum = (vxtemp*vxtemp + vytemp*vytemp + vztemp*vztemp);
					engy= nimiu[ii][k]*Escale[i]*dum;
					vel = sqrt(dum)+1.e-14;   
					//  sigma_total = Nmax[k]*max_sigmav_i[ii][k]/vel/vscale;      
					tempsigma1 = Nneut[k][s]/(Nmax[k]*max_sigmav_i[ii][k]/vel/vscale);

					random= frand();

					/*********************************************************************/
					/* determine the type of collision and calculate new velocities, etc */
					/*******************************/
					/* if the collision is elastic */

					if ((random -= tempsigma1*sigma4(ii, k, engy)) <= 0.0) {
						vxtemp /= vel;
						vytemp /= vel;
						vztemp /= vel;

						newvel(engy, vel, &vxtemp, &vytemp, &vztemp, k, 0, sctang[ii][k]);
						/*
	vel=(vxtemp*vxtemp + vytemp*vytemp + vztemp*vztemp);
	printf("E0= %e ,  E1= %e, vx= %e, vy= %e, vz= %e\n", dum, vel, vxtemp, vytemp, vztemp);
						 */
						vx[i][j] = vx[i][j] + nimiu[ii][k]*(vxtemp+vneutx-vx[i][j]);
						vy[i][j] = vy[i][j] + nimiu[ii][k]*(vytemp+vneuty-vy[i][j]);
						vz[i][j] = vz[i][j] + nimiu[ii][k]*(vztemp+vneutz-vz[i][j]);
					} else if ((random - tempsigma1*sigma5(ii, k, engy)) <= 0.0) {
						/***************************************/
						/* if the collision is charge exchange */
						if (slow) {
							del= x[i][j] - s;
							chrgxrate[ii][s]  += (!s) ? 2*(1-del) : 1-del;
							chrgxrate[ii][s+1]+= (s== nc) ? 2*del : del;
						}

						if (ii==k) {
							/* symmetric collisions, i.e. HH, DD, ... */
							vx[i][j] = vneutx;
							vy[i][j] = vneuty;
							vz[i][j] = vneutz;
						} else {
							/* asymmetric collisions, i.e. HD, DT, ... */
							h=icolsp[k];     /* new ion */
							p = np[h];
							np[h]++;
							x[h][p] = x[i][j];
							vx[h][p] = vneutx;
							vy[h][p] = vneuty;
							vz[h][p] = vneutz;

							if (i_t[h]==2) {
								c_omt[h][p]=0;
							}

							p=nnp+N-1;       /* new neutral */
							vx[i][j] = vx[i][p];
							vy[i][j] = vy[i][p];
							vz[i][j] = vz[i][p];
							
							if (i_t[i]==2) {
								c_omt[i][j]=0;
							}
							
							np[i]--;
							N--;
							j--;
						}
					} /* end of charge-exchange collisions */
				} /* end of for(j=nnp; j<nnp+N; j++) loop */
			} /* end of if(icollisional[ii][k]) loop */
		} /* end of for (k=0; k<Nnum; k++) loop */
	} /* end of for (ii=0; ii<Nnum; ii++) loop */
} /* End of mcc() */

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