#include "pdefs.h"
#include <malloc.h>

/****************************************************************/
/* Routine to adjust (initialize, re-pack and inject) particles */
/*	to the desired boundary conditions                           */

adjust()
{
	static float extra[NSMAX][2], iv[NSMAX][2], eVscale, iVscale, ionVscale;
	int i, ii, isp, k, nnp, j, secountl, secountr, ien, s;
	static int init_flag=1;
	float	sig1, sig2, stot, random, prob, dum, vel, engy, rengy, eengy, del;
	float	phi1, cosphi, sinphi, coschi, sinchi, tempvr, tempvth, tempvph;
	float vnext(), del_t, vmag, theta;
	int nreflux[NSMAX], ninjected;
	
	/* INITIALIZE array for computing positions of injected particles */
	if (init_flag)
	{
		/* "enter" is now the no. of particles injected each "dt" */
		for (isp=0; isp<nsp; isp++)
		{
			if (fabs(enter[isp][0]) > 0.) extra[isp][0] = 0.5123123;
			else extra[isp][0] = 0.;
			if (fabs(enter[isp][1]) > 0.) extra[isp][1] = 0.5123123;
			else extra[isp][1] = 0.;
		}
		if(icollisional)
      {
			achrgx *= gden*dr;
			bchrgx *= gden*dt*sqrt(2*1.602E-19/m[icolsp]);
			ascat  *= gden*dr;
			bscat  *= gden*dt*sqrt(2*1.602E-19/m[icolsp]);
		}
		init_flag = 0;
	}

	for (isp=0; isp <nsp; isp++) jwall[isp] = 0.0;

	if(secondary) secountl = secountr = 0;
	
	for (isp=0; isp<nsp; isp++)
	{
		nreflux[isp]= 0;	/* nreflux[isp] equals zero */
		if (np[isp] > 0)
		{
			nnp = np[isp] - 1;
			i = 0;

			/* eliminate "outsiders", allow for secondary electron emission, 
				and if it left thru LH plate, 
				ADD charge there to sigma (plate surface density). */
			do
			{ 
				if (r[isp][i] <= r0) 
				{
					dum = (vr[isp][i]*vr[isp][i] + vth[isp][i]*vth[isp][i]
							+vph[isp][i]*vph[isp][i] - emin[isp])/de[isp];
					s = dum;
					if (s<nbin[isp]-1 && dum>=0)
					{
						dum -= s;
						
						fe[isp][s]  += (!s) ? 2*(1-dum) : 1-dum;
						fe[isp][s+1]+= (s==nbin[isp]-1) ? 2*dum : dum;
					}

         		dum = -atan(sqrt(vth[isp][i]*vth[isp][i]
											+vph[isp][i]*vph[isp][i])
											/vr[isp][i])/dtheta[isp];
					s = dum;
					dum -= s;
					ftheta[isp][s]  += (!s) ? 2*(1-dum) : 1-dum;
					ftheta[isp][s+1]+= (s==nbin[isp]-1) ? 2*dum : dum;

					r[isp][i] = r[isp][nnp];
					vr[isp][i] = vr[isp][nnp];
					vth[isp][i] = vth[isp][nnp];
					vph[isp][i] = vph[isp][nnp];
					nnp--;
					
					jwall[isp] += q[isp]/area;
					
					if (secondary)
					{
						if (frand() < seec[isp])
						{
						    secountl++;
							 jwall[secsp] -= q[secsp]/area;
					   }
					}
				}

				else if (r[isp][i] >= r1)
				{
					r[isp][i] = r[isp][nnp];
					vr[isp][i] = vr[isp][nnp];
					vth[isp][i] = vth[isp][nnp];
					vph[isp][i] = vph[isp][nnp];
					nnp--;
					if(reflux)	nreflux[isp]++;
					else if (secondary)
						if (frand() < seec[isp])	secountr++;
				}
				else i++;
			} while (i < nnp);
			
			np[isp] = nnp +1;
		}	/* end of "if(np[isp]>0...)" loop */

	}	/* end of "for( isp..)"  loop */

	/* INJECT new particles at walls, one species at a time */
	
	for (isp=0; isp<nsp; isp++)
	{
		if (inject[isp] || nreflux[isp])
		{
			for (k=0; k<2; k++)
			{
				extra[isp][k] += enter[isp][k];
				if (reflux && k == 1) extra[isp][k] += nreflux[isp];
				/* "extra" = the fraction of "dt" when emitted */
				ninjected = extra[isp][k];

				while (extra[isp][k] >= 1.)
				{
					del_t= ((int)extra[isp][k])/(ninjected +1.);
					ii = np[isp];
					np[isp]++;
				
					if (ii > MAXLEN) 	 /* Move array boundaries here */
					{
						printf("ADJUST: too many particles, species %d",isp);
						exit(1);
					}

					/* Choose Vr in bit-reversed order */
					vr[isp][ii] = vnext(k,isp);
					vmag = vtpt[isp]*sqrt(2*fabs(log(frand()+DBL_MIN)));
					theta = TWOPI*frand();
					vth[isp][ii] = vmag*sin(theta);
					vph[isp][ii] = vmag*cos(theta);
					
					/* Adjust Vr for effect of electric field */
					vr[isp][ii] += qm[isp]*(del_t-0.5)*e[k*nc]*dttx;
					r[isp][ii] = r0 +k*nc + del_t*vr[isp][ii];
					if(!k)	jwall[isp] -= q[isp]/area;
					extra[isp][k] -= 1;
				}	/* end "while extra..)" loop */

			}	/* end "for k= 0,1" loop */

		}	/* end "if( inject..)"  loop */

	}  /* end "for( isp...)"  loop */
	
	if(ecollisional)
	{
		nnp = np[ecolsp];
	  	for(j=0; j<nnp; j++)
	  	{
	  		/* determine if a collision occurs */
	         dum = (vr[ecolsp][j]*vr[ecolsp][j] + vth[ecolsp][j]*vth[ecolsp][j]
	  					+vph[ecolsp][j]*vph[ecolsp][j]);
	  			engy= Escale[ecolsp]*dum;
	  			vel = sqrt(dum);
	         ien = engy +0.5;
	   		if(ien >= NEMAX) ien = NEMAX - 1;
				stot	= sels[ien] + sext[ien] + sion[ien];
				prob   = 1 - exp(-dr*vel*gden*stot);
				random = frand();

				if( random < prob ) 
				{
					/* determine the type of collision and calculate new velocities, etc */
					random /= prob;
					sig1 = sels[ien]/(stot +DBL_MIN);
      	      sig2 = (sels[ien] + sext[ien])/(stot +DBL_MIN);
					
					/******************************************/
					/* normalize velocity to get direction */
					if (vel >1e-30)
					{
						vr[ecolsp][j]  /= vel;
						vth[ecolsp][j] /= vel;
						vph[ecolsp][j] /= vel;
					}	
					/******************************************/
					/* if the collision is elastic */
		           if (random <= sig1)
					{
 	               /* engy = engy*(1 - (1 - coschi)*2*m[ecolsp]/m[ecolsp+1]);*/
						newvel(engy,vel,&vr[ecolsp][j],&vth[ecolsp][j],&vph[ecolsp][j], 1);
					}
					/***********************************/
					/* if the collision is excitation */
	            else if ((random > sig1) && (random <= sig2))
		 			{
		            engy -= extengy0;
		            vel	= sqrt(fabs(engy)/Escale[ecolsp]);
		 				newvel(engy,vel,&vr[ecolsp][j],&vth[ecolsp][j],&vph[ecolsp][j], 0);
						 
				  	}
				  	/************************************/
				  	/* if the collision is ionization, 
				  		add the released electron first */
		         else
				  	{
		            engy -= ionengy0 - gtemp*log(frand()+DBL_MIN);
				  		rengy = 10.0*tan(frand()*atan(engy/20.0));
                  engy -= rengy;
		            
					   eengy = m[ionsp]*rengy/(m[ionsp] +m[ecolsp]);
                  vel = sqrt(fabs(eengy)/Escale[ecolsp]);
					 	k = np[ecolsp];
					 	vr[ecolsp][k]	= vr[ecolsp][j];
					  	vth[ecolsp][k] = vth[ecolsp][j];
					  	vph[ecolsp][k]	= vph[ecolsp][j];
                  newvel(eengy, vel, &vr[ecolsp][k], &vth[ecolsp][k], &vph[ecolsp][k], 0);
					 	r[ecolsp][k]	= r[ecolsp][j];
		              
						rengy -= eengy;
                  vel = sqrt(fabs(rengy)/Escale[ionsp]);
					   k = np[ionsp];
					   vr[ionsp][k] = frand();
					   vth[ionsp][k]= frand();
					   vph[ionsp][k]= frand();
         		   vel /= sqrt(vr[ionsp][k]*vr[ionsp][k] +vth[ionsp][k]*vth[ionsp][k]
						     +vph[ionsp][k]*vph[ionsp][k]);

						vr[ionsp][k]  *= vel;
						vth[ionsp][k] *= vel;
						vph[ionsp][k] *= vel;
	               r[ionsp][k] = r[ecolsp][j];

						if(++np[ionsp] >MAXLEN || ++np[ecolsp] >MAXLEN)
							ExitMsg("ADJUST(Ionization): too many particles. MUST EXIT!");

						/* then scatter the incident electron */
                  vel = sqrt(fabs(engy)/Escale[ecolsp]);
					 	newvel(engy,vel,&vr[ecolsp][j],&vth[ecolsp][j],&vph[ecolsp][j], 0);
				
						if ( r[ionsp][k] >= r0 && r[ionsp][k] <= r1)
						{
						   s= r[ionsp][k] -r0;

							if (r[ionsp][k] < r0+.5)
					         del = ((r[ionsp][k]+.5)*(r[ionsp][k]+.5)*(r[ionsp][k]+.5)-rg3[0])/
						            ((r[ionsp][k]+.5)*(r[ionsp][k]+.5)*(r[ionsp][k]+.5)-r0*r0*r0);
					   	else
						   {
				            del =((r[ionsp][k]+.5)*(r[ionsp][k]+.5)*(r[ionsp][k]+.5)-rg3[s])
							        /(3.*r[ionsp][k]*r[ionsp][k]+0.25);
						   }

							if (r[ionsp][k] > r1-.5)
					         del =(r1*r1*r1-r[ionsp][k]*r[ionsp][k]*r[ionsp][k])/
									  (r1*r1*r1-(r[ionsp][k]-.5)*(r[ionsp][k]-.5)*(r[ionsp][k]-.5));
						}
					   iontemp[s] += 1 - del;
					   iontemp[s+1] += del; 
					  
					}
			 	}	/* end if(electron-collision) */
			 } 		/* end for(electrons) 	   */
	}					/* end of if(collision)    */
	
	if(icollisional)
	{
		nnp = np[icolsp];
		for(j=0; j<nnp; j++)
		{
         /* determine if a collision occurs */
		   vel = sqrt(vr[icolsp][j]*vr[icolsp][j]+vth[icolsp][j]*vth[icolsp][j]
						 +vph[icolsp][j]*vph[icolsp][j]) +DBL_MIN;
			sig1= ascat*vel  +bscat;
			sig2= achrgx*vel +bchrgx;
			stot= sig1 +sig2;
			prob= 1 - exp(-stot);
			random= frand();