#include "peg_pic.h"


/************** Maxwelian distribution *********************************/
/************** Used for initial loading *******************************/

void maxwelli(float *vx, float *vy, float *vz, float vtpar, float vtperp,
		float v0x, float v0z, int R, float Rcs, float Rsn) {
	float        theta, vtemp1, vtemp2;
	double       vmag;
	static int   init_flag=1;
	static float *parveli;

	if (init_flag) {
		/* do initial memory allocation and fill with values */
		parveli= (float *)malloc(VVMAX*sizeof(float));
		rdist(-5., 5., VVMAX, parveli, 1, 0.);       /* create distribution */
		init_flag =0;
		/* prtest(parveli, VVMAX);                        * testing */
	}

	vmag= vtperp*sqrt(2.*fabs(log(frand()+DBL_MIN)));
	theta= TWOPI*frand();

	if (bmag > 0.) {
		vtemp1 = vtpar*parveli[R];
		vtemp2 = vmag*cos(theta);

		*vx = vtemp1*Rcs - vtemp2*Rsn + v0x;
		*vy = vmag*sin(theta);
		*vz = vtemp1*Rsn + vtemp2*Rcs + v0z;
	} else {
		*vx = vtpar*parveli[R] + v0x;
		*vy = vmag*sin(theta);
		*vz = vmag*cos(theta) + v0z;
	}
}      /* END of maxwelli */


/******* Maxwelian distribution (better centered with respect to 0) ****/
/*********** Used for for particle source and for "heating" ************/ 
/*********** "ii" defines the sighn of the parallel (or if b=0, ********/
/*********** the x) component of the velocity. If ii= -1, 1, then ******/
/*********** it is negative, or positive, respectively. If ii=0, *******/
/*********** then it is defined randomly *******************************/

void maxwell(float *vx, float *vy, float *vz, float vtpar, float vtperp,
		float v0x, float v0z, int R, float Rcs, float Rsn, int ii) {
	float        theta, vtemp1, vtemp2;
	double       vmag;
	static int   init_flag=1;
	static float *parveli;


	if (init_flag) {
		/* do initial memory allocation and fill with values */
		parveli= (float *)malloc(VVMAX*sizeof(float));
		rdist(-5., 0., VVMAX, parveli, 1, 0.);       /* create distribution */
		init_flag =0;
		/*
    for(ii= VVMAX; ii < TVMAX; ii++)
     parveli[ii] = -parveli[TVMAX-ii-1];
		 */
	}

	vmag= vtperp*sqrt(2.*fabs(log(frand()+DBL_MIN)));
	theta= TWOPI*frand();

	if (!ii) {
		if (frand()>0.5) {
			ii= -1;
		} else {
			ii= 1;
		}
	}

	if (bmag > 0.) {
		vtemp1 = -ii*vtpar*parveli[R];
		vtemp2 = vmag*cos(theta);

		*vx = vtemp1*Rcs - vtemp2*Rsn + v0x;
		*vy = vmag*sin(theta);
		*vz = vtemp1*Rsn + vtemp2*Rcs + v0z;
	} else {
		*vx = -ii*vtpar*parveli[R] + v0x;
		*vy = vmag*sin(theta);
		*vz = vmag*cos(theta) + v0z;
	}
}      /* END of maxwell */


/************** Maxwelian distribution *********************************/
/*********** The fast one, used for neutral distribution ***************/

void fmaxwell(float *vx, float *vy, float *vz, float vtpar) {
	int          ii, R1, R2, R3;
	static int   init_flag=1;
	static float *parveli;


	if (init_flag) {
		/* do initial memory allocation and fill with values */
		parveli= (float *)malloc(1024*sizeof(float));
		rdist(-5., 0., 512, parveli, 1, 0.);       /* create distribution */

		for(ii= 512; ii < 1024; ii++) {
			parveli[ii] = -parveli[1023-ii];
		}

		init_flag =0;
	}

	R1 = 1024*frand();
	R2 = 1024*frand();
	R3 = 1024*frand();
	*vx = vtpar*parveli[R1];
	*vy = vtpar*parveli[R2];
	*vz = vtpar*parveli[R3];

}     /* END of fmaxwell */


/*********** (VxMaxwelian) distribution *********************************/
/******* Used for particle injection from the LHS wall ******************/

