#include "peg_pic.h"

float t_src(int);                     // Time dependent source

/**************** Particle source (first one) ****************/

void psrc() {
	int    j, i, k, isp, R;
	float  vmag, vtempd;
	double maxLen, newLen, midPointPot;
	int realNsp;

#ifdef REFLECTED_OWN_SPECIES
	realNsp = (int) nsp / 2;
#else
	realNsp = nsp;
#endif


	if( sourcD != 4  ) {
		/* current solution.... */
		for (isp=0; isp<realNsp; isp++) {
			/* insert particles for all given particle types... */
			Stemp[isp] += Src_con[isp];

			j = Stemp[isp];

			if (j) {
				/* we have to generate new particles */
				Stemp[isp] -=j;
			
				while (j >= 1) {
					i = np[isp];
					np[isp]++;

					if (sourcD==1) {
						x[isp][i] = nc2 + L_s*frand();
					} else {
						x[isp][i] = nc2 + L_s*(0.5 + asin(2.*frand()-1.)/3.1415926535898);
					}

					/* R = TVMAX*frand();  */
					R = VVMAX*frand();

					if (tas) {
						vtempd = sqrt(fabs(1.-tasr*(x[isp][i]-nc05)/L_s));
						vmag = vtempd*vsperpt[isp];
						vtempd *= vspart[isp];
						maxwell(&vx[isp][i], &vy[isp][i], &vz[isp][i], vtempd, 
								vmag,
								vspar0[isp], vsperp0[isp], R, cs, sn, 0);
					} else {
						maxwell(&vx[isp][i], &vy[isp][i], &vz[isp][i], vspart[isp],
								vsperpt[isp], vspar0[isp], vsperp0[isp], R, cs, sn, 0);
					}

					j --;

				}   /* end of while (j > 1) */
			}    /* end of if (i)... */
		}     /* end of for (isp=0.. */
		
	} else {
		/* new possibility for source Si(x) = exp(alpha * phi(x)) */
		/* first get potential at middle point of the system -> this is the shift value for all potential values */
		if( nc % 2 == 0 ) {
			midPointPot = phi[(nc+1)/2];
		} else {
			midPointPot = phi[nc/2];
		}
	
		for( i=0 ; i < nc ; i++ ){
			phi_six[i].y = phi[i] - midPointPot;
		}
	
		for( i=(int)(0.5+nc2), maxLen=0 ; i < L_s ; i++ ) {
			/* calculate area below phi-x graph - or get the length of all bars under phi-x curve */
			/* this length will be area from which we will pick location for new particles */
			kLen[i] = exp(alpha_Six * phi_six[i].y);
			maxLen += kLen[i];
			totalLen[i] = maxLen;
		}

		for (isp=0; isp<realNsp; isp++) {
			/* insert particles for all given particle types... */
//			printf( "Stemp[%d]=%f --- Src_con[%d]=%f --- MAXLEN=%d\n", isp, Stemp[isp], isp, Src_con[isp], MAXLEN );
			Stemp[isp] += Src_con[isp];

			j = Stemp[isp];

			if (j) {
				/* we have to generate new particles */
				Stemp[isp] -=j;

				while (j >= 1) {
					i = np[isp];
					np[isp]++;

					if( fabs(maxLen) < 1e-5 ) {
						/* potential was 0 - act as sourcD == 1 and evenly distribute new particles in source area */
						x[isp][i] = nc2 + L_s*frand();
						
					} else {
						/* get location for new particle */
						newLen = frand() * maxLen;

						/* find cell in which new particle lies */
						for( k=(int)(0.5+nc2) ; k < L_s && totalLen[k] < newLen ; k++) ;

						if( k == 0 ) {
							x[isp][i] = phi_six[k].x + newLen/kLen[k];
						} else {
							/* position inside cell is defined as ratio between newLen and star length of current cell */
							x[isp][i] = phi_six[k].x + (newLen-totalLen[k-1])/kLen[k];
						}
								
//						printf("totalLen[k]=%g, newLen=%g, totalLen[k-1]=%g, phi_prev_step[k].y=%g --> rez=%g\n", totalLen[k], newLen, totalLen[k-1], phi_prev_step[k].y, (newLen-totalLen[k-1])/exp(alpha_Six*phi_prev_step[k].y));
//						printf( "tstep=%d, j=%d, i=%d, k=%d, phi[k].x=%d, phi[k].y=%g --> x[%d][%d]=%g\n", tstep, j, i, k, phi_prev_step[k].x, phi_prev_step[k].y, isp, i, x[isp][i]);
					}
				
					/* R = TVMAX*frand();  */
					R = VVMAX*frand();

					if (tas) {
						vtempd = sqrt(fabs(1.-tasr*(x[isp][i]-nc05)/L_s));
						vmag = vtempd*vsperpt[isp];
						vtempd *= vspart[isp];
						maxwell(&vx[isp][i], &vy[isp][i], &vz[isp][i], vtempd, 
								vmag,
								vspar0[isp], vsperp0[isp], R, cs, sn, 0);
					} else {
						maxwell(&vx[isp][i], &vy[isp][i], &vz[isp][i], vspart[isp],
								vsperpt[isp], vspar0[isp], vsperp0[isp], R, cs, sn, 0);
					}

					j --;

				}
			}    /* end of if (i)... */
		}     /* end of for (isp=0.. */

//		scanf("Press any key....%d", &k );

	}

}
/* END of psrc() */

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

