#include "peg_pic.h"


/******* Calculation of time history of different values **********/

void chist(float *jt) {
	int i, isp;

	for (isp=0; isp<nsp; isp++) {
		/* caclulate all values that are species dependent */
		np_hist[isp][hist_hi] = np[isp];
		gridref_l_hist[isp][hist_hi] = gridref_l[isp];
		gridref_r_hist[isp][hist_hi] = gridref_r[isp];

		if (slow) {
			jwall_hist[isp][hist_hi] = jt[isp];
			jt[isp]= 0.;

			if(hist_hi) {
				jwall_hist[isp][hist_hi] +=jwall_hist[isp][hist_hi-1];
			}

			F_hist[isp][hist_hi] = FLHS[isp];
			Eflhs_hist[isp][hist_hi] = Eflhs[isp];
			F_histRHS[isp][hist_hi] = FRHS[isp];
			Efrhs_hist[isp][hist_hi] = Efrhs[isp];
			Tlhs_hist[isp][hist_hi] = Tlhs[isp];
			Trhs_hist[isp][hist_hi] = Trhs[isp];

			if (T_diag) {
				Tx_hist[isp][hist_hi] = 0.;
				T_dif_hist[isp][hist_hi] = 0.;
				Vxt[isp][hist_hi] = 0.;

				for (i=0; i<=nc; i++) {
					Tx_hist[isp][hist_hi] += Tx[isp][i];
					T_dif_hist[isp][hist_hi] += T_dif[isp][i];
					Vxt[isp][hist_hi] +=Vx[isp][i];
				}

				Tx_hist[isp][hist_hi] /= (xnc +1.);
				T_dif_hist[isp][hist_hi] /= (xnc +1.);
				Vxt[isp][hist_hi] /= (xnc +1.);

			}     /* end of if (T_diag) */

			kes_hist[isp][hist_hi] = 0;

			for (i=0; i<np[isp]; i++) {
				kes_hist[isp][hist_hi] += vx[isp][i]*vx[isp][i] +vy[isp][i]*vy[isp][i]+vz[isp][i]*vz[isp][i];
			}

			kes_hist[isp][hist_hi] = kes_hist[isp][hist_hi]*Escale[isp]/(np[isp]+1e-12) + DBL_MIN;
			/*  end of if (slow)  */

		} else if (Ef_flag) {
			F_hist[isp][hist_hi] = FLHS[isp];
			Eflhs_hist[isp][hist_hi] = Eflhs[isp];
			F_histRHS[isp][hist_hi] = FRHS[isp];
			Efrhs_hist[isp][hist_hi] = Efrhs[isp];
			Tlhs_hist[isp][hist_hi] = Tlhs[isp];
			Trhs_hist[isp][hist_hi] = Trhs[isp];

		} else {
			F_hist[isp][hist_hi] = FLHS[isp];
			F_histRHS[isp][hist_hi] = FRHS[isp];
		}
	}    /* end of for (isp=0...*/


	/***  
	  Parameters not associated with one particular particle species
	 ***/

	if ((T_diag) && (ccc > 1)) {
		Teq[hist_hi] = Tx_hist[0][hist_hi] - Tx_hist[1][hist_hi];
	}

	t_array[hist_hi] = t;
	com_phi_hist[0][hist_hi] = phi_av0;
	com_phi_hist[1][hist_hi] = phi_av05;

	if (slow) {
		com_cur_hist[hist_hi] = exti;
		ese_hist[hist_hi] = 0.5*(rho[0]*phi[0] + rho[nc]*phi[nc]);

		for (i=1; i<nc; i++) {
			ese_hist[hist_hi] += rho[i]*phi[i];
		}

		ese_hist[hist_hi] *= 0.5*area*dx;
		ese_hist[hist_hi] -= epsilon*e[nc]*area*phi[nc];
		ese_hist[hist_hi] /= 1.6022e-19;
		ese_hist[hist_hi] = fabs(ese_hist[hist_hi]+DBL_MIN);

		for(i=0; i<Nnum; i++) {
			Halhs_hist[i][hist_hi] = Halhs[i];
			Harhs_hist[i][hist_hi] = Harhs[i];
		}
	} else if (Ef_flag) {
		for(i=0; i<Nnum; i++) {
			Halhs_hist[i][hist_hi] = Halhs[i];
			Harhs_hist[i][hist_hi] = Harhs[i];
		}
	}
}  /* END of chist() */

/***************************************************************/
/* Rescaling of time histories after 4^n time steps, n=1,2 ... */
/* So that history represents sequence of values Q(i*dt/4^m),  */
/* i= 0, 1,..., m= max(n). */

void newhist() {
	int i, k, isp;

	for (isp=0; isp<nsp; isp++) {
		for (i=1, k=4; i<HISTMAX/4; i++, k+=4) {
			np_hist[isp][i] = np_hist[isp][k];
			F_hist[isp][i] = F_hist[isp][k];
			F_histRHS[isp][i] = F_histRHS[isp][k];
			gridref_l_hist[isp][i] = gridref_l_hist[isp][k];
			gridref_r_hist[isp][i] = gridref_r_hist[isp][k];

			if (slow) {
				kes_hist[isp][i] = kes_hist[isp][k];
				jwall_hist[isp][i]= jwall_hist[isp][k];
				Eflhs_hist[isp][i] = Eflhs_hist[isp][k];
				Efrhs_hist[isp][i] = Efrhs_hist[isp][k];
				Tlhs_hist[isp][i] = Tlhs_hist[isp][k];
				Trhs_hist[isp][i] = Trhs_hist[isp][k];

				if (T_diag) {
					Tx_hist[isp][i] = Tx_hist[isp][k];
					T_dif_hist[isp][i] = T_dif_hist[isp][k];
					Vxt[isp][i] = Vxt[isp][k];
				}
				/*  end of (slow)  */
			} else if (Ef_flag) {
				Eflhs_hist[isp][i] = Eflhs_hist[isp][k];
				Efrhs_hist[isp][i] = Efrhs_hist[isp][k];
				Tlhs_hist[isp][i] = Tlhs_hist[isp][k];
				Trhs_hist[isp][i] = Trhs_hist[isp][k];
			}
		}
	}   /* end of for (isp=0... */

	/***  
	  Parameters not associated with one particular particle species       
	 ***/
	for (i=1, k=4; i<HISTMAX/4; i++, k+=4) {
		com_phi_hist[0][i] = com_phi_hist[0][k];
		com_phi_hist[1][i] = com_phi_hist[1][k];
		com_cur_hist[i] = com_cur_hist[k];
		t_array[i] = t_array[k];

		if (slow) {
			ese_hist[i] = ese_hist[k];

			if ((T_diag) && (ccc > 1)) {
				Teq[i] = Teq[k];
			}

			for(isp=0; isp<Nnum; isp++) {
				Halhs_hist[isp][i] =  Halhs_hist[isp][k];
				Harhs_hist[isp][i] =  Harhs_hist[isp][k];
			}
		} else if (Ef_flag) {
			for(isp=0; isp<Nnum; isp++) {
				Halhs_hist[isp][i] =  Halhs_hist[isp][k];
				Harhs_hist[isp][i] =  Harhs_hist[isp][k];
			}
		}
	}

	hist_hi = HISTMAX/4;
	interval *= 4;

}  /* END of newhist() */

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