#include "peg_pic.h"


/********* Calculation of n, V, T, Q profiles ******************/
/*** prtq flag *******************  calculated are *************/
/***   0       ***  no profiles                              ***/
/***   1       ***  all profiles                             ***/
/***   2       ***  all profiles except T_par-T_perp         ***/
/***   3       ***  temperature and heat flux profiles       ***/
/***   4       ***  temperature profile                      ***/
/***************************************************************/
/***   If prtq=0 and Ef_flag=1, routine still calculates     ***/
/*** temperatures at two specified points.                   ***/


void avq(int iisp, int prtd) {
	int   i, j=0, k;
	float s, Vx_temp, Vy_temp, Vz_temp, Tx_temp, Qh_temp, Qtot_temp;
	float srho_t;

	
	if (prtd==1) {
		/* All profiles */
#pragma omp parallel for private(i)
		for (i=0; i<=nc; i++) {
			srho[iisp][i] = 0.;
			Vx[iisp][i] = 0.;
			Vy[iisp][i] = 0.;
			Vz[iisp][i] = 0.;
			Tx[iisp][i] = 0.;
			Qh[iisp][i] = 0.;
			Qtot[iisp][i] = 0.;
			T_dif[iisp][i] = 0.;
		}
	} else if (prtd==2) {
		/* All profiles except T_par-T_perp */
#pragma omp parallel for private(i)
		for (i=0; i<=nc; i++) {
			srho[iisp][i] = 0.;
			Vx[iisp][i] = 0.;
			Vy[iisp][i] = 0.;
			Vz[iisp][i] = 0.;
			Tx[iisp][i] = 0.;
			Qh[iisp][i] = 0.;
			Qtot[iisp][i] = 0.;
		}
	} else if (prtd==3) {
		/* Temperature and heat flux profiles */
#pragma omp parallel for private(i)
		for (i=0; i<=nc; i++) {
			srho[iisp][i] = 0.;
			Vxef[i] = 0.;
			Vyef[i] = 0.;
			Vzef[i] = 0.;
			Tx[iisp][i] = 0.;
			Qtot[iisp][i] = 0.;
			Qh[iisp][i] = 0.;
		}
	} else if (prtd==4) {
		/* Temperature profile */
#pragma omp parallel for private(i)
		for (i=0; i<=nc; i++) {
			srho[iisp][i] = 0.;
			Vxef[i] = 0.;
			Vyef[i] = 0.;
			Vzef[i] = 0.;
			Tx[iisp][i] = 0.;
		}
	}

	if (prtd==0) {
		/* Accumulate n */
		for (i=np[iisp]-1; i>=0; i--) {
			if (x[iisp][i] > 0. && x[iisp][i] < xnc) {
				j = x[iisp][i];
				s = x[iisp][i] - j;
				srho[iisp][j] += 1. - s;
				srho[iisp][j+1] += s;
			}
		}

		if (reflux>1) {
			/* case with reflection, or reinjection */
			srho[iisp][0] += srho[iisp][nc];
			srho[iisp][nc] = srho[iisp][0];
		} else {
			/* adjust the 2 HALF-cells at ends */
			srho[iisp][0]  *= 2.;
			srho[iisp][nc] *= 2.;
		}
	}     /* end of if (prtd==0) */


	if ((prtd==1)||(prtd==2)) {
		for (i=np[iisp]-1; i>=0; i--) {
			/* Accumulate n and V */
			if (x[iisp][i] > 0. && x[iisp][i] < xnc) {
				j = x[iisp][i];
				s = x[iisp][i] - j;
				srho[iisp][j] += 1. - s;
				srho[iisp][j+1] += s;

				Vx[iisp][j] += (1. - s)*vx[iisp][i];
				Vx[iisp][j+1] += s*vx[iisp][i];
				Vy[iisp][j] += (1. - s)*vy[iisp][i];
				Vy[iisp][j+1] += s*vy[iisp][i];
				Vz[iisp][j] += (1. - s)*vz[iisp][i];
				Vz[iisp][j+1] += s*vz[iisp][i];
			}
		}

#pragma omp parallel for shared(srho,Vx,Vy,Vz,iisp,nc) private(j)
		for (j=0; j<=nc; j++) {
			if (srho[iisp][j] > 1.0e-10) {
				Vx[iisp][j] /= srho[iisp][j];
				Vy[iisp][j] /= srho[iisp][j];
				Vz[iisp][j] /= srho[iisp][j];
			} else {
				Vx[iisp][j] = 0.;
				Vy[iisp][j] = 0.;
				Vz[iisp][j] = 0.;
			}
		}

		for (i=np[iisp]-1; i>=0; i--)  {
			/* Accumulate T, Q */
			if (x[iisp][i] > 0. && x[iisp][i] < xnc) {
				j = x[iisp][i];
				s = x[iisp][i] - j;

				Vx_temp= vx[iisp][i] - (1.-s)*Vx[iisp][j] - s*Vx[iisp][j+1];
				Vy_temp= vy[iisp][i] - (1.-s)*Vy[iisp][j] - s*Vy[iisp][j+1];
				Vz_temp= vz[iisp][i] - (1.-s)*Vz[iisp][j] - s*Vz[iisp][j+1];
				Tx_temp= Vx_temp*Vx_temp + Vy_temp*Vy_temp + Vy_temp*Vy_temp;
				Qtot_temp= vx[iisp][i]*(vx[iisp][i]*vx[iisp][i]+vy[iisp][i]*vy[iisp][i]+vz[iisp][i]*vz[iisp][i]);
				Qh_temp= Vx_temp*Tx_temp;
				Tx[iisp][j] += (1. - s)*Tx_temp;
				Tx[iisp][j+1] += s*Tx_temp;
				Qh[iisp][j] += (1. - s)*Qh_temp;
				Qh[iisp][j+1] += s*Qh_temp;
				Qtot[iisp][j] += (1. - s)*Qtot_temp;
				Qtot[iisp][j+1] += s*Qtot_temp;

				if (prtd==1) {
					Tx_temp= Vx_temp*Vx_temp - (Vy_temp*Vy_temp + Vz_temp*Vz_temp)/2.;
					T_dif[iisp][j] += (1. - s)*Tx_temp;
					T_dif[iisp][j+1] += s*Tx_temp;
				}
			}
		}

		for (j=nc; j>=0; j--) {
			if (bmag>0) {
				Vpar[iisp][j] = Vx[iisp][j]*cs + Vz[iisp][j]*sn;
			}
			
			Qtot[iisp][j] *= Qscale[iisp];
			Qh[iisp][j] *= Qscale[iisp];

			if (srho[iisp][j] > 1.0e-10) {
				Tx[iisp][j] *= Tscale[iisp]/srho[iisp][j];
				
				if (prtd==1)  {
					T_dif[iisp][j] *= Escale[iisp]*2./srho[iisp][j];
				}
			} else {
				Tx[iisp][j] = 0.;
				
				if (prtd==1) {
					T_dif[iisp][j] = 0.;
				}
			}
		}

		if (reflux>1) {
			/* case with reflection, or reinjection */
			srho[iisp][0] += srho[iisp][nc];
			srho[iisp][nc] = srho[iisp][0];
		} else {
			/* adjust the 2 HALF-cells at ends */
			srho[iisp][0]  *= 2.;
			srho[iisp][nc] *= 2.;
		}

		Tlhs_temp[iisp] += Tx[iisp][potpos[0]];
		Trhs_temp[iisp] += Tx[iisp][potpos[1]];

	/* end of if (prtd<3) */
	} else if (prtd==3 || prtd==4) {
		for (i=np[iisp]-1; i>=0; i--) {
			/* Accumulate n and V */
			if (x[iisp][i] > 0. && x[iisp][i] < xnc) {
				j = x[iisp][i];
				s = x[iisp][i] - j;
				srho[iisp][j] += 1. - s;
				srho[iisp][j+1] += s;

				Vxef[j] += (1. - s)*vx[iisp][i];
				Vxef[j+1] += s*vx[iisp][i];
				Vyef[j] += (1. - s)*vy[iisp][i];
				Vyef[j+1] += s*vy[iisp][i];
				Vzef[j] += (1. - s)*vz[iisp][i];
				Vzef[j+1] += s*vz[iisp][i];
			}
		}

#pragma omp parallel for shared(srho,Vxef,Vyef,Vzef,iisp,nc) private(j)
		for (j=0; j<=nc; j++) {
			if (srho[iisp][j] > 1.0e-10) {
				Vxef[j] /= srho[iisp][j];
				Vyef[j] /= srho[iisp][j];
				Vzef[j] /= srho[iisp][j];
			} else {
				Vxef[j] = 0.;
				Vyef[j] = 0.;
				Vzef[j] = 0.;
			}
		}

		for (i=np[iisp]-1; i>=0; i--) {
			/* Accumulate T and Q */
			if (x[iisp][i] > 0. && x[iisp][i] < xnc) {
				j = x[iisp][i];
				s = x[iisp][i] - j;

				Vx_temp= vx[iisp][i] - (1.-s)*Vxef[j] - s*Vxef[j+1];
				Vy_temp= vy[iisp][i] - (1.-s)*Vyef[j] - s*Vyef[j+1];
				Vz_temp= vz[iisp][i] - (1.-s)*Vzef[j] - s*Vzef[j+1];
				Tx_temp= Vx_temp*Vx_temp + Vy_temp*Vy_temp + Vy_temp*Vy_temp;
				Tx[iisp][j] += (1. - s)*Tx_temp;
				Tx[iisp][j+1] += s*Tx_temp;

				if (prtd==3) {
					Qtot_temp= vx[iisp][i]*(vx[iisp][i]*vx[iisp][i]+vy[iisp][i]*vy[iisp][i]+vz[iisp][i]*vz[iisp][i]);
					Qh_temp= Vx_temp*Tx_temp;
					Qh[iisp][j] += (1. - s)*Qh_temp;
					Qh[iisp][j+1] += s*Qh_temp;
					Qtot[iisp][j] += (1. - s)*Qtot_temp;
					Qtot[iisp][j+1] += s*Qtot_temp;
				}
			}
		}

		if (prtd==3) {
			for (j=nc; j>=0; j--) {
				if (bmag>0) {
					Vpar[iisp][j] = (Vxef[j]*cs + Vzef[j]*sn)*vscale;
				} else {
					Vpar[iisp][j] = Vxef[j]*vscale;
				}
				
				Qtot[iisp][j] *= Qscale[iisp];
				Qh[iisp][j] *= Qscale[iisp];

				if (srho[iisp][j] > 1.0e-10) {
					Tx[iisp][j] *= Tscale[iisp]/srho[iisp][j];
				} else {
					Tx[iisp][j] = 0.;
				}
			}
		} else {
			if (srho[iisp][j] > 1.0e-10) {
				Tx[iisp][j] *= Tscale[iisp]/srho[iisp][j];
			} else {
				Tx[iisp][j] = 0.;
			}
		}

		if (reflux>1) {
			/* case with reflection, or reinjection */
			srho[iisp][0] += srho[iisp][nc];
			srho[iisp][nc] = srho[iisp][0];
		} else {
			/* adjust the 2 HALF-cells at ends */
			srho[iisp][0]  *= 2.;
			srho[iisp][nc] *= 2.;
		}

		if (Ef_flag) {
			Tlhs_temp[iisp] += Tx[iisp][potpos[0]];
			Trhs_temp[iisp] += Tx[iisp][potpos[1]];
		}
		/* end of else if (prtd==3 || prtd==4) */
	} else if (Ef_flag) {
		/* Themperature at two specified points */
		for (k=0; k<2; k++) {
			Vx_temp= 0.;
			Vy_temp= 0.;
			Vz_temp= 0.;
			srho_t= 0.;
			j= potpos[k];

			for (i=np[iisp]-1; i>=0; i--) {
				/* Accumulate V */
				if (j==(int)x[iisp][i]) {
					Vx_temp += vx[iisp][i];
					Vy_temp += vy[iisp][i];
					Vz_temp += vz[iisp][i];
					srho_t += 1.;
				}
			}

			if (srho_t > 0.5) {
				Vx_temp /= srho_t;
				Vy_temp /= srho_t;
				Vz_temp /= srho_t;
				Tx_temp =0.;

				for (i=np[iisp]-1; i>=0; i--) {
					/* Accumulate T */
					if (j==(int)x[iisp][i]) {
						Qh_temp = vx[iisp][i] - Vx_temp;
						Qtot_temp = vy[iisp][i] - Vy_temp;
						s = vz[iisp][i] - Vz_temp;
						Tx_temp += Qh_temp*Qh_temp + Qtot_temp*Qtot_temp + s*s;
					}
				}

				if (k==0) {
					Tlhs_temp[iisp] += (Tscale[iisp]*Tx_temp/srho_t);
				} else {
					Trhs_temp[iisp] += (Tscale[iisp]*Tx_temp/srho_t);
				}

			} /* end of if (srho_t > 0.5) */
		} /* end of for (k=0; k<2; k++) loop */
	}  /* end of if (Ef_flag) */

#pragma omp parallel for shared(srho, wt, nc, iisp) private(i)
	for (i=0; i<=nc; i++) {
		srho[iisp][i] /= wt[iisp];
	}

}  /* END of avq() */

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