#include "peg_pic.h"

/********* distribution functions used in the simulation  *************/

double fmaxw(double, float);
double vfmaxw(double, float);


/**************************************************************************/
/** Calculates distribution of "vx" from "vmin" to "vmax" using the *******/
/** distribution function, "vn" is the size of "vx" array (i.e. number of */
/** grids in vx space), "k" describes the distribution function, if k=1 ***/
/** then, distribution is maxwellian, if k=2 then distribution is *********/
/** Vxmaxwellian, and so on. "v0" is some parameter, for example the ******/
/** shift velocity. *******************************************************/

void rdist(float vmin, float vmax, int vn, float *vx, int k, float v0) {
	int    ii;
	float  v1;
	double vv, ddv, sf, ff, ddf, sumf, ss, ds;
	double (*dfptr)(double, float);

	if (k==1) {
		dfptr = fmaxw;                       /* Maxwellian */
	} else {
		dfptr = vfmaxw;                      /* VxMaxwellian */
	}
	/* ...  other functions */

	v1 = vmin;
	ddv = 0.5*(vmax-v1)/vn;
	vv = v1;
	sf = 0.5*(*dfptr)(vv, v0);

	for (ii=0; ii<vn; ii++) {
		vv += ddv;
		ddf = (*dfptr)(vv, v0);
		sf += ddf + ddf;
		vv += ddv;
		ddf = (*dfptr)(vv, v0);
		sf += ddf;
	}
	
	ddf = 0.5*(*dfptr)(vmax, v0);
	sf -= ddf;
	sf = (sf + sf)*ddv/3.;             /* Cumulative function */

	sumf = 0.;
	ii = 0;
	vv = v1;
	ds = sf/vn;
	ss = ds/2.;

	L12:    ff = (*dfptr)(vv, v0);
	ddv = 1./128.;
	
	if (ff > 0.) {
		ddv = min(ddv, .45*ds/ff);
	}
	
	do {
		ddf = (*dfptr)(vv, v0) + 4.*(*dfptr)(vv+ddv, v0) + (*dfptr)(vv+ddv+ddv, v0);
		ddf *= ddv/3.;
		if (ddf >= ds) ddv *= 0.5;
	} while (ddf >= ds);


	vv += ddv+ddv;
	sumf += ddf;
	
	if (sumf < ss) {
		goto L12;
	}


	/* When a desired velocity has been passed, interpolate
           to find it, and start looking for the next velocity */

	vx[ii] = vv - (sumf - ss)*ddv/ddf;
	ii++;
	ss += ds;
	
	if (ii < vn) {
		goto L12;
	}

}  /* END of rdist() */


/***********************************************************************/
/** Maxwell distribution used in maxwelli(), maxwell() and fmaxwell() **/

double fmaxw(double vv, float vemty) {
	return(exp(-vv*vv*0.5));
}     
/* END of fmaxw() */    


/**********************************************************************/
/********** VxMaxwell function used in vmaxwell()  ********************/

double vfmaxw(double vv, float v0) {
	double u;

	u = vv - v0;
	return(vv*exp(-u*u*0.5));
}    
/* END of vfmaxw() */   


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