/**************** Particle source (second one) ****************/

void psrc1() {
	int    j, i, isp, R;
	float  vmag, vtempd;
	int realNsp;

#ifdef REFLECTED_OWN_SPECIES
	realNsp = (int) nsp / 2;
#else
	realNsp = nsp;
#endif

	for (isp=0; isp<realNsp; isp++) {
		Stemp1[isp] += Src_con1[isp]*t_src(isp);
		//printf("%f\n", Stemp1[isp]);
		j = Stemp1[isp];

		if (j) {
			Stemp1[isp] -=j;

			while (j >= 1) {
				i = np[isp];
				np[isp]++;

				if (sourcD1==1) {
					x[isp][i] = nc21 + L_s1*frand();
				} else {
					x[isp][i] = nc21 + L_s1*(0.5 + asin(2.*frand()-1.)/3.1415926535898);
				}

				/* R = TVMAX*frand();  */
				R = VVMAX*frand();

				if (tas1) {
					vtempd = sqrt(fabs(1.-tasr1*(x[isp][i]-nc051)/L_s1));
					vmag = vtempd*vsperpt[isp];
					vtempd *= vspart[isp];
					maxwell(&vx[isp][i], &vy[isp][i], &vz[isp][i], vtempd, vmag,
							vspar01[isp], vsperp01[isp], R, cs, sn, 0);
				} else {
					maxwell(&vx[isp][i], &vy[isp][i], &vz[isp][i], vspart1[isp],
							vsperpt1[isp], vspar01[isp], vsperp01[isp], R, cs, sn, 0);
				}

				j --;

			}   /* end of while (j > 1) */
		}    /* end of if (i)... */
	}     /* end of for (isp=0.. */
}
/* END of psrc1() */

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

/**** Calculation of the time dependent source strength *****/
/**** The relative strength is taken from the input file: ***/
/**** <input file>.src **************************************/

float t_src(int iisp) {
	int		i;
	float		ttmp, htmp;
	static int 	init=1, j, nrow;
	static float 	last;
	FILE   	*mytxt;

	if (timed) {
		if (!iisp) {
			if (init) {
				if ((mytxt = fopen(dargv4, "r")) == NULL) {
					puts("The file for cross-field heat flux is missing!");
					exit(1);
				}

				fscanf(mytxt, "%d\n", &nrow);

				tsrc= (int *)malloc(nrow*sizeof(int));
				hfsrc= (float *)malloc(nrow*sizeof(float));	

				for(i=0; i<nrow; i++) {
					fscanf(mytxt, "%e %e\n", &ttmp, &htmp);
					if (dshort) tsrc[i] = ttmp/dt/dshort;
					else tsrc[i]= ttmp/dt;
					hfsrc[i]= htmp;
				}

				fclose(mytxt);
				j= 0;
				init= 0;
			}

			while((tsrc[j]<tstep) && (j<nrow)) {
				j++;
			}

			if (j==0) {
				last= (tsrc[0]-tstep + hfsrc[0]*tstep)/tsrc[0];   // initial relative strength = 1
			} else if (j>=nrow) {
				last= hfsrc[nrow-1];
			} else {
				last= (hfsrc[j-1]*(tsrc[j]-tstep) + hfsrc[j]*(tstep - tsrc[j-1]))/(tsrc[j]-tsrc[j-1]);
			}
		} // end of if (!iisp) loop

		//	printf("it= %d, j=%d t=%d S=%e\n", tstep, j, tsrc[j], last);
		return(last);
	} else {
		return(1.);
	}

} // End of t_src()

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

/********** Particle heating in the source region ***********/
/********** corresponds to the first source *****************/

void heatup(int iisp, int j, int k) {
	int    R;
	float  vmag, vtempd;

	/* R = TVMAX*frand(); */
	R = VVMAX*frand();

 	if (tas) {
		vtempd = sqrt(fabs(1.-tasr*(x[iisp][j]-nc05)/L_s));
		vmag = vtempd*vsperpt[iisp];
		vtempd *= vspart[iisp];
		maxwell(&vx[iisp][j], &vy[iisp][j], &vz[iisp][j], vtempd, vmag,
				vspar0[iisp], vsperp0[iisp], R, cs, sn, k);
	} else {
		maxwell(&vx[iisp][j], &vy[iisp][j], &vz[iisp][j], vspart[iisp],
				vsperpt[iisp], vspar0[iisp], vsperp0[iisp], R, cs, sn, k);
	}

} /* END of heat() */

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

/********** Particle heating in the source region ***********/
/********** corresponds to the second source *****************/

void heatup1(int iisp, int j, int k) {
	int    R;
	float  vmag, vtempd;

	/* R = TVMAX*frand(); */
	R = VVMAX*frand();

	if (tas1) {
		vtempd = sqrt(fabs(1.-tasr1*(x[iisp][j]-nc051)/L_s1));
		vmag = vtempd*vsperpt[iisp];
		vtempd *= vspart[iisp];
		maxwell(&vx[iisp][j], &vy[iisp][j], &vz[iisp][j], vtempd, vmag,
				vspar01[iisp], vsperp01[iisp], R, cs, sn, k);
	} else {
		maxwell(&vx[iisp][j], &vy[iisp][j], &vz[iisp][j], vspart1[iisp],
				vsperpt1[iisp], vspar01[iisp], vsperp01[iisp], R, cs, sn, k);
	}

} /* END of heat1() */

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


