#include "peg_pic.h"


/* Calculation, or input from the external files <input file>X.txt,  */
/* X= -, 1, 2, of the neutral particles density and temperature      */
/* and some other parameters for charged-neutral particle collisions */

void neutinp() {
	int    i, j, jj, mm, s, l;
	float  *vn, *vtd, vnt, vtt, trz;
	double *vx, vxt, newdx;
	char   txtname[80];
	FILE   *mytxt;

	ecolsp= ecollisional-1;

	for (j=0; j < Nnum; j++) {
		icolsp[j]--;

		if (n_neut_L[j] <3) {
			for(i = 0; i <= nc; i++) {
				Nneut[j][i] = Nneut0[j];

				if (n_neut_L[j] == 1) {
					/* linear profiles */
					if (i<=nc2) {
						Nneut[j][i] *= (1.-i/nc2);
					} else if (i>=nc3) {
						Nneut[j][i] *= (1. - (nc-i)/(nc-nc3));
					} else {
						Nneut[j][i] = 0.;
					}
 				} else if (n_neut_L[j] == 2)  {
 					/* exponential profiles */
					if (i<nc05) {
						Nneut[j][i] *= exp(-i*dlabda[j]);
					} else if (i>nc05) {
						Nneut[j][i] *= exp((i-nc)*dlabda[j]);
					} else {
						Nneut[j][i] = 0.;
					}
				}

				if (slow) {
					Tneutshow[j][i]= Tneut0[j];
				}

				Tneut[j][i] = sqrt(fabs(Tneut0[j]*1.6022e-19/mn[j]))/vscale;

			} /* end of for(i = 0; i ... */
		} else {
			/* input from the external file */
			if (dshort1) {
				newdx =dx*dshort;
			} else {
				newdx =dx;
			}

			strcpy(txtname, dargv2);
			
			if (j==0) {
				strcat(txtname, ".txt");
			} else if (j==1) {
				strcat(txtname, "1.txt");
			} else {
				strcat(txtname, "2.txt");
			}

			if ((mytxt = fopen(txtname, "r")) == NULL) {
				puts("txt file for neutral density and temperature profiles is missing!");
				exit(1);
			}

			fscanf(mytxt, "%d\n", &mm);
			vx= (double *)malloc((mm+1)*sizeof(double));
			vn= (float *)malloc((mm+1)*sizeof(float));
			vtd= (float *)malloc((mm+1)*sizeof(float));

			for(i=0; i<=mm; i++) {
				fscanf(mytxt, "%lf %e %e\n", &vxt, &vnt, &vtt);

				vx[i] =vxt/newdx;
				vn[i] =vnt;
				vtd[i] =vtt;
				//  printf("x= %lf, n= %e, T= %e\n", vx[i], vn[i], vtd[i]);

				if ((i) && (vx[i]<=vx[i-1])) {
					puts("x_neutral is not monotonic!");
					exit(1);
				}
			}

			if (vx[mm] - xnc > 1.) {
				if (dshort1) {
					printf("vx= %lf\n", vx[i]);
					puts("x_neutral > L_x !");
					exit(1);
				} else {
					i= 0;
					l= 0;
					trz= vx[mm] - dlabda[j];
					
					while (vx[i] < dlabda[j]) {
						i++;
					}
					
					while (vx[l] < trz) {
						l++;
					}

					trz += dlabda[j] - xnc;
					
					for (jj= l+1, s= i+1; jj <= mm; jj++, s++) {
						vx[s] = vx[jj] - trz;
						vn[s] = vn[jj];
						vtd[s] = vtd[jj];
						// printf(" s=%d, jj= %d, %e \n", s, jj, vx[s]);
					}

					mm= s-1;
				}
			}

			jj= vx[0];               /* first point */
			
			for(i=0; i<=jj; i++) {
				Nneut[j][i]=vn[0];
				Tneut[j][i]=vtd[0];
			}  

			jj= vx[mm];             /* last point */
			
			for(i=jj; i<ng; i++) {  
				Nneut[j][i]=vn[mm];
				Tneut[j][i]=vtd[mm];
			}

			for(i=1; i<=mm; i++) {
				vxt= vx[i] - vx[i-1];
				jj = vx[i-1];
				s = vx[i];

				for(l=jj+1; l<=s; l++) {
					vnt= (l - vx[i-1])/vxt;
					vtt= 1. - vnt;
					Nneut[j][l]=vn[i-1]*vtt + vn[i]*vnt;
					Tneut[j][l]=vtd[i-1]*vtt + vtd[i]*vnt;
					//printf("mm=%d, jj=%d, l= %d, x= %d, n= %e, T= %e\n", mm, jj, s, l, Nneut[j][l], Tneut[j][l]);
				}
			}

			for(i=0; i<ng; i++) {
				if (slow) {
					Tneutshow[j][i]= Tneut[j][i];
				}

				Tneut[j][i]=sqrt(fabs(Tneut[j][i]*1.6022e-19/mn[j]))/vscale;
			}

			fclose(mytxt);

		} /* end of else */
	} /* end of "for (j=0, j < Nnum..." */

	mccpar();                                 /* some parameters for mcc */

}  /* End of neutinp */

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