void vmaxwell0(float *vx, float *vy, float *vz, int iisp) {
	int          ii, k;
	float        theta, v0;
	double       vmag;
	static float *parveli;
	static int   init_flag=1;


	if (vnpart[iisp][0] > 0.0) {                      /* Hot beam */
		if (init_flag) {
			parveli= (float *)malloc(VVMAX*sizeof(float));
			
			for (k=0; k<nsp; k++) {
				if (vnpart[k][0] > 0.0) {
					v0= vnpar0[k][0]/vnpart[k][0];                  /* shift velocity */
					rdist(0., 5.+v0, VVMAX, parveli, 2, v0);       /* create distribution */

					for(ii= 0; ii < VVMAX; ii++) {
						parvell[k][ii] = vnpart[k][0]*parveli[ii];
					}
				}
			}
			
			init_flag =0;
			free(parveli);              /* was problem with saving of nonstatic float! */
		}

		k = VVMAX*frand();
		*vx = parvell[iisp][k];
		
		/* end of if(vnpart[iisp][0] > 0.0) loop */
	} else { 
        /* Cold beam */
		*vx = vnpar0[iisp][0];
	}

	vmag= vnperpt[iisp][0]*sqrt(2.*fabs(log(frand()+DBL_MIN)));
	theta= TWOPI*frand();

	*vy = vmag*sin(theta);
	*vz = vmag*cos(theta) + vnperp0[iisp][0];

}      /* END of vmaxwell0 */


/*********** (VxMaxwelian) distribution *********************************/
/******* Used for particle injection from the RHS wall ******************/

void vmaxwell1(float *vx, float *vy, float *vz, int iisp){
	int          ii, k;
	float        theta, v0;
	double       vmag;
	static int   init_flag=1;
	static float *parveli;


	if (vnpart[iisp][1] > 0.0) {
		/* Hot beam */
		if (init_flag) {
			parveli= (float *)malloc(VVMAX*sizeof(float));
			
			for (k=0; k<nsp; k++) {
				if (vnpart[k][1] > 0.0) {
					v0= vnpar0[k][1]/vnpart[k][1];                 /* shift velocity */
					rdist(0., 5.+v0, VVMAX, parveli, 2, v0);       /* create distribution */

					for(ii= 0; ii < VVMAX; ii++) {
						parvelr[k][ii] = -vnpart[k][1]*parveli[ii];
					/* prtest(parveli, VVMAX);                        * testing */
					}
				}
			}

			init_flag =0;
			free(parveli);              /* was problem with saving of nonstatic float! */
		}

		k = VVMAX*frand();
		*vx = parvelr[iisp][k];

		/* end of if(vnpart[iisp][0] > 0.0) loop */
	} else {  
        /* Cold beam */
		*vx = -vnpar0[iisp][1];
	}

	vmag= vnperpt[iisp][1]*sqrt(2.*fabs(log(frand()+DBL_MIN)));
	theta= TWOPI*frand();

	*vy = vmag*sin(theta);
	*vz = vmag*cos(theta) + vnperp0[iisp][1];

}      /* END of vmaxwell1 */

/**************************************************************************/
/************** Velocity distribution for the secondary particles *********/

void secdist(float *vx, float *vy, float *vz, int injtype) {
	int          ii, k;
	float        theta;
	double       vmag;
	static int   init_flag=1;
	static float *parveli;


	if (init_flag) {
		/* do initial memory allocation and fill with values */
		for (k=0; k<inj_max; k++) {
			parveli= (float *)malloc(VVMAX*sizeof(float));
			rdist(0., 5., VVMAX, parveli, injtype+2, 0.); /* create distribution */

			for(ii= 0; ii < VVMAX; ii++) {
				parvels[k][ii] = parveli[ii];
			}
		}

		init_flag =0;
		free(parveli);             /* was problem with saving of nonstatic float! */
	}

	k = VVMAX*frand();
	*vx = parvels[injtype][k];

	vmag= sqrt(2.*fabs(log(frand()+DBL_MIN)));
	theta= TWOPI*frand();

	*vy = vmag*sin(theta);
	*vz = vmag*cos(theta);
}
/* END of secdist() */                                                                                                                                                             
/************************************************************************/
