#include "peg_pic.h"


/**** Read and check input parameters for given species (isp).     *****/
/**** "initnp" is the initial number of particles (of given species) ***/
/* "fj0[][k]" is the particle flux from RHS (k=1) and LHS (k=0) walls  */

void species(int *initnp, float *fj0[2], int isp) {
	int   i, nb1, nb2, k, loop;
	char aachar[512];
	float iinitn, j0l, j0r, vtpar, v0x, v0z, vtperp, temin, temax;
	float vtlpar, vtlperp, vlx, vlz, vx_max, con_sr, con_sr1;
	float vtrpar, vtrperp, vrx, vrz,  vtspar, vtsperp, vsx, vsz, xd1;
	float vtspar1, vtsperp1, vsx1, vsz1;

#ifdef REFLECTED_OWN_SPECIES
	int real_nsp = (int) nsp / 2;
	int vsp = isp + real_nsp;			/* index of "virtual species" index */

	if( vsp > nsp - 1 ) {
		/* "virtual species" (one defined because of reflection), nothing to be done; return back */
		return;
	}
#endif

	/*** Read in SPECIES parameters... ****/

	while (fscanf(InputDeck, "%e %e %d %d %d", &q[isp], &m[isp], &nstep[isp], &wt[isp], &ht_sp[isp]) < 5) {
		fscanf(InputDeck, "%s", aachar);
	}

#ifdef REFLECTED_OWN_SPECIES
	/* copy values to "virtual" species */
	q[vsp] 		= q[isp];
	m[vsp] 		= m[isp];
	nstep[vsp] 	= nstep[isp];
	wt[vsp] 	= wt[isp];
	ht_sp[vsp] 	= ht_sp[isp];
#endif

	while (fscanf(InputDeck, "%e %e %e %e %e", &iinitn, &vtpar, &vtperp, &v0x, &v0z) < 5) {
		fscanf(InputDeck, "%s", aachar);
	}

	while (fscanf(InputDeck, "%e %e %e %e %e", &j0l, &vtlpar, &vtlperp, &vlx, &vlz) < 5) {
		fscanf(InputDeck, "%s", aachar);
	}

	while (fscanf(InputDeck, "%e %e %e %e %e %d", &j0r, &vtrpar, &vtrperp, 	&vrx, &vrz, &i_t[isp]) < 6) {
		fscanf(InputDeck, "%s", aachar);
	}

#ifdef REFLECTED_OWN_SPECIES
	/* copy values to "virtual" species */
	i_t[vsp] = i_t[isp];
#endif

	while (fscanf(InputDeck, "%e %e %e %e %e", &con_sr, &vtspar, &vtsperp, &vsx, &vsz) < 5) {
		fscanf(InputDeck, "%s", aachar);
	}

	loop = 1;
	/* default is to use same factor for Source1 as is defined in Source */
	con_sr1 = con_sr;

	do {
		/* next line could have 4 or 5 values to read (4=old file, 5=new) */
		/* first read whole line as string and then try to read values from line */
		fgets(aachar, 511, InputDeck);

		if( sscanf(aachar, "%e %e %e %e %e", &con_sr1, &vtspar1, &vtsperp1, &vsx1, &vsz1) == 5) {
			/* all OK, new input file */
			loop = 0;

		} else if( sscanf(aachar, "%e %e %e %e", &vtspar1, &vtsperp1, &vsx1, &vsz1) == 4) {
			/* all OK, old input file = use same factor for Source1 as is defined in Source */
			loop = 0;
		}
	} while( loop );

	while (fscanf(InputDeck, "%d %e %e %d %e %e %e %e %e %e", &nb1, &temin,
			&temax, &nb2, &xd1, &vx_max, &gridref_l0[isp], &gridref_lsramp[isp], &gridref_r0[isp], 
			&gridref_rsramp[isp] ) <10) {
		fscanf(InputDeck, "%s", aachar);
	}

#ifdef REFLECTED_OWN_SPECIES
	/* copy values to "virtual" species */
	gridref_l0[vsp]		= gridref_l0[isp];
	gridref_lsramp[vsp] = gridref_lsramp[isp];
	gridref_r0[vsp] 	= gridref_r0[isp];
	gridref_rsramp[vsp] = gridref_rsramp[isp];
#endif

	while (fscanf(InputDeck, "%d %d %d %e %e %e %e", &Nsec[isp], &secsp[isp][0],
			&sectype[isp][0], &seec[isp][0][0], &seec[isp][0][1],
			&Ew[isp][0][0], &Ew[isp][0][1]) <7) {
		fscanf(InputDeck, "%s", aachar);
	}

#ifdef REFLECTED_OWN_SPECIES
	/* copy values to "virtual" species */
	Nsec[vsp]		= Nsec[isp];
	secsp[vsp][0] 	= secsp[isp][0];
	sectype[vsp][0]	= sectype[isp][0];
	seec[vsp][0][0] = seec[isp][0][0];
	seec[vsp][0][1] = seec[isp][0][1];
	Ew[vsp][0][0] 	= Ew[isp][0][0];
	Ew[vsp][0][1] 	= Ew[isp][0][1];
#endif

	while (fscanf(InputDeck, "%e %e %d %d %e %e", &E0[isp][0][0], &E0[isp][0][1],      
			&injtype[isp][0][0], &injtype[isp][0][1],
			&Tseec[isp][0][0], &Tseec[isp][0][1]) <6) {
		fscanf(InputDeck, "%s", aachar);
	}

#ifdef REFLECTED_OWN_SPECIES
	/* copy values to "virtual" species */
	E0[vsp][0][0] 		= E0[isp][0][0];
	E0[vsp][0][1] 		= E0[isp][0][1];
	injtype[vsp][0][0] 	= injtype[isp][0][0];
	injtype[vsp][0][1] 	= injtype[isp][0][1];
	Tseec[vsp][0][0] 	= Tseec[isp][0][0];
	Tseec[vsp][0][1] 	= Tseec[isp][0][1];
#endif
	for (i=1; i < Nsec[isp]; i++) {
		while (fscanf(InputDeck, "%d %d %e %e %e %e", &secsp[isp][i],
				&sectype[isp][i], &seec[isp][i][0], &seec[isp][i][1],
				&Ew[isp][i][0], &Ew[isp][i][1]) <6) {
			fscanf(InputDeck, "%s", aachar);
		}

#ifdef REFLECTED_OWN_SPECIES
		/* copy values to "virtual" species */
		secsp[vsp][i] 	= secsp[isp][i];
		sectype[vsp][i]	= sectype[isp][i];
		seec[vsp][i][0] = seec[isp][i][0];
		seec[vsp][i][1] = seec[isp][i][1];
		Ew[vsp][i][0] 	= Ew[isp][i][0];
		Ew[vsp][i][1] 	= Ew[isp][i][1];
#endif

		while (fscanf(InputDeck, "%e %e %d %d %e %e", &E0[isp][i][0], &E0[isp][i][1],      
				&injtype[isp][i][0], &injtype[isp][i][1],
				&Tseec[isp][i][0], &Tseec[isp][i][1]) <6) {
			fscanf(InputDeck, "%s", aachar);
		}

#ifdef REFLECTED_OWN_SPECIES
		/* copy values to "virtual" species */
		E0[vsp][i][0] 		= E0[isp][i][0];
		E0[vsp][i][1]		= E0[isp][i][1];
		injtype[vsp][i][0]	= injtype[isp][i][0];
		injtype[vsp][i][1]	= injtype[isp][i][1];
		Tseec[vsp][i][0] 	= Tseec[isp][i][0];
		Tseec[vsp][i][1] 	= Tseec[isp][i][1];
#endif
	}

	/* Assign input parameters to arguments */
	fj0[0][isp] = fabs(j0l);
	fj0[1][isp] = fabs(j0r);

#ifdef REFLECTED_OWN_SPECIES
	fj0[0][vsp] = fj0[0][isp];
	fj0[1][vsp] = fj0[1][isp];
#endif

	vpart[isp]  = sqrt(fabs(vtpar*1.6022e-19/m[isp]));
	vperpt[isp] = sqrt(fabs(vtperp*1.6022e-19/m[isp]));
	vnpart[isp][0]  = sqrt(fabs(vtlpar*1.6022e-19/m[isp]));
	vnperpt[isp][0] = sqrt(fabs(vtlperp*1.6022e-19/m[isp]));
	vnpart[isp][1]  = sqrt(fabs(vtrpar*1.6022e-19/m[isp]));
	vnperpt[isp][1] = sqrt(fabs(vtrperp*1.6022e-19/m[isp]));

#ifdef REFLECTED_OWN_SPECIES
	vpart[vsp]  = sqrt(fabs(vtpar*1.6022e-19/m[vsp]));
	vperpt[vsp] = sqrt(fabs(vtperp*1.6022e-19/m[vsp]));
	vnpart[vsp][0]  = sqrt(fabs(vtlpar*1.6022e-19/m[vsp]));
	vnperpt[vsp][0] = sqrt(fabs(vtlperp*1.6022e-19/m[vsp]));
	vnpart[vsp][1]  = sqrt(fabs(vtrpar*1.6022e-19/m[vsp]));
	vnperpt[vsp][1] = sqrt(fabs(vtrperp*1.6022e-19/m[vsp]));
#endif

	vpar0[isp] 		= v0x;
	vperp0[isp] 	= v0z;
	vnpar0[isp][0]  = vlx;
	vnperp0[isp][0] = vlz;
	vnpar0[isp][1]  = -vrx;
	vnperp0[isp][1] = vrz;
	initnp[isp] 	= iinitn*length*area*wt[isp]/nc2p +.5;
	Fscale[isp] 	= nc2p/(nfft1*nstep[isp]*dt*wt[isp]*area);

#ifdef REFLECTED_OWN_SPECIES
	vpar0[vsp] 		= vpar0[isp];
	vperp0[vsp] 	= vperp0[isp];
	vnpar0[vsp][0]  = vnpar0[isp][0];
	vnperp0[vsp][0] = vnperp0[isp][0];
	vnpar0[vsp][1]  = vnpar0[isp][1];
	vnperp0[vsp][1] = vnperp0[isp][1];
	initnp[vsp] 	= iinitn*length*area*wt[vsp]/nc2p +.5;
	Fscale[vsp] 	= nc2p/(nfft1*nstep[vsp]*dt*wt[vsp]*area);
#endif


	for (i=0; i < Nsec[isp]; i++) {
		if (Nsec[isp]>nsp) {
			puts("Number of secondaries exceeds number of species!\n");
			exit(1);
		}

		if (!secsp[isp][i]) {
			puts("SecSp = 0 !\n");
			exit(1);
		}

		if (secsp[isp][i] > nsp) {
			puts("SecSp exceeds number of species!\n");
			exit(1);
		}

		secsp[isp][i] --;

#ifdef REFLECTED_OWN_SPECIES
		secsp[vsp][i] --;
#endif

		for (k=0; k<2; k++) {
			/* T --> V */
			Tseec[isp][i][k] = sqrt(fabs(Tseec[isp][i][k]*1.6022e-19/m[isp]));
#ifdef REFLECTED_OWN_SPECIES
			Tseec[vsp][i][k] = sqrt(fabs(Tseec[vsp][i][k]*1.6022e-19/m[vsp]));
#endif
		}
	}

	if (sourcD) {
		/* The first source */
		Src_con[isp] = con_sr*Src*wt[isp];
		vspart[isp]  = sqrt(fabs(vtspar*1.6022e-19/m[isp]));
		vsperpt[isp] = sqrt(fabs(vtsperp*1.6022e-19/m[isp]));
		vspar0[isp]  = vsx;
		vsperp0[isp] = vsz;

#ifdef REFLECTED_OWN_SPECIES
		Src_con[vsp] = con_sr*Src*wt[vsp];
		vspart[vsp]  = sqrt(fabs(vtspar*1.6022e-19/m[vsp]));
		vsperpt[vsp] = sqrt(fabs(vtsperp*1.6022e-19/m[vsp]));
		vspar0[vsp]  = vsx;
		vsperp0[vsp] = vsz;
#endif

	} else if (ccc) {
		/* For Landau log. */
		vspart[isp]  = sqrt(fabs(vtspar*1.6022e-19/m[isp]));
#ifdef REFLECTED_OWN_SPECIES
		vspart[vsp]  = sqrt(fabs(vtspar*1.6022e-19/m[vsp]));
#endif
	}

	if (sourcD1) {
		/* The second source */
		Src_con1[isp] = con_sr1*Src1*wt[isp];
		vspart1[isp]  = sqrt(fabs(vtspar1*1.6022e-19/m[isp]));
		vsperpt1[isp] = sqrt(fabs(vtsperp1*1.6022e-19/m[isp]));
		vspar01[isp]  = vsx1;
		vsperp01[isp] = vsz1;

#ifdef REFLECTED_OWN_SPECIES
		Src_con1[vsp] = con_sr1*Src1*wt[vsp];
		vspart1[vsp]  = sqrt(fabs(vtspar1*1.6022e-19/m[vsp]));
		vsperpt1[vsp] = sqrt(fabs(vtsperp1*1.6022e-19/m[vsp]));
		vspar01[vsp]  = vsx1;
		vsperp01[vsp] = vsz1;
#endif
	}


	if (slow || stopfl) {
		if (nb1<=1) {
			nb1 = 2;
		}
		
		if (nb2<=1) {
			nb2 = 2;
		}
		
		if (temin >= temax) {
			puts("Emin >= Emax!\n");
			exit(1);
		}
		
		if (xd1 >= vx_max) {
			puts("vxmin >= vxmax!\n");
			exit(1);
		}
		
		nbin[isp]   = nb1;
		nbin_fvx[isp] = nb2;
		vx_min[isp] = xd1;
		emin[isp]   = temin;
		de[isp]     = (temax-temin)/(nb1-1);
		dtheta[isp] = PI/2/(nb1-1);
		de_fvx[isp] = (vx_max-xd1)/nb2;

#ifdef REFLECTED_OWN_SPECIES
		nbin[vsp]   = nbin[isp];
		nbin_fvx[vsp] = nbin_fvx[isp];
		vx_min[vsp] = vx_min[isp];
		emin[vsp]   = emin[isp];
		de[vsp]     = de[isp];
		dtheta[vsp] = dtheta[isp];
		de_fvx[vsp] = de_fvx[isp];
#endif
	}
}   /* End of SPECIES */

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