#include "peg_pic.h"


/* Calculation of parameters (for "isp" species) used in the simulation */
/* "fj0[][k]" is the particle flux from RHS (k=1) and LHS (k=0) walls   */

void specpar(int isp, float *fj0[2]) {
	float ttt, s;

	vxscale= max(vxscale, fabs(vpar0[isp])+vpart[isp]);
	vxscale= max(vxscale, fabs(vnpar0[isp][0])+vnpart[isp][0]);
	vxscale= max(vxscale, fabs(vnpar0[isp][1])+vnpart[isp][1]);
	
	if (sourcD) {
		vxscale= max(vxscale, fabs(vspar0[isp])+vspart[isp]);
	}


	/** Normalize all velocities **/
	vpart[isp] /= vscale;
	vpar0[isp] /= vscale;
	vnpart[isp][0] /= vscale;
	vnpar0[isp][0] /= vscale;
	vnpart[isp][1] /= vscale;
	vnpar0[isp][1] /= vscale;
	vperpt[isp]/= vscale;
	vperp0[isp]/= vscale;
	vnperpt[isp][0]/= vscale;
	vnperp0[isp][0]/= vscale;
	vnperpt[isp][1]/= vscale;
	vnperp0[isp][1]/= vscale;

	if ((!ht)&&(!ht1)) {
		ht_sp[isp]= 0;
	}

	if (slow || stopfl) {
		vx_min[isp] /= vscale;
		de_fvx[isp] /= vscale;
	}

	if (bmag>0.0) {
		/* If B!=0, particles injected at the wall have T_par=T_perp */
		if (((fj0[0][isp]>0.)&&(vnperpt[isp][0]!=vnpart[isp][0])) || ((fj0[1][isp]>0.)&&(vnperpt[isp][1]!=vnpart[isp][1]))) {
			puts("For magnetized particles injected at the wall T_par=T_perp\n");
			exit(1);
		}
	}

	/** Normalize other parameters **/
	Escale[isp] = .5*vscale*vscale*m[isp]/(1.6022e-19);

	if (slow || stopfl) {
		Qscale[isp] = .5*vscale*vscale*vscale*m[isp]*nc2p/wt[isp]/area/dx; /* energy flux */  
		Qwscale[isp] = Escale[isp]*nc2p*(1.6022e-19)/(nfft1*nstep[isp]*wt[isp]*dt*area);
		Tscale[isp] = 2*Escale[isp]/3.;
		emin[isp]  /= Escale[isp];
		de[isp]    /= Escale[isp];
	} else {
		if (Ef_flag || (qtp==3) || (qtp==4)) {
			Tscale[isp] = 2*Escale[isp]/3.;
		}
		
		if (Ef_flag) { 
			Qwscale[isp]= Escale[isp]*nc2p*(1.6022e-19)/(nfft1*dt*nstep[isp]*wt[isp]*area);
		}
		
		if (qtp==3) {
			Qscale[isp] = .5*vscale*vscale*vscale*m[isp]*nc2p/wt[isp]/area/dx;
		}
	}


	/** Conversions from physical to code units **/
	qm[isp] = q[isp]/m[isp];
	q[isp] *= nc2p/wt[isp]/area;
	fr_E[isp]= qm[isp]*dttx;               /* code will be faster */

	if (bmag>0.0) {
		fr_E[isp] *=0.5;
	}

	if (Ey_flag) {
		/* for Ey */
		Ey[isp] = 0.5*E_y*qm[isp]*dttx;

		if (x_E) {
			E1y[isp]= Ey[isp]/x_E;
		}
	}

	if (sourcD) {
		/** Parameters for the first source  **/
		vspart[isp] /= vscale;
		vspar0[isp] /= vscale;
		vsperpt[isp]/= vscale;
		vsperp0[isp]/= vscale;

		if (ht>0.) {
			if (bmag>0.0) {
				heat[isp] = exp(-vspart[isp]*nstep[isp]*ht*cs/L_s);
			} else {
				heat[isp] = exp(-vspart[isp]*nstep[isp]*ht/L_s);
			}
		}
	} else {
		ht= 0.;
		if (ccc) {
			/* Just for Landau log. */
			vspart[isp] /= vscale;
		}
	}

	if (sourcD1) {
		/** Parameters for the second source  **/
		vspart1[isp] /= vscale;
		vspar01[isp] /= vscale;
		vsperpt1[isp]/= vscale;
		vsperp01[isp]/= vscale;

		if (ht1>0.) {
			if (bmag>0.0) {
				heat1[isp] = exp(-vspart1[isp]*nstep[isp]*ht1*cs/L_s1);
			} else {
				heat1[isp] = exp(-vspart1[isp]*nstep[isp]*ht1/L_s1);
			}
		}
	} else {
		ht1= 0.;
	}

	if (isp < ccc) {
		ccpar(isp);                 /** Parameters for Coulomb collisions **/
	}

	/**  calculation of number of the injected particles **/
	enter[isp][0] = fj0[0][isp]*nstep[isp]*dt*area*wt[isp]/nc2p;
	enter[isp][1] = fj0[1][isp]*nstep[isp]*dt*area*wt[isp]/nc2p;

	if ((enter[isp][0]>0.0001) || (enter[isp][1]>0.0001)) {
		inject[isp] = 1;
	} else {
		inject[isp] = 0;
	}


	if (bmag>0.) {
		/** Parameters for the mover in presence of a B field **/
		ttt= tan(.5*bmag*qm[isp]*nstep[isp]*dt);
		s= 2*ttt/(1+ttt*ttt);
		tx[isp]= ttt*cs;
		tz[isp]= ttt*sn;
		sx[isp]= s*cs;
		sz[isp]= s*sn;
		tkv[isp]= (1.-ttt*ttt)/(1.+ttt*ttt);          /*  Let's make code faster */

		if (reflux > 1) {
			/* there is reflection or reinjection at the RHS wall */    
			i_t[isp]= 0;
		} else {
			if ((i_t[isp]==2) && (fabs(qm[isp])>0.)) {
				/* For "outsiders" */
				count_max[isp]=abs(2.*PI/(qm[isp]*bmag*dt*nstep[isp]));
			} else {
				count_max[isp]=0;
			}
		}
	} else {
		i_t[isp]= 0;       /* If B=0, a half-bounded plama model is simpler! */
	}

	srfpar(isp);                                  /** For surface physics **/

} /* End of SPECPAR */

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