#include "peg_pic.h"


/**************************************************/
/** Calculation of smoothed n, V, T,... profiles **/ 

void accum() {
	int isp, i;


#pragma omp parallel for private(i)
	for (i=0; i<=nc; i++) {
		rho[i] = rhoback;
	}

	if (slow && dwindow) {
		/* For the X-Vx plot */
#pragma omp parallel for private(isp, i)
		for (isp=0; isp<nsp; isp++) {
			for (i=0; i<np[isp]; i++) {
				xd[isp][i] = x[isp][i];
			}
		}
	}


#pragma omp parallel for private(isp)
	for (isp=0; isp<nsp; isp++) {
		/* Calcul. of  n, V, T, Q profiles */
		avq(isp, qtp);
	}
	
	for (isp=0; isp<nsp; isp++) {
		if (reflux>1) {
			/* the case with reflection, or reinjection */
			One_2_One_d(srho[isp], nc, nsmoothing);
		} else { 
			One_2_One(srho[isp], nc, nsmoothing);
		}

#ifdef INTEL_CC
		if ((qtp==1) || (qtp==2)) {
			if (bmag>0) {
				One_2_One(Vpar[isp], nc, nsmoothing);
			}
			
			if (qtp==1) {
				One_2_One(T_dif[isp], nc, nsmoothing);
			}
			
			One_2_One(Vx[isp], nc, nsmoothing);
			One_2_One(Tx[isp], nc, nsmoothing);
			One_2_One(Qh[isp], nc, nsmoothing);
			One_2_One(Qtot[isp], nc, nsmoothing);
			
		} else if (qtp==3) {
			/* For movie */
			One_2_One(Vpar[isp], nc, nsmoothing);
			One_2_One(Tx[isp], nc, nsmoothing);
			One_2_One(Qh[isp], nc, nsmoothing);
			One_2_One(Qtot[isp], nc, nsmoothing);

		} else if (qtp==4) {
			/* Temperature for effective ionization */
			One_2_One(Tx[isp], nc, nsmoothing);
		}
#endif
#ifdef GCC
		if ((qtp==1) || (qtp==2)) {
	#pragma omp parallel sections
			{
		#pragma omp section 
				{
					if (bmag>0) {
						One_2_One(Vpar[isp], nc, nsmoothing);
					}
				}
				
		#pragma omp section 
				{
					if (qtp==1) {
						One_2_One(T_dif[isp], nc, nsmoothing);
					}
				}
	
		#pragma omp section 
				{
					One_2_One(Vx[isp], nc, nsmoothing);
				}
				
		#pragma omp section 
				{
					One_2_One(Tx[isp], nc, nsmoothing);
				}

		#pragma omp section 
				{
					One_2_One(Qh[isp], nc, nsmoothing);
				}
				
		#pragma omp section 
				{
					One_2_One(Qtot[isp], nc, nsmoothing);
				}
			}
		} else if (qtp==3) {
			/* For movie */
	#pragma omp parallel sections 
			{
		#pragma omp section 
				{
					One_2_One(Vpar[isp], nc, nsmoothing);
				}
	
		#pragma omp section 
				{
					One_2_One(Tx[isp], nc, nsmoothing);
				}
	
		#pragma omp section 
				{
					One_2_One(Qh[isp], nc, nsmoothing);
				}
	
		#pragma omp section 
				{
					One_2_One(Qtot[isp], nc, nsmoothing);
				}
			}
		} else if (qtp==4) {
			/* Temperature for effective ionization */
			One_2_One(Tx[isp], nc, nsmoothing);
		}
#endif
		
#pragma omp parallel for private(i)
		for (i=0; i<=nc; i++) {
			rho[i] += srho[isp][i]*q[isp]*wt[isp]/dx;
		}
	} /* end "for (isp=0, ...  */
} 
/* END of accum() */

/***************************************************************/
/* Smoothing theArray using the 1-2-1 method with the proper   */
/* boundary conditions to conserve charge. */

void One_2_One(float *theArray, int nmax, int loops) {
	int i, j, srcArray;
	float *temp[2], *src, *trg;
	
	temp[0] = (float *)malloc((nmax+1)*sizeof(float));
	temp[1] = (float *)malloc((nmax+1)*sizeof(float));
	
	/* define source and target arrays */
	src = theArray;
	srcArray = 0;
	trg = temp[1-srcArray];
	
	for( j=0 ; j < loops; j++) {
		/* d orequirent numer of loop */
		trg[0]= (3*src[0] + src[1])/4;
	
		for(i=1; i< nmax; i++) {  
			trg[i]= (src[i-1] +6*src[i] +src[i+1])/8;
		}
		
		trg[nmax]= (3*src[nmax] + src[nmax-1])/4;
		
		/* now exchange source and target arrays for next loop */
		srcArray = 1-srcArray;
		src = temp[srcArray];
		trg = temp[1-srcArray];
	}
	
	
	for(i=0; i<= nmax; i++) {
		theArray[i]= src[i];
	}
	
	free(temp[0]);
	free(temp[1]);
}


/***************************************************************/
/* Smoothing theArray using the 1-2-1 method */
/* For infinite plasma */

void One_2_One_d(float *theArray, int nmax, int loops) {
	int i, j, srcArray;
	float *temp[2], *src, *trg;
	
	temp[0] = (float *)malloc((nmax+1)*sizeof(float));
	temp[1] = (float *)malloc((nmax+1)*sizeof(float));
	
	/* define source and target arrays */
	src = theArray;
	srcArray = 0;
	trg = temp[1-srcArray];
	
	for( j=0 ; j < loops; j++) {
		/* d orequirent numer of loop */
		trg[0]= 0.5*(src[0] + src[1]);
	
		for(i=1; i< nmax; i++) {  
			trg[i]= (src[i-1] +2*src[i] +src[i+1])/4;
		}
		
		trg[nmax]= 0.5*(3*src[nmax] + src[nmax-1]);
		
		/* now exchange source and target arrays for next loop */
		srcArray = 1-srcArray;
		src = temp[srcArray];
		trg = temp[1-srcArray];
	}
	
	for(i=0; i<= nmax; i++) {
		theArray[i]= src[i];
	}
	
	free(temp[0]);
	free(temp[1]);
}

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